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
Depends on:
Reported: 2019-05-08 16:53 UTC by Karl Williamson
Modified: 2019-05-09 18:35 UTC (History)
4 users (show)

See Also:

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


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 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

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

 ~/work/x/bin/gcc -o z a.c -fno-builtin -lm && ./z
0x1.431e0fae6d722p+96  <-- libm pow() called
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
$ 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)


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

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.

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 .