[Numpy-discussion] A weekend floating point/compiler question

Fernando Perez Fernando.Perez at colorado.edu
Sat Apr 29 16:16:10 CDT 2006


Charles R Harris wrote:

>>I don't see why the answer should be 99. The number .99 can not be exactly
>>represented in IEEE floating point, in fact it is ~
>>0.9899999999999999911182. So as you can see the result is perfectly
>>correct given the standard conversion to int by truncation. IMHO, this is
>>programmer error, not a compiler problem and should be fixed in the code.
>>Now you may get slightly different results depending on roundoff error if
>>you indulge in such things as (.5 + .49)*100 vs (.33 + .17 + .49)*100, and
>>since these numbers are constants they may also be precomputed by the
>>compiler and the results will depend on the accuracy of the compiler's
>>computation. The whole construction is ambiguous.
>>
>>Chuck
>>
> 
> 
> As an example:

[...]

Thanks to yours and the other replies.  I did try resetting the FPU control 
word as suggested to only 64 bits, and in fact the 'problem' does disappear, 
and I suspect that's also why Robert sees differences in CPUs without the 
extra 16 internal FPU bits.

I do agree that I don't like code like this, but unfortunately this one is 
outside of my control.  For the sake of completeness (since this thread has 
some educational value on the vagaries of FP arithmetic), I've slightly 
extended your example to:

abdul[f77bug]> cat print99.c
#include <stdio.h>

int main(int argc, char** argv)
{
    int x = 100;

    float fy = .49;
    float fz = .50;
    float fw = (fy + fz)*x;
    int ifw = fw;


    double y = .49;
    double z = .50;
    double w = (y + z)*x;
    int iw = w;

    long double ly = .49;
    long double lz = .50;
    long double lw = (ly + lz)*x;
    int ilw = lw;

    printf("floats:\n");
    printf("w=%25.22f, iw=%d\n", fw,ifw);

    printf("doubles:\n");
    printf("w=%25.22f, iw=%d\n", w,iw);

    printf("long doubles:\n");
    printf("w=%25.22Lf, iw=%d\n", lw,ilw);

    return 0;
}
// EOF

which gives on my box (AMD chip, running 32-bit fedora3):

abdul[f77bug]> ./print99.gcc
floats:
w=99.0000000000000000000000, iw=99
doubles:
w=99.0000000000000000000000, iw=99
long doubles:
w=98.9999999999999991118216, iw=98


This is consitent with the calculations done in 80 bits giving also different 
results.

One of the nice things about this community is precisely this kind of friendly 
expertise.  Many thanks to all.

Cheers,

f




More information about the Numpy-discussion mailing list