[PATCH] implementation for expm1l()

Steve Kargl sgk at troutmask.apl.washington.edu
Wed Sep 26 03:14:09 UTC 2012


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


More information about the freebsd-numerics mailing list