From owner-freebsd-numerics@FreeBSD.ORG Wed Sep 26 03:14:09 2012 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id C3A08106564A for ; Wed, 26 Sep 2012 03:14:09 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) by mx1.freebsd.org (Postfix) with ESMTP id 5FBD08FC08 for ; Wed, 26 Sep 2012 03:14:09 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.5/8.14.5) with ESMTP id q8Q3E8eX033714; Tue, 25 Sep 2012 20:14:08 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.5/8.14.5/Submit) id q8Q3E8u9033713; Tue, 25 Sep 2012 20:14:08 -0700 (PDT) (envelope-from sgk) Date: Tue, 25 Sep 2012 20:14:07 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20120926031407.GA33693@troutmask.apl.washington.edu> References: <20120925154606.GA28919@troutmask.apl.washington.edu> <20120926021900.R2081@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20120926021900.R2081@besplex.bde.org> User-Agent: Mutt/1.5.21 (2010-09-15) Cc: numerics@freebsd.org Subject: Re: [PATCH] implementation for expm1l() X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 26 Sep 2012 03:14:09 -0000 On Wed, Sep 26, 2012 at 04:39:36AM +1000, Bruce Evans wrote: > On Tue, 25 Sep 2012, Steve Kargl wrote: > > It's expm1(x)/x that must be approximated by p(x)/x, not expm1(x) by > p(x). This combined with the nonstandard comment style confused me for > a while. (Plain expm1(x) is just as easy to approximate as exp(x) and > doesn't require any long double coefficients.) > > > + * rounded to long double precision where B5 through B14 can be used with > > + * a double type. > > + */ > > +static const long double > > +B3 = 0x1.5555555555555556p-3L, /* 1.66666666666666666671e-01 */ > > +B4 = 0x1.5555555555555556p-5L; /* 4.16666666666666666678e-02 */ > > These 2 aren't actually rounded to long double precision on i386, so the > point of having them in more precision than the others is completely lost. > gcc rounds them to nearest, so the result is the same as: > > static const non_long double > B3 = 0x1.5555555555555800p-3, /* some wrong value */ > B4 = 0x1.5555555555555800p-5; /* some wrong value */ > > This is the only obvious problem. > You nailed it on the head. I haven't tried your coefficient, yet. I just regenerated a new set and split B3 and B4 into high and low parts. With these coefficients: /* * Remes polynomial coefficients on the interval [T1:T2] with an error * estimate of log2(|f(x)-p(x)|) ~ -68 with f(x)=(expm1(x)-1-x-x**2/2)/x**3. */ static const long double B3hi = 0x1.5555555555555p-3, B3lo = 0x1.5554d1fb0e211p-57, B4hi = 0x1.5555555555555p-5, B4lo = 0x1.66bec322050bep-59, B5 = 0x1.1111111111111p-7, B6 = 0x1.6c16c16c16b69p-10, B7 = 0x1.a01a01a019ee2p-13, B8 = 0x1.a01a01a10d6f9p-16, B9 = 0x1.71de3a568cf37p-19, B10 = 0x1.27e4f2c784f73p-22, B11 = 0x1.ae644762b882bp-26, B12 = 0x1.1f3300cf9538ap-29, B13 = 0x1.6181ba21eb17dp-33; I needed to re-arrange the relevant section of code to add terms in the correct order: x2 = x * x; t1 = (long double)B3lo + x * B4lo + x2 * B5 + x * B4hi + B3hi; With these changes, I see % ./testl -m -0.25 -M 0.25 Interval tested: [-0.250000:0.250000] Max ULP: 0.54516 x Max ULP: -2.45493495493495494e-01, -0x1.f6c54b3433c96bfap-3 Count: 1000000 I test with coefficients in a bit. -- Steve