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

Aleksej Lebedev root at zta.lk
Fri Sep 23 00:07:07 PDT 2016


Actually, it doesn't matter much. I just noticed that fpgetprec is 
deprecated while new API fesetround doesn't allow to set up 80-bit 
precision (63-bit mantissa).

So, I fixed the test by simply switching to powl.

--
Aleksej Lebedev

On 22/09/16 20:18, Aleksej Lebedev wrote:
> 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



More information about the Users mailing list