Bug 237800 - pow(3) returns inaccurate results
Summary: pow(3) returns inaccurate results
Status: Open
Alias: None
Product: Base System
Classification: Unclassified
Component: misc (show other bugs)
Version: CURRENT
Hardware: Any Any
: --- Affects Many People
Assignee: freebsd-numerics mailing list
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2019-05-08 16:53 UTC by Karl Williamson
Modified: 2019-11-05 22:42 UTC (History)
6 users (show)

See Also:


Attachments
Reproducer .c (86 bytes, text/plain)
2019-05-08 16:53 UTC, Karl Williamson
no flags Details

Note You need to log in before you can comment on or make changes to this bug.
Description Karl Williamson 2019-05-08 16:53:46 UTC
Created attachment 204271 [details]
Reproducer .c

The attached essentially one-liner, program prints 0x1.431e0fae6d722p+96
on FreeBSD 13 clang version 8.0.0 (tags/RELEASE_800/final 356365) (based on LLVM 8.0.0)

On Linux, it prints 0x1.431e0fae6d721p+96
a difference of 1 in the final hex digit of the mantissa.  The Linux version is the correct value.
We are getting identical failures on 11.2, 12.0, and Open BSD 6.4

http://www.cpantesters.org/distro/M/Math-Clipper.html?grade=3&perlmat=1&patches=2&oncpan=2&distmat=2&perlver=ALL&osname=ALL&version=1.27

I work on maintaining and extending the Perl 5 programming language.  We made a change earlier in our current development cycle that improved the accuracy of converting a string representing a floating point number into an equivalent double precision C value.  We now use strtod() instead of atof().  That broke the module Math::Clipper.  It turns out it was relying on the imprecision of atof() to work around this pow() bug.  That work around was added in May 2018.

The above smoke reports rely on a volunteer network of people to do the testing.  That means the coverage of a module such as Math::Clipper is spotty.  Probably it's failing Open BSD in other versions of Perl, but those just didn't get tested.
Comment 1 Xin LI freebsd_committer 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.
Comment 2 sgk 2019-05-08 19:00:50 UTC
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 freebsd_committer 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.
Comment 4 sgk 2019-05-08 20:16:50 UTC
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.
Comment 5 sgk 2019-05-08 20:24:32 UTC
(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
Comment 6 Karl Williamson 2019-05-09 17:09:00 UTC
(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.
Comment 7 Karl Williamson 2019-05-09 17:12:35 UTC
(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.
Comment 8 sgk 2019-05-09 17:30:58 UTC
(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
Comment 9 Dima Pasechnik 2019-05-09 17:42:31 UTC
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.
Comment 10 sgk 2019-05-09 18:35:43 UTC
(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 .
Comment 11 me 2019-11-05 20:26:10 UTC
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 freebsd_committer 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.
Comment 13 me 2019-11-05 20:51:18 UTC
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 freebsd_committer 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.
Comment 15 me 2019-11-05 22:21:18 UTC
> 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 freebsd_committer 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!