Re: Possible Mistype in exp.c

[Damian McGuckin <damianm@esi.com.au>]

Just for reference, I have summarized the errors from a refactored version
of exp(x) from MUSL 1.1.21, the SUN fdlibm variant. Your table driven code
will run faster than my rework and have smaller errors than this rewrite of
exp.c.  But I thought I would just put this error summary out there.

I took the original, made the polynomial calculation more super-scalar
friendly, improved the comparisons at the start to reduce their average
impact which also means I could trash using GET_HIGH_WORD, and finally
discarded the scalbn at the end of exp(x), and computed the final result
of exp(x) as

 	return h * f * (one + y)

where y is the variable in the existing MUSL exp(x). After your wise
suggestion, I have changed it to

 	return h * (f + f * y)

FMA or non-FMA seems to make little difference.

The 'f' is a constant of the form 2^j, either 2^(+(p+2)) or 2^(-(p+2))
where 'p' is the precision, depending solely on the sign of 'k' from
the MUSL routine so either 0x1.0p+55 or 0x1.0p-55. No branching is
used in its computation. That 'h' is a normal number of the form 2^i
where i is equal to k-j, i.e. k as per MUSL and 'f' from above, with
the exponent j laying within the range [-1020,+969].

I compare against the results of GLIBC's expl(x). I hope that is OK.

 	exp(x) for x in [-745.0..600.0] by (1.25e-06) 1073741824 cases
 	*epsilon  frequency  %of-total   seen within
 	[0.0,0.1) 307561520  28.64390%  -745 .. +600
 	[0.1,0.2) 304143157  28.32554%  -744 .. +600
 	[0.2,0.3) 231766165  21.58491%  -744 .. +600
 	[0.3,0.4) 146794197  13.67127%  -744 .. +600
 	[0.4,0.5)  80455276   7.49298%  -744 .. +600
 	[0.5,0.6)   2932759   0.27313%  -723 .. +600
 	[0.6,0.7)     80629   0.00751%  -710 .. +600
 	[0.7,0.8)      8121   0.00076%  -709 .. -708
 	0.8+above         0   0.00000%  +N/A .. -N/A
 	WORST ERROR < 0.8*ULP

There is no problem in the error with exp(x) when x > +600 so I
stopped including it in my tests.

I never compare exp(x) for x < -745. In fact, exp(-745.00) returns
0x1.0p-1024 whereas expl(x) returns 0.57125014747105418 times that.
As 0.57123 will round up to 1.0, I treat them as matching anyway.
Note that my revised exp(x) underflows destructively to zero at a
value of x = -745.133219101941165264, i.e. calling the refactored
version exp', I see

 	exp'(-745.133219101941165263) = 4.94065646e-324 = 0x1.0p-1074;
 	exp'(-745.133219101941165264) = 0.0 (oops, too tiny)
while
 	expl(-745.133219101941165263) = 0.50000000000002121 * 0x1.0p-1074
 	expl(-745.133219101941165264) = 0.50000000000002121 * 0x1.0p-1074

This is better than in the original version, but still not quite at the
theoretical limit which is -745.133219101941207623 = ln(2.0) * 1075.

Even though I see correct rounding at the lower limit, I do not explicitly 
try and perform correct rounding as the old one did not either.

Also, I found that the computation times varies little between the early
model Xeon E5-1660 @ 3.3Ghz and a Xeon E5-2650-v4 @ 2.1Ghz. Interesting.

Regards - Damian

Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer