From mboxrd@z Thu Jan 1 00:00:00 1970 X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/13709 Path: news.gmane.org!.POSTED.blaine.gmane.org!not-for-mail From: Damian McGuckin Newsgroups: gmane.linux.lib.musl.general Subject: Re: Possible Mistype in exp.c Date: Sun, 3 Feb 2019 14:07:58 +1100 (AEDT) Message-ID: References: <20190129110135.GC21289@port70.net> <20190129114308.GD21289@port70.net> <20190130093738.GE21289@port70.net> <20190130121947.GF21289@port70.net> <20190130141040.GG21289@port70.net> <20190202033929.GK21289@port70.net> Reply-To: musl@lists.openwall.com Mime-Version: 1.0 Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII Injection-Info: blaine.gmane.org; posting-host="blaine.gmane.org:195.159.176.226"; logging-data="63945"; mail-complaints-to="usenet@blaine.gmane.org" User-Agent: Alpine 2.02 (LRH 1266 2009-07-14) To: musl@lists.openwall.com Original-X-From: musl-return-13725-gllmg-musl=m.gmane.org@lists.openwall.com Sun Feb 03 04:08:18 2019 Return-path: Envelope-to: gllmg-musl@m.gmane.org Original-Received: from mother.openwall.net ([195.42.179.200]) by blaine.gmane.org with smtp (Exim 4.89) (envelope-from ) id 1gq88n-000GWF-7n for gllmg-musl@m.gmane.org; Sun, 03 Feb 2019 04:08:17 +0100 Original-Received: (qmail 21654 invoked by uid 550); 3 Feb 2019 03:08:15 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Original-Received: (qmail 21636 invoked from network); 3 Feb 2019 03:08:14 -0000 X-Authentication-Warning: key0.esi.com.au: damianm owned process doing -bs In-Reply-To: Xref: news.gmane.org gmane.linux.lib.musl.general:13709 Archived-At: 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