Bug 237800

Comment 1 Xin LI 2019-05-08 18:02:13 UTC

Confirmed on amd64 -CURRENT.

Probably related: Android's assembler version (https://android.googlesource.com/platform/bionic/+/refs/heads/master/libm/x86_64/e_pow.S) would give the correct result.  These routines are 3-clause BSD licensed and we should probably take a look as well.

troutmask:kargl[219] ./testd -a 10 29
libm =  1.0000000000000001e+29, /* 0x45f431e0, 0xfae6d722 */
mpfr =  9.9999999999999991e+28, /* 0x45f431e0, 0xfae6d721 */
 ulp = 0.51303

A 1/2 ULP error is certainly not a bad result given the
nuances of floating point arithmetic.

Comment 3 Peter Jeremy 2019-05-08 19:49:48 UTC

We could special case "y is a small positive integer" - that might give a better result though it's unclear in this case.

Alternatively, we could look at what perl is doing - converting between ASCII strings and floating point values is full of dragons and evaluating pow(10,i) anywhere in that process suggests it's not being done optimally.

When you say you get the correct value on Linux, are you using clang 8.0.0 on Linux or some version of gcc?  gcc will constant fold pow(10,29) using MPFR.  Clang will generated a call to pow() in libm.

In any event, the library is computing pow(10,29) as exp(29*log(10))
where special polynomial approximations are being used.  With a ULP
of 0.51, I suppose reading Goldberg's paper is merited.

(In reply to Peter Jeremy from comment #3)

Peter, special casing small positive y values is still 
fraught with floating rounding errors.  Expanding on the
toy program a bit.

#include <stdio.h>
#include <math.h>

int
main()
{
   int n;
   double x, x2, x3, a, b;

   printf("%a\n", 1.e29);        /* check printf */
   printf("%a\n", pow(10, 29));  /* gcc constant-folds */

   x = 10;
   printf("%a\n", pow(x, 29));   /* need library call */

   x3 = x * x * x;   /* x**3 */
   a = x3 * x3 * x;  /* x**7 */
   a *= a;           /* x**14 */
   b = x * a * a;    /* x**29 */
   printf("%a\n", b);

   x2 = x * x;    /* x**2  */
   a = x2 * x2;   /* x**4  */
   b = a * a;     /* x**8  */
   b = x * b * b * b * a;
   printf("%a\n", b);

   x = frexp(x, &n);
   x2 = x * x;    /* x**2  */
   a = x2 * x2;   /* x**4  */
   a *= a;        /* x**8  */
   b = x * a * a * a * (x2 * x2);
   b = ldexp(b, n*29);
   printf("%a\n", b);
}

With clang,
cc -o z a.c -lm && ./z
0x1.431e0fae6d721p+96
0x1.431e0fae6d722p+96
0x1.431e0fae6d722p+96
0x1.431e0fae6d721p+96
0x1.431e0fae6d722p+96
0x1.431e0fae6d722p+96

With top-of-tree GCC,
 ~/work/x/bin/gcc -o z a.c -lm && ./z
0x1.431e0fae6d721p+96
0x1.431e0fae6d721p+96   <-- constant-folding with mpfr
0x1.431e0fae6d722p+96
0x1.431e0fae6d721p+96
0x1.431e0fae6d722p+96
0x1.431e0fae6d722p+96

 ~/work/x/bin/gcc -o z a.c -fno-builtin -lm && ./z
0x1.431e0fae6d721p+96
0x1.431e0fae6d722p+96  <-- libm pow() called
0x1.431e0fae6d722p+96
0x1.431e0fae6d721p+96
0x1.431e0fae6d722p+96
0x1.431e0fae6d722p+96

(In reply to sgk from comment #4)

I sent an email reply yesterday, but it seems that only logging in works, so I'm doing that now

$ clang -lm pow.c && ./a.out
0x1.431e0fae6d721p+96
$ clang -v
clang version 6.0.0-1ubuntu2 (tags/RELEASE_600/final)

We have advised the module maintainer that perhaps they were expecting too much precision.  But since the Free BSD implementation varied from others, I thought I'd file a ticket.

It is understandable that you would close this as not a bug.

(In reply to Peter Jeremy from comment #3)
There's nothing in the ticket that indicates perl is using pow() to convert between strings and floating point values.  To do that, it instead uses a call to a C library function: some strtod() flavor, which particular one depends on the desired precision, and what is available on the platform, determined by configuration options at the time the perl executable is compiled.  If those options indicate to use long doubles, it will use strtodl(), for example.

However, it used to be that perl used atof() in some circumstances instead of when it otherwise could use plain strtod().

Perl does use pow() when performing an exponentiation operation and regular doubles were requested at perl compilation time.  This module does exponentiation.  The maintainer of the module discovered that the results of that operation were sometimes wrong, and further somehow figured out that if the floating point value was stringified and then reparsed by perl, that the value magically was corrected.  Statements s/he wrote in the perl language expand to do a C library sprintf and then an atof on that result.  Except that now it's sprintf followed by strtod.

But it turns out that the magic was only because of the poorer precision offered by atof(), and when perl converted to not ever use atof, the workaround failed.

(In reply to Karl Williamson from comment #6)

Karl,

Thanks for the follow-up.  FreeBSD setting up bugzilla
to kill the ability to reply via email is rather dumb
IMHO.

To the matter at hand, checking for equality between
computed floating point values is not a Good Thing.
One normally checks either the ULP as I have done, or
a relative error such as "abs(x - y) / abs(x)" is
less than some tiny value.

It seems that whomever is developing the perl
module should read Goldberg's paper.

https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

https://www.bugzilla.org/features/#email-in is probably easy to set up, but hard to defend against spam. Or it's just the lack of resources to deal with it.

(In reply to Dima Pasechnik from comment #9)
> https://www.bugzilla.org/features/#email-in is probably easy to set up, but
> hard to defend against spam. Or it's just the lack of resources to deal with
> it.

Don't know how hard it is to deal with spam.  I know that I have 
had a GCC bugzilla account for years and have fixed a thousand
or so bugs in gfortran.  GCC bugzilla manages to allow for email
replies.  A simple filter to check if the email is coming from
someone with a valid FreeBSD bugzilla account would likely kill
all spam .

I've just run into this issue while working on some fp issues in V8. I'm not skilled enough to improve the algorithm that freebsd uses here, but if I can help in any other ways let me know.

Comment 12 Steve Kargl 2019-11-05 20:44:00 UTC

(In reply to me from comment #11)
>
> I've just run into this issue while working on some fp issues in V8.
> I'm not skilled enough to improve the algorithm that freebsd uses
> here, but if I can help in any other ways let me know.
>

What do you mean that you ran into this problem?  The bugzilla
thread indicates that the OP was getting 0.51 ULP error.  That
is not inaccurate.  That is not understanding how floating 
point math works.

The test case I have at the moment is 0.6840442338072671 ** 2.4. As I understand it, it should be 0.4019777798321958 but I am getting 0.40197777983219574.

Comment 14 Steve Kargl 2019-11-05 21:51:38 UTC

(In reply to me from comment #13)
>
> The test case I have at the moment is 0.6840442338072671 ** 2.4.
> As I understand it, it should be 0.4019777798321958 but I am
> getting 0.40197777983219574.

The first number contains 16 digits. The second contains
17 digits, so cannot make a judgement on rounding.  But,

% ./testd -a 0.6840442338072671 2.4
libm =  4.0197777983219574e-01, /* 0x3fd9ba01, 0x02864520 */
mpfr =  4.0197777983219579e-01, /* 0x3fd9ba01, 0x02864521 */
 ulp = 0.51122

MPFR is an arbitrary precision library where I have it set
to use 212 bits of precision to compute pow().  pow() is
FreeBSD's libm value with 53-bits of precision.  You have
a ULP of 0.51122.  That is about as accurate as one might
expect.

See comment #8 for a pointer to Goldberg's paper.

> You have a ULP of 0.51122.  That is about as accurate as one might
expect.

I get that there's some variance here, but everything else I've tested on gives the precise result. I can loosen the tests to check with epsilon, but I thought I'd reach out here first to see if it might be changed.

Comment 16 Steve Kargl 2019-11-05 22:42:47 UTC

(In reply to me from comment #15)
>> You have a ULP of 0.51122.  That is about as accurate as one might
>> expect.
>
> I get that there's some variance here, but everything else I've
> tested on gives the precise result. I can loosen the tests to
> check with epsilon, but I thought I'd reach out here first to
> see if it might be changed.

Are you familiar with Goldberg's paper?  If not, you need to
read it.  If you are, then you need to re-read it. A ULP of
0.51122 is precise!