
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 superscalar 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 nonFMA 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.0p55. No branching is used in its computation. That 'h' is a normal number of the form 2^i where i is equal to kj, 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.25e06) 1073741824 cases *epsilon frequency %oftotal 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.0p1024 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.94065646e324 = 0x1.0p1074; exp'(745.133219101941165264) = 0.0 (oops, too tiny) while expl(745.133219101941165263) = 0.50000000000002121 * 0x1.0p1074 expl(745.133219101941165264) = 0.50000000000002121 * 0x1.0p1074 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 E51660 @ 3.3Ghz and a Xeon E52650v4 @ 2.1Ghz. Interesting. Regards  Damian Pacific Engineering Systems International, 277279 Broadway, Glebe NSW 2037 Ph:+61285710847 .. Fx:+61296929623  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.