Implementation of expl()

Steve Kargl sgk at troutmask.apl.washington.edu
Fri Nov 2 17:13:53 PDT 2007


With all the success of developing several other missing
C99 math routines, I decided to tackle expl() (which I
need to tackle tanhl()).  As with the other functions, 
this implementation is known to work on amd64 and should
be valid on other 64-bit architectures.

As with the previous function, I can only do a limited
amount of statistical testing of expl().  The reference
values in testing are computing using GMP/MPFR with
265 bits of precision.  The reference value is then 
converted to a long double for comparison to my expl().

Testing shows:
Special cases:
  exp(0)   = 1e+00
  exp(inf) = inf
  exp(-inf)= 0e+00
  exp(nan) = nan

With x in the non-exceptional range of [ln(LDBL_MIN):ln(LDBL_MAX)],
I computed 22.7 million values of expl() that span the range.  The
maximum observed ULP was 1 where only 7378 values (0.0325%) had a
ULP exceeding 0.5.

In the exceptional range leading to gradual underflow, testing was
much more difficult.  Computing a ULP seemed troublesome because of
the loss of precision as the subnormal numbers approach zero.  So,
I broke the range into 4 segments and compute the relative maximum
error in expl() compared to a value computed by GMP/MPFR. 

x in [-11366.00:-11355.14]
  RME > 1e-18 --> (73524/1086300) = 6.77%   
  RME > 1e-10 --> (69314/1086300) = 6.38%
  RME > 1e-8  --> (69314/1086300) = 6.38%
  RME > 1e-4  --> (69314/1086300) = 6.38%

x in [-11377.00:-11366.00]
  RME > 1e-18 --> (2/1100000) = 0.00%

x in [-11388.00:-11377.00]
  RME > 1e-18 --> (0/1100000) = 0.00%

x in [-11399.50:-11388.00]
  RME > 1e-18 --> (156152/1149852) = 13.58%
  RME > 1e-10 --> (156152/1149852) = 13.58%
  RME > 1e-8  --> (156152/1149852) = 13.58%
  RME > 1e-4  --> (136328/1149852) = 11.86%

The above suggests that in the middle on exceptional range, my
expl() gives an accurate anwer.  In the last segment, the RME
significantly increases and is attributed to the loss of precision.

Although I'm not completely satisfied with the gradual underflow
behavior, I'm providing the code for comments and possible inclusion
into src/msun.  Note, a diff against cvs is not provided because of
the entanglement of my other patches to math.h, Symbol.map, and 
Makefile.

Enjoy.

-- 
Steve


More information about the freebsd-standards mailing list