[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Date: Sun, 3 Feb 2019 14:07:58 +1100 (AEDT)
From: Damian McGuckin <damianm@....com.au>
To: musl@...ts.openwall.com
Subject: Re: Possible Mistype in exp.c
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
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
