Re: [lmi] [lmi-commits] odd/expm1_log1p 59fff893 4/4: Test gcc builtins for certain transcendentals

[Vadim Zeitlin]

On Thu, 19 May 2022 20:11:53 +0000 Greg Chicares <gchicares@sbcglobal.net> 
wrote:

GC> On 5/19/22 16:32, Vadim Zeitlin wrote:
GC> > On Thu, 19 May 2022 06:37:38 -0400 (EDT) Greg Chicares 
<gchicares@sbcglobal.net> wrote:
GC> > 
GC> > GC> branch: odd/expm1_log1p
GC> > GC> commit 59fff893554eff66f5dd097477a204f218692d7a
GC> [...]
GC> > GC>     Test gcc builtins for certain transcendentals
GC> [...] 
GC> >  I've tried expm1(1.01) with gcc versions from 8 to 11 under Linux and 
they
GC> > all output the same value, which is the same one as you obtain
GC> > (which is 3FFBEDFB5475CB01 as raw IEEE-754 64-bit or also 
0x1.bedfb5475cb01
GC> > using %a printf() specifier), so it doesn't seem like there is any
GC> > difference between gcc versions,
GC> 
GC> That is good to know. Thanks.
GC> 
GC> > but there is indeed a difference of one
GC> > ULP between Linux and MSW. Interestingly, MSVC also returns the same value
GC> > (0x1.bedfb5475cb00) as MinGW, even though, AFAIK, MinGW does not use its
GC> > CRT for the maths functions implementation.

 I should have thought to have done it earlier but, just to close this
discussion by the always useful method of appeal to authority, according to
Wolfram Alpha[1], the value of e^1.01-1 is (after clicking "More digits")

    1.7456010150169164939897763166603876240737508195959622916673980879...

which, according to my favourite online IEEE 754 calculator[2], is
BEDFB5475CB00B6F569D6F31D238 in hexadecimal and rounds to 3FFBEDFB5475CB01
and so glibc is right here and MinGW and MSVC are both off by <1 ULP.

[1]: https://www.wolframalpha.com/input?i2d=true&i=Power%5Be%2C1.01%5D-1
[2]: https://babbage.cs.qc.cuny.edu/IEEE-754/


 FWIW I've also checked this under Mac using M1 CPU and the result is
correct there too, even though the implementation is completely different
(I don't read ARM disassembly well, but expm1() there is clearly is quite
different from x86 code).

GC> I carefully followed, and sometimes participated in, the development
GC> of MinGW's C99 additions. If they're using code copied from illustrious
GC> specialists of the 1970s, then their implementation may be highly
GC> accurate, and perhaps even more accurate than glibc's. Later development
GC> became less fastidious than I thought appropriate for standard library
GC> implementation.
GC> 
GC> Here are some people who do seem fastidious:
GC>   https://hal.inria.fr/hal-03141101/document
GC> When they say, for example, that glibc-2.34's expm1(long double) is
GC> accurate within 0.914 ulp, I would be inclined to suspect they're right,
GC> at least if their final paper includes all the code needed to reproduce
GC> their findings.

 They don't quite say this, they say that it's accurate to no more than
0.914 ULP because they didn't do an exhaustive search except in the single
precision case ("academic resources are insufficient") and so there is no
guarantee that this is really the maximal error. But at least our first
test doesn't invalidate this hypothesis.

 In any case, thanks for this link, this seems useful and I might return to
it in the future and try reproducing/comparing with their results.

 Regards,
VZ

pgpd9Jx96Hvxv.pgp
Description: PGP signature