[Bug 237800] pow(3) returns inaccurate results

[bugzilla-noreply at freebsd.org]

https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=237800

--- Comment #5 from sgk at troutmask.apl.washington.edu ---
(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

-- 
You are receiving this mail because:
You are the assignee for the bug.