floating point rounding problem

[Thorsten Glaser <tg@mirbsd.de>]

Hi,

I’m working on porting Python 2.7 and 3.2 to Debian/m68k and
have found out that it contains dtoa code that wants “53 bit”
IEEE 754 precision, i.e. double instead of extended. Now m68k
like i386 has… issues in that regard.

I’ve found several ways now to change the fpcr, but none of
them made Python’s test programme return 0 instead of 1 (I
normally don’t do floating point, so I don’t know if their
test is right).

i386:

tg@frozenfish:~ $ rm -f a.out; gcc -O2 t.c; ./a.out ; echo $?
0000037F 0000127F 0000127F
0

m68k:

root@ara5:~ # rm -f a.out; gcc -O2 t.c; ./a.out ; echo $?
00000000 00000080 00000080
1

The code:

root@ara5:~ # cat t.c
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <fpu_control.h>
int main() {
    volatile double x, y, z;
    fpu_control_t cwold, cwnew, cwgot;

    _FPU_GETCW(cwold);
#ifdef __i386__
    cwnew = 0x127f;
#else
    cwnew = _FPU_RC_NEAREST | _FPU_DOUBLE;
#endif
    _FPU_SETCW(cwnew);
    _FPU_GETCW(cwgot);
    printf("%08X %08X %08X\n",cwold,cwnew,cwgot);

    /* 1./(1-2**-53) -> 1+2**-52 (correct), 1.0 (double rounding) */
    x = 0.99999999999999989; /* 1-2**-53 */
    y = 1./x;
    if (y != 1.)
        exit(0);
    /* 1e16+2.99999 -> 1e16+2. (correct), 1e16+4. (double rounding) */
    x = 1e16;
    y = 2.99999;
    z = x + y;
    if (z != 1e16+4.)
        exit(0);
    /* both tests show evidence of double rounding */
    exit(1);
}


Any help appreciated. For now, I can just let Python use their
fallback (2.6) code, but this is of course not optimal. (The
test programme doesn’t even work with -ffloat-store -mieee-fp…)

bye,
//mirabilos
-- 
I believe no one can invent an algorithm. One just happens to hit upon it
when God enlightens him. Or only God invents algorithms, we merely copy them.
If you don't believe in God, just consider God as Nature if you won't deny
existence.		-- Coywolf Qi Hunt