From owner-freebsd-standards@FreeBSD.ORG Sun Dec 9 20:35:31 2007 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id EFB9216A41A; Sun, 9 Dec 2007 20:35:31 +0000 (UTC) (envelope-from das@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id C81A213C45D; Sun, 9 Dec 2007 20:35:31 +0000 (UTC) (envelope-from das@FreeBSD.org) Received: from freefall.freebsd.org (das@localhost [127.0.0.1]) by freefall.freebsd.org (8.14.2/8.14.2) with ESMTP id lB9KZVuO078079; Sun, 9 Dec 2007 20:35:31 GMT (envelope-from das@freefall.freebsd.org) Received: (from das@localhost) by freefall.freebsd.org (8.14.2/8.14.1/Submit) id lB9KZVkM078075; Sun, 9 Dec 2007 20:35:31 GMT (envelope-from das) Date: Sun, 9 Dec 2007 20:35:31 GMT Message-Id: <200712092035.lB9KZVkM078075@freefall.freebsd.org> To: vincent@vinc17.org, das@FreeBSD.org, freebsd-standards@FreeBSD.org, das@FreeBSD.org From: das@FreeBSD.org Cc: Subject: Re: standards/85080: output of long double subnormals (with printf) is wrong by a factor 2 X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 09 Dec 2007 20:35:32 -0000 Synopsis: output of long double subnormals (with printf) is wrong by a factor 2 State-Changed-From-To: open->patched State-Changed-By: das State-Changed-When: Sun Dec 9 20:34:28 UTC 2007 State-Changed-Why: Fixed in HEAD Responsible-Changed-From-To: freebsd-standards->das Responsible-Changed-By: das Responsible-Changed-When: Sun Dec 9 20:34:28 UTC 2007 Responsible-Changed-Why: Fixed in HEAD http://www.freebsd.org/cgi/query-pr.cgi?pr=85080 From owner-freebsd-standards@FreeBSD.ORG Sun Dec 9 21:25:33 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 3054116A418 for ; Sun, 9 Dec 2007 21:25:33 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id CFBCE13C44B for ; Sun, 9 Dec 2007 21:25:32 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.2/8.14.1) with ESMTP id lB9LP5uu009758; Sun, 9 Dec 2007 16:25:05 -0500 (EST) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.2/8.14.1/Submit) id lB9LP5DB009757; Sun, 9 Dec 2007 16:25:05 -0500 (EST) (envelope-from das@FreeBSD.ORG) Date: Sun, 9 Dec 2007 16:25:05 -0500 From: David Schultz To: Steve Kargl Message-ID: <20071209212505.GA9698@VARK.MIT.EDU> Mail-Followup-To: Steve Kargl , freebsd-standards@FreeBSD.ORG References: <20071103001305.GA82775@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071103001305.GA82775@troutmask.apl.washington.edu> Cc: freebsd-standards@FreeBSD.ORG Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 09 Dec 2007 21:25:33 -0000 On Fri, Nov 02, 2007, Steve Kargl wrote: > With all the success of developing several other missing > C99 math routines, I decided to tackle expl() (which I > need to tackle tanhl()). Hmm, great, but where's the patch? :) Maybe the mailing list software ate it. From owner-freebsd-standards@FreeBSD.ORG Sun Dec 9 21:39:17 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 2287416A46E for ; Sun, 9 Dec 2007 21:39:17 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.freebsd.org (Postfix) with ESMTP id EABE713C4D1 for ; Sun, 9 Dec 2007 21:39:16 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.1/8.14.1) with ESMTP id lB9LYoZ4096090 for ; Sun, 9 Dec 2007 13:34:50 -0800 (PST) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.1/8.14.1/Submit) id lB9LYo1F096089 for freebsd-standards@FreeBSD.ORG; Sun, 9 Dec 2007 13:34:50 -0800 (PST) (envelope-from sgk) Date: Sun, 9 Dec 2007 13:34:50 -0800 From: Steve Kargl To: freebsd-standards@FreeBSD.ORG Message-ID: <20071209213450.GA95726@troutmask.apl.washington.edu> References: <20071103001305.GA82775@troutmask.apl.washington.edu> <20071209212505.GA9698@VARK.MIT.EDU> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071209212505.GA9698@VARK.MIT.EDU> User-Agent: Mutt/1.4.2.3i Cc: Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 09 Dec 2007 21:39:17 -0000 On Sun, Dec 09, 2007 at 04:25:05PM -0500, David Schultz wrote: > On Fri, Nov 02, 2007, Steve Kargl wrote: > > With all the success of developing several other missing > > C99 math routines, I decided to tackle expl() (which I > > need to tackle tanhl()). > > Hmm, great, but where's the patch? :) Maybe the mailing list > software ate it. This is the current version. I need to revise how I computed the ploynomial coefficient. In short, I mapped r in [-ln(2)/2:ln(2)/2] into the range x in [-1,1] for the Chebyshev interpolation, but I never scaled x back into r. This is the reason why the lines "r = r * TOLN2;" exists. I don't remember if bde sent me comments on this code. I sure he has plenty. :) steve /*- * Copyright (c) 2007 Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Use argument reduction to compute exp(x). The reduction writes * x = k * ln(2) + r with r in the range [0, ln(2)]. This then * gives exp(x) = 2**k * exp(r), and exp(r) is evaluated via a nearly * minimax polynomial approximation such that r is mapped into [-1:1]. */ #include "math.h" #include "math_private.h" #include "fpmath.h" /* ln(LDBL_MAX) = 11356.523406294144 */ #define XMAX 0x2.c5c85fdf473de6ap12L /* ln(LDBL_MIN) = -11355.137111933024 */ #define XMIN -0x2.c5b2319c4843accp12L /* | ln(smallest subnormal) | = 11399.498531488861 */ #define GRAD 0x2.c877f9fc278aeaap12L #define LN2 0xb.17217f7d1cf79acp-4L /* ln(2) */ #define LN2HI 0xb.17217f7d2000000p-4L #define LN2LO -0x3.08654361c4c67fcp-44L #define TOLN2 0x2.e2a8eca5705fc2fp0L /* 2/ln(2) */ #define ILN2 0x1.71547652b82fe17p0L /* 1/ln(2) */ #define ILN2HI 0x1.71547652b800000p0L #define ILN2LO 0x2.fe1777d0ffda0d2p-44L #define LN2O2 0x5.8b90bfbe8e7bcd6p-4L /* ln(2)/2 */ #define ZERO 0.L /* * This set of coefficients is used in the polynomial approximation * for exp(r) where r is in [-ln(2)/2:ln(2)/2], and the r comes from * the argument reduction of x. */ #define C0 0x1.000000000000000p0L #define C1 0x5.8b90bfbe8e7bcd6p-4L #define C2 0xf.5fdeffc162c7543p-8L #define C3 0x1.c6b08d704a0bf8bp-8L #define C4 0x2.76556df749cee54p-12L #define C5 0x2.bb0ffcf14ce6221p-16L #define C6 0x2.861225f0d8f0edfp-20L #define C7 0x1.ffcbfc588b0c687p-24L #define C8 0x1.62c0223a5c823fdp-28L #define C9 0xd.a929e9caf3e1ed2p-36L #define C10 0x7.933d4562e3b2cd7p-40L #define C11 0x3.d1958e6a3764b64p-44L #define C12 0x1.c3bd650fc1e343ap-48L #define C13 0xc.0b0c98b3649ff26p-56L #define C14 0x4.c525936609b02cfp-60L #define C15 0x1.c36e84400493e74p-64L long double expl(long double x) { union IEEEl2bits z; int k, s; long double r; z.e = x; s = z.bits.sign; z.bits.sign = 0; /* x is either 0 or a subnormal number. */ if (z.bits.exp == 0) { if ((z.bits.manl | z.bits.manh) == 0) return (1); else return (1 + x); } if (XMIN <= x && x <= XMAX) { /* Argument reduction. */ k = (int) (z.e * ILN2HI + z.e * ILN2LO); r = z.e - k * LN2HI - k * LN2LO; if (r > LN2O2) { r -= LN2; k++; } /* Compute exp(r) via the polynomial approximation. */ r = r * TOLN2; z.e = C0 + r * (C1 + r * (C2 + r * (C3 + r * (C4 + r * (C5 + r * (C6 + r * (C7 + r * (C8 + r * (C9 + r * (C10 + r * (C11 + r * (C12 + r * (C13 + r * (C14 + r * C15)))))))))))))); if (s) { z.e = 1 / z.e; z.bits.exp -= k; } else z.bits.exp += k; return (z.e); } /* * If x = +Inf, then exp(x) = Inf. * If x = -Inf, then exp(x) = 0. * If x = NaN, then exp(x) = NaN. */ if (z.bits.exp == 32767) { mask_nbit_l(z); if (!(z.bits.manh | z.bits.manl)) return (s ? ZERO : 1.L / ZERO); return ((x - x) / (x - x)); } /* If x > 0 then, overflow to +Inf. */ if (!s) return (1.L / ZERO); /* For x < 0, check if gradual underflow is needed. */ if (z.e > GRAD) return (ZERO); /* Argument reduction. */ k = (int) (z.e * ILN2); r = z.e - k * LN2HI - k * LN2LO; if (r > LN2O2) { r -= LN2; k++; } /* Compute exp(r) via the polynomial approximation. */ r *= TOLN2; z.e = C0 + r * (C1 + r * (C2 + r * (C3 + r * (C4 + r * (C5 + r * (C6 + r * (C7 + r * (C8 + r * (C9 + r * (C10 + r * (C11 + r * (C12 + r * (C13 + r * (C14 + r * C15)))))))))))))); z.e = 1 / z.e; /* * FIXME: There has to be a better way to handle gradual underflow * because the relative absolute error is fairly large for numerous * vlaues of x as exp(x) goes to zero. */ z.e = scalbnl(z.e, -k); return (z.e); } From owner-freebsd-standards@FreeBSD.ORG Sun Dec 9 21:49:55 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 7CBD216A417 for ; Sun, 9 Dec 2007 21:49:55 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id 3D69F13C467 for ; Sun, 9 Dec 2007 21:49:55 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.2/8.14.1) with ESMTP id lB9LnRf6009881; Sun, 9 Dec 2007 16:49:27 -0500 (EST) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.2/8.14.1/Submit) id lB9LnRXH009880; Sun, 9 Dec 2007 16:49:27 -0500 (EST) (envelope-from das@FreeBSD.ORG) Date: Sun, 9 Dec 2007 16:49:27 -0500 From: David Schultz To: Bruce Evans Message-ID: <20071209214927.GB9698@VARK.MIT.EDU> Mail-Followup-To: Bruce Evans , Steve Kargl , freebsd-standards@FreeBSD.ORG References: <20071012180959.GA36345@troutmask.apl.washington.edu> <20071206090833.GA95428@VARK.MIT.EDU> <20071206231143.GA63969@troutmask.apl.washington.edu> <20071207173222.D702@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071207173222.D702@delplex.bde.org> Cc: freebsd-standards@FreeBSD.ORG, Steve Kargl Subject: Re: [PATCH] hypotl, cabsl, and code removal in cabs X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 09 Dec 2007 21:49:55 -0000 On Fri, Dec 07, 2007, Bruce Evans wrote: > >As to the sqrtl question, I have an implementation that supposely > >does correct rounding in all rounding modes. It is restricted to > >64-bit significand long doubles. The code does not use bit twiddle; > >instead, it uses fenv. > > This I haven't looked at closely. I fear extreme slowness. [...] > Anyway, the software version of sqrtl is irrelevant on > athlon-xp, since athlon-xp has sqrtl in hardware (takes 35 cycles). > Similarly for amd64, ia64 and possibly sparc64 (sparc64 has sqrt in > hardware so it hopefully has sqrtl in hardware). arm and powerpc > apparently have long double == double, so the software version of sqrtl > is apparently only needed on ia64. In general, even if we don't use it on any architectures that FreeBSD currently supports, I'd like to have a working MI implementation checked in before we add a bunch of MD versions. Otherwise we risk overburdening people trying to port FreeBSD to new architectures. Having an MI version makes cross-testing easier, too. In the specific case of sqrt, pretty much everyone who claims to support IEEE 754 floating point is going to have hardware support for it, so maybe we can get away with a bunch of MD versions alone if we have to... For ia64, we can use Intel's BSD-licensed math lib. From owner-freebsd-standards@FreeBSD.ORG Sun Dec 9 22:39:45 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id E507916A418 for ; Sun, 9 Dec 2007 22:39:45 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id A8E8613C468 for ; Sun, 9 Dec 2007 22:39:45 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.2/8.14.1) with ESMTP id lB9MdI17010110; Sun, 9 Dec 2007 17:39:18 -0500 (EST) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.2/8.14.1/Submit) id lB9MdIxx010109; Sun, 9 Dec 2007 17:39:18 -0500 (EST) (envelope-from das@FreeBSD.ORG) Date: Sun, 9 Dec 2007 17:39:18 -0500 From: David Schultz To: Steve Kargl Message-ID: <20071209223918.GA9920@VARK.MIT.EDU> Mail-Followup-To: Steve Kargl , freebsd-standards@FreeBSD.ORG References: <20071103001305.GA82775@troutmask.apl.washington.edu> <20071209212505.GA9698@VARK.MIT.EDU> <20071209213450.GA95726@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071209213450.GA95726@troutmask.apl.washington.edu> Cc: freebsd-standards@FreeBSD.ORG Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 09 Dec 2007 22:39:46 -0000 On Sun, Dec 09, 2007, Steve Kargl wrote: > On Sun, Dec 09, 2007 at 04:25:05PM -0500, David Schultz wrote: > > On Fri, Nov 02, 2007, Steve Kargl wrote: > > > With all the success of developing several other missing > > > C99 math routines, I decided to tackle expl() (which I > > > need to tackle tanhl()). > > > > Hmm, great, but where's the patch? :) Maybe the mailing list > > software ate it. > > This is the current version. I need to revise how I computed > the ploynomial coefficient. In short, I mapped r in > [-ln(2)/2:ln(2)/2] into the range x in [-1,1] for the Chebyshev > interpolation, but I never scaled x back into r. This is the > reason why the lines "r = r * TOLN2;" exists. Some minor fixes (mostly whitespace) are below. Also, don't you lose precision when you do stuff like this? z.e = 1 / z.e; In any case, if you can get me the appropriate constants for the 128-bit format, I'll clean everything up and check it in, which will make a lot of ports maintainers happy. That will pave the way to other stuff (e.g., MD versions of this) as well. We can worry about subnormals later. From owner-freebsd-standards@FreeBSD.ORG Sun Dec 9 22:42:01 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 798CC16A420 for ; Sun, 9 Dec 2007 22:42:01 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id 4E82413C455 for ; Sun, 9 Dec 2007 22:42:01 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.2/8.14.1) with ESMTP id lB9MfXv2010168; Sun, 9 Dec 2007 17:41:33 -0500 (EST) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.2/8.14.1/Submit) id lB9MfXKH010167; Sun, 9 Dec 2007 17:41:33 -0500 (EST) (envelope-from das@FreeBSD.ORG) Date: Sun, 9 Dec 2007 17:41:33 -0500 From: David Schultz To: Steve Kargl , freebsd-standards@FreeBSD.ORG Message-ID: <20071209224133.GA10128@VARK.MIT.EDU> Mail-Followup-To: Steve Kargl , freebsd-standards@FreeBSD.ORG References: <20071103001305.GA82775@troutmask.apl.washington.edu> <20071209212505.GA9698@VARK.MIT.EDU> <20071209213450.GA95726@troutmask.apl.washington.edu> <20071209223918.GA9920@VARK.MIT.EDU> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071209223918.GA9920@VARK.MIT.EDU> Cc: Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 09 Dec 2007 22:42:01 -0000 On Sun, Dec 09, 2007, David Schultz wrote: > On Sun, Dec 09, 2007, Steve Kargl wrote: > > On Sun, Dec 09, 2007 at 04:25:05PM -0500, David Schultz wrote: > > > On Fri, Nov 02, 2007, Steve Kargl wrote: > > > > With all the success of developing several other missing > > > > C99 math routines, I decided to tackle expl() (which I > > > > need to tackle tanhl()). > > > > > > Hmm, great, but where's the patch? :) Maybe the mailing list > > > software ate it. > > > > This is the current version. I need to revise how I computed > > the ploynomial coefficient. In short, I mapped r in > > [-ln(2)/2:ln(2)/2] into the range x in [-1,1] for the Chebyshev > > interpolation, but I never scaled x back into r. This is the > > reason why the lines "r = r * TOLN2;" exists. > > Some minor fixes (mostly whitespace) are below. Oops, here it is. --- /usr/src/lib/msun/src/e_expl.c.bak 2007-12-09 17:30:35.000000000 -0500 +++ /usr/src/lib/msun/src/e_expl.c 2007-12-09 17:34:14.000000000 -0500 @@ -38,45 +38,46 @@ #define XMAX 0x2.c5c85fdf473de6ap12L /* ln(LDBL_MIN) = -11355.137111933024 */ -#define XMIN -0x2.c5b2319c4843accp12L +#define XMIN -0x2.c5b2319c4843accp12L /* | ln(smallest subnormal) | = 11399.498531488861 */ #define GRAD 0x2.c877f9fc278aeaap12L -#define LN2 0xb.17217f7d1cf79acp-4L /* ln(2) */ -#define LN2HI 0xb.17217f7d2000000p-4L -#define LN2LO -0x3.08654361c4c67fcp-44L +#define LN2 0xb.17217f7d1cf79acp-4L /* ln(2) */ +#define LN2HI 0xb.17217f7d2000000p-4L +#define LN2LO -0x3.08654361c4c67fcp-44L -#define TOLN2 0x2.e2a8eca5705fc2fp0L /* 2/ln(2) */ +#define TOLN2 0x2.e2a8eca5705fc2fp0L /* 2/ln(2) */ -#define ILN2 0x1.71547652b82fe17p0L /* 1/ln(2) */ -#define ILN2HI 0x1.71547652b800000p0L -#define ILN2LO 0x2.fe1777d0ffda0d2p-44L +#define ILN2 0x1.71547652b82fe17p0L /* 1/ln(2) */ +#define ILN2HI 0x1.71547652b800000p0L +#define ILN2LO 0x2.fe1777d0ffda0d2p-44L -#define LN2O2 0x5.8b90bfbe8e7bcd6p-4L /* ln(2)/2 */ +#define LN2O2 0x5.8b90bfbe8e7bcd6p-4L /* ln(2)/2 */ + +#define ZERO 0.L -#define ZERO 0.L /* * This set of coefficients is used in the polynomial approximation * for exp(r) where r is in [-ln(2)/2:ln(2)/2], and the r comes from * the argument reduction of x. */ -#define C0 0x1.000000000000000p0L -#define C1 0x5.8b90bfbe8e7bcd6p-4L -#define C2 0xf.5fdeffc162c7543p-8L -#define C3 0x1.c6b08d704a0bf8bp-8L -#define C4 0x2.76556df749cee54p-12L -#define C5 0x2.bb0ffcf14ce6221p-16L -#define C6 0x2.861225f0d8f0edfp-20L -#define C7 0x1.ffcbfc588b0c687p-24L -#define C8 0x1.62c0223a5c823fdp-28L -#define C9 0xd.a929e9caf3e1ed2p-36L -#define C10 0x7.933d4562e3b2cd7p-40L -#define C11 0x3.d1958e6a3764b64p-44L -#define C12 0x1.c3bd650fc1e343ap-48L -#define C13 0xc.0b0c98b3649ff26p-56L -#define C14 0x4.c525936609b02cfp-60L -#define C15 0x1.c36e84400493e74p-64L +#define C0 0x1.000000000000000p0L +#define C1 0x5.8b90bfbe8e7bcd6p-4L +#define C2 0xf.5fdeffc162c7543p-8L +#define C3 0x1.c6b08d704a0bf8bp-8L +#define C4 0x2.76556df749cee54p-12L +#define C5 0x2.bb0ffcf14ce6221p-16L +#define C6 0x2.861225f0d8f0edfp-20L +#define C7 0x1.ffcbfc588b0c687p-24L +#define C8 0x1.62c0223a5c823fdp-28L +#define C9 0xd.a929e9caf3e1ed2p-36L +#define C10 0x7.933d4562e3b2cd7p-40L +#define C11 0x3.d1958e6a3764b64p-44L +#define C12 0x1.c3bd650fc1e343ap-48L +#define C13 0xc.0b0c98b3649ff26p-56L +#define C14 0x4.c525936609b02cfp-60L +#define C15 0x1.c36e84400493e74p-64L long double expl(long double x) @@ -90,15 +91,10 @@ z.bits.sign = 0; /* x is either 0 or a subnormal number. */ - if (z.bits.exp == 0) { - if ((z.bits.manl | z.bits.manh) == 0) - return (1); - else - return (1 + x); - } + if (z.bits.exp == 0) + return (1 + x); if (XMIN <= x && x <= XMAX) { - /* Argument reduction. */ k = (int) (z.e * ILN2HI + z.e * ILN2LO); r = z.e - k * LN2HI - k * LN2LO; @@ -117,11 +113,13 @@ if (s) { z.e = 1 / z.e; z.bits.exp -= k; - } else + } else { z.bits.exp += k; + } return (z.e); } + /* * If x = +Inf, then exp(x) = Inf. * If x = -Inf, then exp(x) = 0. @@ -157,10 +155,11 @@ (C6 + r * (C7 + r * (C8 + r * (C9 + r * (C10 + r * (C11 + r * (C12 + r * (C13 + r * (C14 + r * C15)))))))))))))); z.e = 1 / z.e; + /* * FIXME: There has to be a better way to handle gradual underflow * because the relative absolute error is fairly large for numerous - * vlaues of x as exp(x) goes to zero. + * values of x as exp(x) goes to zero. */ z.e = scalbnl(z.e, -k); From owner-freebsd-standards@FreeBSD.ORG Mon Dec 10 05:18:51 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 7330216A419; Mon, 10 Dec 2007 05:18:51 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail01.syd.optusnet.com.au (mail01.syd.optusnet.com.au [211.29.132.182]) by mx1.freebsd.org (Postfix) with ESMTP id E915013C458; Mon, 10 Dec 2007 05:18:50 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from c211-30-219-213.carlnfd3.nsw.optusnet.com.au (c211-30-219-213.carlnfd3.nsw.optusnet.com.au [211.30.219.213]) by mail01.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id lBA5IdoH025443 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Mon, 10 Dec 2007 16:18:43 +1100 Date: Mon, 10 Dec 2007 16:18:39 +1100 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: David Schultz In-Reply-To: <20071209223918.GA9920@VARK.MIT.EDU> Message-ID: <20071210153413.W3836@delplex.bde.org> References: <20071103001305.GA82775@troutmask.apl.washington.edu> <20071209212505.GA9698@VARK.MIT.EDU> <20071209213450.GA95726@troutmask.apl.washington.edu> <20071209223918.GA9920@VARK.MIT.EDU> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@FreeBSD.org, Steve Kargl Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Dec 2007 05:18:51 -0000 On Sun, 9 Dec 2007, David Schultz wrote: > Some minor fixes (mostly whitespace) are below. > > Also, don't you lose precision when you do stuff like this? > z.e = 1 / z.e; > > In any case, if you can get me the appropriate constants for the > 128-bit format, I'll clean everything up and check it in, which > will make a lot of ports maintainers happy. That will pave the way > to other stuff (e.g., MD versions of this) as well. We can worry > about subnormals later. Why not convert the fdlibm algorithm for exp() as is done for expf()? It is much better (*), doesn't need to be debugged (except for the conversion), and would be easier to maintain. Better algorithms exist, like someone named das@ used for exp2(), but would be harder to debug and maintain. amd64 should probably use hardware for expl(), since long doubles are much harder to handle efficiently than doubles. Thus the software algorithm is needed mainly for ia64 and sparc64. (*) Starting with the transformation of exp(x) to x*(exp(x)+1)/(exp(x)-1). The power series for the latter converges much faster and has better numeric properties. Polynomial approximations inherit these properties. Bruce From owner-freebsd-standards@FreeBSD.ORG Mon Dec 10 05:29:50 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 7340516A41A; Mon, 10 Dec 2007 05:29:50 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail05.syd.optusnet.com.au (mail05.syd.optusnet.com.au [211.29.132.186]) by mx1.freebsd.org (Postfix) with ESMTP id E9A1213C447; Mon, 10 Dec 2007 05:29:49 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from c211-30-219-213.carlnfd3.nsw.optusnet.com.au (c211-30-219-213.carlnfd3.nsw.optusnet.com.au [211.30.219.213]) by mail05.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id lBA5TjlP011340 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Mon, 10 Dec 2007 16:29:47 +1100 Date: Mon, 10 Dec 2007 16:29:45 +1100 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: David Schultz In-Reply-To: <20071209214927.GB9698@VARK.MIT.EDU> Message-ID: <20071210162247.O3964@delplex.bde.org> References: <20071012180959.GA36345@troutmask.apl.washington.edu> <20071206090833.GA95428@VARK.MIT.EDU> <20071206231143.GA63969@troutmask.apl.washington.edu> <20071207173222.D702@delplex.bde.org> <20071209214927.GB9698@VARK.MIT.EDU> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@FreeBSD.org, Steve Kargl Subject: Re: [PATCH] hypotl, cabsl, and code removal in cabs X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Dec 2007 05:29:50 -0000 On Sun, 9 Dec 2007, David Schultz wrote: > On Fri, Dec 07, 2007, Bruce Evans wrote: > [...] >> Anyway, the software version of sqrtl is irrelevant on >> athlon-xp, since athlon-xp has sqrtl in hardware (takes 35 cycles). >> Similarly for amd64, ia64 and possibly sparc64 (sparc64 has sqrt in >> hardware so it hopefully has sqrtl in hardware). arm and powerpc >> apparently have long double == double, so the software version of sqrtl >> is apparently only needed on ia64. > > In general, even if we don't use it on any architectures that > FreeBSD currently supports, I'd like to have a working MI > implementation checked in before we add a bunch of MD > versions. Otherwise we risk overburdening people trying to port > FreeBSD to new architectures. Having an MI version makes > cross-testing easier, too. I agree. I could write amd64 and i386 versions in asm the most important functions in about 5 minutes each, but have refrained from doing so since I think this would get in the way of general support. > For ia64, we can use Intel's BSD-licensed math lib. Sure. I use it a hacked up version of the glibc version of it locally as a gold standard, but don't want to maintain it. Is is fully BSD-licensed? Bruce From owner-freebsd-standards@FreeBSD.ORG Mon Dec 10 07:17:20 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 4945C16A468 for ; Mon, 10 Dec 2007 07:17:20 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id 0957A13C4E1 for ; Mon, 10 Dec 2007 07:17:19 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.2/8.14.1) with ESMTP id lBA7GoaX012987; Mon, 10 Dec 2007 02:16:50 -0500 (EST) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.2/8.14.1/Submit) id lBA7GoHT012986; Mon, 10 Dec 2007 02:16:50 -0500 (EST) (envelope-from das@FreeBSD.ORG) Date: Mon, 10 Dec 2007 02:16:50 -0500 From: David Schultz To: Bruce Evans Message-ID: <20071210071650.GA12886@VARK.MIT.EDU> Mail-Followup-To: Bruce Evans , Steve Kargl , freebsd-standards@FreeBSD.ORG References: <20071103001305.GA82775@troutmask.apl.washington.edu> <20071209212505.GA9698@VARK.MIT.EDU> <20071209213450.GA95726@troutmask.apl.washington.edu> <20071209223918.GA9920@VARK.MIT.EDU> <20071210153413.W3836@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071210153413.W3836@delplex.bde.org> Cc: freebsd-standards@FreeBSD.ORG, Steve Kargl Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Dec 2007 07:17:20 -0000 On Mon, Dec 10, 2007, Bruce Evans wrote: > On Sun, 9 Dec 2007, David Schultz wrote: > > >Some minor fixes (mostly whitespace) are below. > > > >Also, don't you lose precision when you do stuff like this? > > z.e = 1 / z.e; > > > >In any case, if you can get me the appropriate constants for the > >128-bit format, I'll clean everything up and check it in, which > >will make a lot of ports maintainers happy. That will pave the way > >to other stuff (e.g., MD versions of this) as well. We can worry > >about subnormals later. > > Why not convert the fdlibm algorithm for exp() as is done for expf()? > It is much better (*), doesn't need to be debugged (except for the > conversion), and would be easier to maintain. > > Better algorithms exist, like someone named das@ used for exp2(), but > would be harder to debug and maintain. I'm not worried about maintenance, since I don't expect God to add any major new features to e^x any time soon. Writing exp2() took a lot of reading papers and tinkering, and it's a pain in the neck to generate the constants and figure out the resulting error for each interval. I seem to recall trying the same tricks for expl(), but there were problems with rescaling in base e instead of base 2 without losing accuracy. In any case, I don't have the kind of time needed to fix all of that stuff now. I don't really care how expl() is implemented; anything that works is better than anything that doesn't. If someone comes up with a better, cleverer scheme later, that's great, but we've been talking about a lot of this stuff forever and there's still nothing in the tree. From owner-freebsd-standards@FreeBSD.ORG Mon Dec 10 07:17:55 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 32F1816A417 for ; Mon, 10 Dec 2007 07:17:55 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id E772D13C469 for ; Mon, 10 Dec 2007 07:17:54 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.2/8.14.1) with ESMTP id lBA7HQgF013000; Mon, 10 Dec 2007 02:17:26 -0500 (EST) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.2/8.14.1/Submit) id lBA7HQ9V012999; Mon, 10 Dec 2007 02:17:26 -0500 (EST) (envelope-from das@FreeBSD.ORG) Date: Mon, 10 Dec 2007 02:17:26 -0500 From: David Schultz To: Bruce Evans Message-ID: <20071210071726.GB12886@VARK.MIT.EDU> Mail-Followup-To: Bruce Evans , Steve Kargl , freebsd-standards@FreeBSD.ORG References: <20071012180959.GA36345@troutmask.apl.washington.edu> <20071206090833.GA95428@VARK.MIT.EDU> <20071206231143.GA63969@troutmask.apl.washington.edu> <20071207173222.D702@delplex.bde.org> <20071209214927.GB9698@VARK.MIT.EDU> <20071210162247.O3964@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071210162247.O3964@delplex.bde.org> Cc: freebsd-standards@FreeBSD.ORG, Steve Kargl Subject: Re: [PATCH] hypotl, cabsl, and code removal in cabs X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Dec 2007 07:17:55 -0000 On Mon, Dec 10, 2007, Bruce Evans wrote: > >For ia64, we can use Intel's BSD-licensed math lib. > > Sure. I use it a hacked up version of the glibc version of it locally > as a gold standard, but don't want to maintain it. Is is fully > BSD-licensed? 3-clause, but yes, that's my understanding. From owner-freebsd-standards@FreeBSD.ORG Mon Dec 10 10:16:33 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 206C216A417; Mon, 10 Dec 2007 10:16:33 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail07.syd.optusnet.com.au (mail07.syd.optusnet.com.au [211.29.132.188]) by mx1.freebsd.org (Postfix) with ESMTP id B12F813C4D1; Mon, 10 Dec 2007 10:16:32 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from besplex.bde.org (c211-30-219-213.carlnfd3.nsw.optusnet.com.au [211.30.219.213]) by mail07.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id lBAAGRli004629 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Mon, 10 Dec 2007 21:16:30 +1100 Date: Mon, 10 Dec 2007 21:16:27 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: David Schultz In-Reply-To: <20071210071650.GA12886@VARK.MIT.EDU> Message-ID: <20071210200058.Y30409@besplex.bde.org> References: <20071103001305.GA82775@troutmask.apl.washington.edu> <20071209212505.GA9698@VARK.MIT.EDU> <20071209213450.GA95726@troutmask.apl.washington.edu> <20071209223918.GA9920@VARK.MIT.EDU> <20071210153413.W3836@delplex.bde.org> <20071210071650.GA12886@VARK.MIT.EDU> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@FreeBSD.ORG, Steve Kargl Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Dec 2007 10:16:33 -0000 On Mon, 10 Dec 2007, David Schultz wrote: > On Mon, Dec 10, 2007, Bruce Evans wrote: >> Better algorithms exist, like someone named das@ used for exp2(), but >> would be harder to debug and maintain. > > I'm not worried about maintenance, since I don't expect God to add > any major new features to e^x any time soon. But every different precision requires a significantly different implementation (mainly for magic number, but table-driven methods give a huge number of those) since we don't have algorithms or tools to generate all the magic numbers from the precision parameter. > Writing exp2() took a > lot of reading papers and tinkering, and it's a pain in the neck > to generate the constants and figure out the resulting error for > each interval. I seem to recall trying the same tricks for expl(), > but there were problems with rescaling in base e instead of base 2 > without losing accuracy. In any case, I don't have the kind of time > needed to fix all of that stuff now. The problems are well understood :-). I didn't understand them when you started working on exp2() but do now. Essentially, they are due to the strange phenomenon that linear functions are usually harder to approximate (to within 1 ULP) than transcendental functions expanded in a Taylor series about 0. The linear scaling step with a not-exactly representable factor like log2(e) ends up being the most difficult point (unless it is done inaccurately). > I don't really care how expl() is implemented; anything that works > is better than anything that doesn't. If someone comes up with a > better, cleverer scheme later, that's great, but we've been > talking about a lot of this stuff forever and there's still > nothing in the tree. I only agree if "working" includes no new bugs and good accuracy at little cost. We still don't support 64-bit precision on i386. Fixing that has higher priority. Bruce From owner-freebsd-standards@FreeBSD.ORG Mon Dec 10 11:07:13 2007 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 5741E16A47D for ; Mon, 10 Dec 2007 11:07:13 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 4601D13C45A for ; Mon, 10 Dec 2007 11:07:13 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.2/8.14.2) with ESMTP id lBAB7D5C073503 for ; Mon, 10 Dec 2007 11:07:13 GMT (envelope-from owner-bugmaster@FreeBSD.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.2/8.14.1/Submit) id lBAB7ClV073499 for freebsd-standards@FreeBSD.org; Mon, 10 Dec 2007 11:07:12 GMT (envelope-from owner-bugmaster@FreeBSD.org) Date: Mon, 10 Dec 2007 11:07:12 GMT Message-Id: <200712101107.lBAB7ClV073499@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: gnats set sender to owner-bugmaster@FreeBSD.org using -f From: FreeBSD bugmaster To: freebsd-standards@FreeBSD.org Cc: Subject: Current problem reports assigned to freebsd-standards@FreeBSD.org X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Dec 2007 11:07:13 -0000 Current FreeBSD problem reports Critical problems Serious problems S Tracker Resp. Description -------------------------------------------------------------------------------- o bin/25542 standards /bin/sh: null char in quoted string o stand/54410 standards one-true-awk not POSIX compliant (no extended REs) o stand/82654 standards C99 long double math functions are missing o stand/94729 standards [libc] fcntl() throws undocumented ENOTTY 4 problems total. Non-critical problems S Tracker Resp. Description -------------------------------------------------------------------------------- o bin/21519 standards sys/dir.h should be deprecated some more o bin/24390 standards Replacing old dir-symlinks when using /bin/ln s stand/24590 standards timezone function not compatible witn Single Unix Spec s stand/36076 standards Implementation of POSIX fuser command o stand/39256 standards snprintf/vsnprintf aren't POSIX-conformant for strings p stand/41576 standards POSIX compliance of ln(1) o stand/44425 standards getcwd() succeeds even if current dir has perm 000. o stand/46119 standards Priority problems for SCHED_OTHER using pthreads o stand/54833 standards [pcvt] more pcvt deficits o stand/54839 standards [pcvt] pcvt deficits p stand/55112 standards glob.h, glob_t's gl_pathc should be "size_t", not "int o stand/56476 standards cd9660 unicode support simple hack o stand/58676 standards grantpt(3) alters storage used by ptsname(3) s stand/62858 standards malloc(0) not C99 compliant s kern/64875 standards [libc] [patch] [feature request] add a system call: fd o stand/66357 standards make POSIX conformance problem ('sh -e' & '+' command- o stand/66531 standards _gettemp uses a far smaller set of filenames than docu o stand/70813 standards [PATCH] ls(1) not Posix compliant o stand/72006 standards floating point formating in non-C locales o stand/79056 standards regex(3) regression tests a stand/80293 standards sysconf() does not support well-defined unistd values o stand/81287 standards [PATCH]: fingerd(8) might send a line not ending in CR o stand/83845 standards [libm] [patch] add log2() and log2f() support for libm o stand/92360 standards [headers] [patch] Missing TAB3 in kernel headers o stand/92362 standards [headers] [patch] Missing SIGPOLL in kernel headers o kern/93705 standards [headers] [patch] ENODATA and EGREGIOUS (for glibc com o stand/96016 standards [headers] clock_getres et al should be in o stand/96236 standards [PATCH] [POSIX] sed.1 incorrectly describes a function p stand/99517 standards Missing SIGRTMIN and SIGRTMAX signals o stand/99960 standards [Patch] make(1): Add -p flag o stand/100017 standards [Patch] Add fuser(1) functionality to fstat(1) o stand/104743 standards [headers] [patch] Wrong values for _POSIX_ minimal lim o stand/104841 standards [libm] [patch] C99 long double square root. o stand/107561 standards [patch] Missing SUS function tcgetsid() o kern/114578 standards [libc] wide character printing using swprintf(dst, n, o stand/114633 standards /etc/rc.subr: line 511: omits a quotation mark: "force o stand/116081 standards make does not work with the directive sinclude o stand/116221 standards SUS issue -- FreeBSD has not flag WNOWAIT for wait*() o stand/116826 standards [PATCH] sh support for POSIX character classes o stand/118047 standards SUGGESTION: /etc/printcap vs mergemaster 40 problems total. From owner-freebsd-standards@FreeBSD.ORG Mon Dec 10 11:22:17 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id BFC5E16A49E for ; Mon, 10 Dec 2007 11:22:17 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail01.syd.optusnet.com.au (mail01.syd.optusnet.com.au [211.29.132.182]) by mx1.freebsd.org (Postfix) with ESMTP id EB84C13C478 for ; Mon, 10 Dec 2007 11:22:16 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from besplex.bde.org (c211-30-219-213.carlnfd3.nsw.optusnet.com.au [211.30.219.213]) by mail01.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id lBABMBnL006034 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Mon, 10 Dec 2007 22:22:13 +1100 Date: Mon, 10 Dec 2007 22:22:11 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl In-Reply-To: <20071209213450.GA95726@troutmask.apl.washington.edu> Message-ID: <20071210213513.Q30644@besplex.bde.org> References: <20071103001305.GA82775@troutmask.apl.washington.edu> <20071209212505.GA9698@VARK.MIT.EDU> <20071209213450.GA95726@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@FreeBSD.org Subject: Re: Implementation of expl() X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Dec 2007 11:22:17 -0000 On Sun, 9 Dec 2007, Steve Kargl wrote: > On Sun, Dec 09, 2007 at 04:25:05PM -0500, David Schultz wrote: >> On Fri, Nov 02, 2007, Steve Kargl wrote: >>> With all the success of developing several other missing >>> C99 math routines, I decided to tackle expl() (which I >>> need to tackle tanhl()). >> >> Hmm, great, but where's the patch? :) Maybe the mailing list >> software ate it. > > This is the current version. I need to revise how I computed > the ploynomial coefficient. In short, I mapped r in > [-ln(2)/2:ln(2)/2] into the range x in [-1,1] for the Chebyshev > interpolation, but I never scaled x back into r. This is the > reason why the lines "r = r * TOLN2;" exists. > > I don't remember if bde sent me comments on this code. I sure > he has plenty. :) I don't think I commented about expl() before. Sorry, it's further from what I would like than the trig functions. Here are some better polynomial approximations for it. ---------------------------------------------------------------------------- (1) 64-bit precision, approximating exp(), saves a term and should be more accurate. Note the difference in the plots between the Chebyshev approximation and the final approximation. The uglyness and extra error in the latter are due to rounding the very-high-precision coeffs used in the former to 64 bits. The final approximation uses a Remez algorithm on rounded coeffs to try to choose the best set of rounded coeffs (!= the best set of coeffs, rounded). Its plot is ugly because my algorithm is unstable for functions that aren't even, so I had to stop it after 1 significant iteration (adjusting the P3 term from 0xaaaaaaaaaaaaaaab.0p-66 to ...aa8.0p-66). If it worked properly then it would push all the extrema back to the top and bottom of the plot, and be slightly more accurate than the Chebyshev approximation. Remez gains very little over Chebyshev in infinite precision (about a factor of 2), but can compensate for huge errors (a factor of 100[0]) due to unfortunate rounding of coeffs. ---------------------------------------------------------------------------- --- Simple Lagrange interpolation at Chebyshev points --- 5.975e-24 |'''''x'''''''''x'''''''''''x''''''''''''x''''''''''_x''''''_''" | : x " " :: : : |_ : : : _ " : " : : : :: : |: : x : : : : : : : : :: : |: " : : : : " : " : x : x : |: : : " : _ : : : " : : : : |: : : : : : : " : : : : : : |: : : : " : : : : : : : : : |: : : : : : x : : : : _ : : |: : : : : : : : " : : : : : |:: : x _ : _ : : : : : : ::| `::`:```:````:`````:`````:``````:`````"`````:`````"````_```:`::` : : : : : _ : : : : : : : :| : : " : : : : x : : : : : :| : : : : : : : : : _ : : : :| : : : : : : _ : : : : : : :| : : : : _ : : : " : _ : : :| : x: " : x : _ : : : : : :| : :: : : : : : : : : _ " :| : : :: : x : _ x : : : _| : _ :_ _ _ _ : | -5.95e-24 x........."...........x............x...........x........._.....| -0.34658 0.34658 --- Lagrange/Chebyshev with coeff adjusted to 64-bit precision --- adjn n = 64 k = 3 j = 120 toter 0 --- rounded coeff --- P0 = 1.000000000000000000000000000000000000000, /* 0x8000000000000000.0p-63 */ P1 = 1.000000000000000000000000000000000000000, /* 0x8000000000000000.0p-63 */ P2 = 0.5000000000000000000000000000000000000000, /* 0x8000000000000000.0p-64 */ P3 = 0.1666666666666666666310000000000000000000, /* 0xaaaaaaaaaaaaaaa8.0p-66 */ P4 = 0.04166666666666666666100000000000000000000, /* 0xaaaaaaaaaaaaaaa9.0p-68 */ P5 = 0.008333333333333338066950000000000000000000, /* 0x8888888888889e5d.0p-70 */ P6 = 0.001388888888888889569550000000000000000000, /* 0xb60b60b60b60cf28.0p-73 */ P7 = 0.0001984126984124767297770000000000000000000, /* 0xd00d00d00c013acd.0p-76 */ P8 = 0.00002480158730156209466100000000000000000000, /* 0xd00d00d00c1851e1.0p-79 */ P9 = 0.000002755731927361941412700000000000000000000, /* 0xb8ef1d304cd05f90.0p-82 */ P10 = 0.0000002755731927069719307580000000000000000000, /* 0x93f27dbffa11c00d.0p-85 */ P11 = 0.00000002505205082467052704390000000000000000000, /* 0xd7320ad84885e745.0p-89 */ P12 = 0.000000002087671071533106249630000000000000000000, /* 0x8f76b2a8ea87f9b3.0p-92 */ P13 = 1.60924147227467723196999999999999999999 E-10, /* 0xb0f01edf2f5fac20.0p-96 */ P14 = 1.14941936278242808142000000000000000000 E-11, /* 0xca353f02fa6b7741.0p-100 */ 7.252e-24 |''''''''''''''''_''''''''''''''''''''_""_'''''''''''''''''''''| | " x x | | x _ x _" | | " x " | --------------_---------------_xxx"--------x------_------------- | _x : " _" _ | | : : _ x x _ : | | : " _ x : | | _ : "__x x _ : x | | : _ : x_ x : | | : " : " : | |_ : xx : : x | |: : _ : : x |: x " : : |:: : x : : |::: :: : _: :| : _ :| : :| : x| : | -2.79e-23 _..............................................................| -0.34658 0.34658 adjn n = 64 k = 3 toter 0 ---------------------------------------------------------------------------- (2) 113-bit precision, approximating exp(). ---------------------------------------------------------------------------- --- Simple Lagrange interpolation at Chebyshev points --- 3.239e-38 x''''x'''''x''''''_''''''''x''''''''"''''''''_''''''"'''''"''''" : : : :x : : x: : : | |: : :: :: : _ _ : :: :: : _| |" : _ : : : " : : " : : : _ : :| |: : : : : : : : : : : : : : : : : :| |: : : : x : : : : : : : : x : : : :| |: : : : : " : : : : : : " : : : : :| | : : _ : : : " : : : : " : : : x : : | | : _ : : : : : : : : : : : : : : _ : | | : : : : : : : _ " " _ : : : : : : : | | : : : : : : : : : : : : : : : : : : | ``x`:`:``:```:``:````:```:```:````:```:```:````:``:```:``:`:`x`` | : : : : : : : : : : : : : : : : : : | | : : : x : : : : : : : : : : x : : : | | : : : : x x : : : : : : x _ : : : : | | :: : : : : " : : : : " : : : : :: | | : : : : : : : x x : : : : : : : | | : :: : : : " : : " : : : :: : | | : _: :: : : : : : : :: :: : | | : :: :: :: :: :: :: :" : | | : : :_ _: :: :_ :: : : | -3.13e-38 |..x....x....."........x.......xx......._......."".....x...._..| -0.34658 0.34658 --- Lagrange/Chebyshev with coeff adjusted to 113-bit precision --- adjn n = 113 k = 1 j = 120 toter 0 --- rounded coeff --- P0 = 1.000000000000000000000000000000000000000, /* 0x10000000000000000000000000000.0p-112 */ P1 = 1.000000000000000000000000000000000000000, /* 0x10000000000000000000000000000.0p-112 */ P2 = 0.4999999999999999999999999999999999500000, /* 0x1ffffffffffffffffffffffffffff.0p-114 */ P3 = 0.1666666666666666666666666666666666600000, /* 0x15555555555555555555555555555.0p-115 */ P4 = 0.04166666666666666666666666666668415500000, /* 0x155555555555555555555555560af.0p-117 */ P5 = 0.008333333333333333333333333333336938300000, /* 0x11111111111111111111111111a6d.0p-119 */ P6 = 0.001388888888888888888888888886437923800000, /* 0x16c16c16c16c16c16c16c15fa9389.0p-122 */ P7 = 0.0001984126984126984126984126980627689100000, /* 0x1a01a01a01a01a01a01a0191e8217.0p-125 */ P8 = 0.00002480158730158730158730175735685579700000, /* 0x1a01a01a01a01a01a01ad93230fd1.0p-128 */ P9 = 0.000002755731922398589065255750604381149600000, /* 0x171de3a556c7338faac28603831fb.0p-131 */ P10 = 0.0000002755731922398589065187949851278033800000, /* 0x127e4fb7789f5c72e69d53c0ac1d3.0p-134 */ P11 = 0.00000002505210838544171877444493611961585000000, /* 0x1ae64567f544e38fdb412a0eb7de0.0p-138 */ P12 = 0.000000002087675698786810064988525091623603400000, /* 0x11eed8eff8d8981cadc22e13e3610.0p-141 */ P13 = 1.60590438368216158655710381729946340000 E-10, /* 0x16124613a86d09fa08f48ca8c4d5e.0p-145 */ P14 = 1.14707455977270929588227069467205860000 E-11, /* 0x193974a8c07640267fec271f370a0.0p-149 */ P15 = 7.6471637318180850420027983995938208999 E-13, /* 0x1ae7f3e733b16c573381f7ce23115.0p-153 */ P16 = 4.77947733504221136080232269675803310000 E-14, /* 0x1ae7f3e773e8b2e2b7a23c9a525ff.0p-157 */ P17 = 2.81145725589067097361019490883767970000 E-15, /* 0x1952c7706c73b79a1577d9ebcf121.0p-161 */ P18 = 1.56191903751316458500805670560905000000 E-16, /* 0x168276d284fb449ca1ff48bfbb891.0p-165 */ P19 = 8.2206265778415557385251561200634072999 E-18, /* 0x12f499f725c9f790053993c3f2b17.0p-169 */ P20 = 4.11616915419431123065482372671396870000 E-19, /* 0x1e5f394471e23048e2a82e8f22b34.0p-174 */ P21 = 1.9600698694144411728713343986 E-20, /* 0x1723f29b62196e8abacb3c5c40000.0p-178 */ 9.874e-37 "''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''| |_ | | | | _ | | : | | : | | "x | | " | | _ | | | | _ | | | | "x_ | | x | | x | | _ _x | x_ x | | "_ xx _" x | | """x " "x | | x x"x___" | | " __ __ _" | 0 ,,,,,,,,,,,,,,,,,,,,,,"xxx",,"x__x",,"x__x,,,,,,,,,,,,,,,,,,,,,, -0.34658 0.34658 adjn n = 113 k = 1 toter 0 ---------------------------------------------------------------------------- (3) 64-bit precision, approximating nqexp() := x * (exp(x) + 1) / (exp(x) - 1) exp(x) can be recovered from nqexp(), and the latter has much better numerical behaviour. E.g., it needs only 7 coeffs instead of 15. The Remez search completed and equalized all the extrema. It's interesting that The Chevyshev approximation equalized the extrema for exp() but not for nqexp(). ---------------------------------------------------------------------------- --- Simple Lagrange interpolation at Chebyshev points --- 2.327e-21 "''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''" : : : _" "_ : : :: :: : : :: :: : : : : "" "" : : : : : : x _ _ x : : : : : x : : : : x : : : : : : : : : : : : : : : _ " __ __ " _ : : : : : : : " "x_ _x" " : : : : `:_```:```:``````x``````"``````""``````"``````x``````:```:```_:` |:: : : x x : : ::| |:: _ _ x _ _ x _ _ ::| |:: : : _ _ _ _ : : ::| |:: : : " " : : ::| |: : : : : :| |: : x x : :| |: x x :| |: x x :| |: :| -2.23e-21 |_............................................................_| -0.34658 0.34658 --- Lagrange/Chebyshev with coeff adjusted to 64-bit precision --- adjn n = 64 k = 13 toter 0 C-1 = 2.000000000000000000000000000000000000000, /* 0x8000000000000000.0p-62 */ C0 = 0.1666666666666666660071100000000000000000, /* 0xaaaaaaaaaaaaaa7a.0p-66 */ C1 = -0.002777777777777671015759280000000000000000, /* -0xb60b60b60b5904a2.0p-72 */ C2 = 0.00006613756613181145420083840000000000000000, /* 0x8ab355dfd4d5cfe3.0p-77 */ C3 = -0.000001653439010798086927524850000000000000000, /* -0xddebbb58747e58c9.0p-83 */ C4 = 0.00000004175172569487852482302880000000000000000, /* 0xb35282048125b16f.0p-88 */ C5 = -0.000000001045797728772232364397340000000000000000, /* -0x8fbbbc85f0e61575.0p-93 */ 1.275e-21 "''"''''''''"''''''''''''"x''''''''''x"''''''''''''"''''''''"''" : : x _ _ x : : : :: " : : " " : : " :: : : :x : : : : : : x: : : : : : : : x x : : : : : : : : : : x " : : " x : : : : : : : _ : : : : : : _ : : : : : : : : : " " : : : : : : : " : : : : : : : : " : : : : : : x x x x : : : : : |:: : : _ : __ : _ : : ::| `::``_```:`````:``````:``````````````````:``````:`````:```_``::` |:: : : : : : : : : ::| |:: : x : : : : x : ::| |: : : : " " : : : :| |: : : : : : : : : :| |: : : " : : " : : :| |: : : : : : : : : :| |: : : : " " : : : :| |: " " : : : : " " :| |_ " : : " _| -1.28e-21 |......_.........._"........................"_.........._......| -0.34658 0.34658 adjn n = 64 ---------------------------------------------------------------------------- (4) 113-bit precision, approximating nqexp(). The Remez search was stable but made negative progress after the C3 (?) term so I stopped it early. It had only almost equalized the extrema when stopped. ---------------------------------------------------------------------------- --- Simple Lagrange interpolation at Chebyshev points --- 4.567e-37 |''''''"''''''''''''''''''''''''''''''''''''''''''''''''"''''''| | " : : " | | : : _ _ : : | | :: : : :_ _: : : :: | | :: : : :: :: : : :: | | x: : " : : x x : : " : :x | | :: _ : : : x x : : : _ :: | | :: : : " : : " " : : " : : :: | | :: : : : _ : __ __ : _ : : : :: | --:-:-:-:---:--:---x---x---_--"xx"--_---x---x---:--:---:-:-:-:-- | : : : : : : : : : : : : : : | | : : : : : : : x " " x : : : : : : : | | : :: : : : : " " : : : : :: : | |: :: : : x " " x : : :: :| |: :: " " " " :: :| |: :: : : " " : : :: :| |: _: : : :_ :| |: : _ _ : :| |: " " :| |: :| |" "| -5.88e-37 _.............................................................._ -0.34658 0.34658 --- Lagrange/Chebyshev with coeff adjusted to 113-bit precision --- adjn n = 113 k = 3 j = 120 toter 0 --- rounded coeff --- C-1 = 2.000000000000000000000000000000000000000, /* 0x10000000000000000000000000000.0p-111 */ C0 = 0.1666666666666666666666666666666662000000, /* 0x15555555555555555555555555542.0p-115 */ C1 = -0.002777777777777777777777777777553326200000, /* -0x16c16c16c16c16c16c16c16b85140.0p-121 */ C2 = 0.00006613756613756613756613752850992610100000, /* 0x11566abc011566abc0114a7e047ee.0p-126 */ C3 = -0.000001653439153439153439150292458494169500000, /* -0x1bbd779334ef0aac6588d54fb28a1.0p-132 */ C4 = 0.00000004175351397573619780544017497215914800000, /* 0x166a8f2bf70ebd9cab06da3dfcbe9.0p-137 */ C5 = -0.000000001056838027737493956831265550261534200000, /* -0x122805d6442668558fc0ca2fcc585.0p-142 */ C6 = 2.67650730612754829770358952948658760000 E-11, /* 0x1d6db2c4e01fe52c6eb7674cf86f8.0p-148 */ C7 = -6.7793605801133740547903998874808967000 E-13, /* -0x17da4e1ebc3556f34306a709a13c6.0p-153 */ C8 = 1.71721130775699965239461146049180330000 E-14, /* 0x1355864cf343fe6e71f6644193ab0.0p-158 */ C9 = -4.34912151080995129543227548203712180000 E-16, /* -0x1f56b690b4b95dd505b8a8e6666b1.0p-164 */ C10 = 1.08203893915267000322982006799182099999 E-17, /* 0x18f33b03a36ff97329d6990c65e9c.0p-169 */ 3.741e-37 |''''''"''''''''''''''''''''''''''''''''''''''''''''''''"''''''| | : _ _ : | | :: : _ " " _ : :: | | :: :_ ": :" _: : : | | _ x : : : :: " " " " :: : : : x _ | | : : : : : : : : : : : : : : : : : : | | :" : : _ : : : : : : : : : : x _ : ": | | :: : " : : : x : : : : x : : : : : :: | | :: : : : : : : : " " : : : : : : : :: | | :: : : : : " : x : : x : " : : : : :: | | :: : : : : : : : :: : : : : :: : :: | --::-:---::---x---:---:---:----""----:---:---:---x---::---:-::-- | :: : :: : : : : : : : : :: : :: | |: : : :" : : : : : : : : "_ : : :| |: :: " : : _ : : _ : : :: :| |: :: : : : : : : : : :: :| |: :: : x : " " : x : :: :| |: :_ : : : : : : : : _: :| |: :: _: :: :: :_ :: :| |: : : _: :_ : : :| |: : _ x x _ : :| -3.37e-37 _x.."......................................................"..x_ -0.34658 0.34658 adjn n = 113 k = 3 toter 0 --- Attempts to use 64-bit coeffs for the 113-bit approximation were unsuccessful -- the error was about 5e-24, but about 2^-(113+5) = 3e-36 is needed. The interval is just too wide for this to work. Probably only the first couple of coeffs need to have full precision. Bruce