pow(2) returns slightly different result compared to other platforms' implementation

[Aleksej Lebedev]

I did some investigation and found that FreeBSD and therefore DragonFly 
allows to set up rounding precision modes via fpgetprec(3), but
the rounding precision that is supposed to be enough (FP_PE) is already 
set by default.

Also, quick search showed that IEEE only requires so called exact 
rounding only for a few basic operations: +,-,*,/ and sqrt, and doesn't 
seem to say anything about pow. But I would expect to get correct (as in 
exact rounding) result when rounding precision is set to FP_PE via 
fpgetprec, which means extended rounding precision, i.e. with 64-bit 
significand.

--
Aleksej Lebedev

On 22/09/16 19:33, Aleksej Lebedev wrote:
> Hi, everyone!
>
> I'm porting ChezScheme that was recently opensourced to DragonFly.
>
> I was able to compile it, but lots of tests are failing. At least some
> of them are due to (minor) bugs in DragonFly
> (http://bugs.dragonflybsd.org/issues/2951).
>
> Also I noticed that one of the tests is failing because function pow(3)
> returns a different (wrong?) value compared to pow on other platforms.
>
> On DragonFly:
> a at kl:~/tmp$ uname -a
> DragonFly kl.zta.lk 4.6-RELEASE DragonFly v4.6.0.10.g16fba-RELEASE #10:
> Wed Aug 17 14:26:31 CEST 2016
> root at kl.zta.lk:/usr/obj/usr/src/sys/X86_64_GENERIC  x86_64
> a at kl:~/tmp$ cat powtest.c
> #include <stdio.h>
> #include <math.h>
>
> void printrep(double x)
> {
>   unsigned long r = *(unsigned long *)&x;
>   for (int i = 0; i != 64; ++i)
>     putchar((r & (1lu << (63 - i))) ? '1' : '0');
>   putchar('\n');
> }
>
> int main()
> {
>   double x, y;
>   scanf("%lf%lf", &x, &y);
>   printrep(pow(x, y));
>   return 0;
> }
> a at kl:~/tmp$ cc powtest.c -lm
> a at kl:~/tmp$ echo 10.0 -20.0 | ./a.out
> 0011101111000111100111001010000100001100100100100100001000100100
>
> On Mac OSX the same program:
> a at zdev:~/tmp$ uname -a
> Darwin zdev.local 15.6.0 Darwin Kernel Version 15.6.0: Mon Aug 29
> 20:21:34 PDT 2016; root:xnu-3248.60.11~1/RELEASE_X86_64 x86_64
> a at zdev:~/tmp$ cc powtest.c
> a at zdev:~/tmp$ echo 10.0 -20.0 | ./a.out
> 0011101111000111100111001010000100001100100100100100001000100011
>
> I am not an expert in floating point arithmetic, more importantly not an
> expert in IEEE754. Maybe this is not a bug. Looks like a different
> rounding, but I'm not sure if it is allowed by IEE754. Does anyone know
> if it's OK? I can simply use powl instead of pow, this will make the
> test pass, but I would like not to create a workaround if this is really
> a bug.
>

-- 
Aleksej Lebedev