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);