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

[Fernando Perez]

Hi Robert and George,

We found a bug in g77 v. 3.4.4 as well as in gcc, which manifests itself in 
the following little snippet:

planck[f77bug]> cat testbug.f
        program  testbug
c
        implicit real *8 (a-h,o-z)
c
        half = 0.5d0
        x    = 0.49d0
        nnx  = 100
        iax  = (x+half)*nnx

        print *, 'Should be 99:',iax

        stop
        end

c EOF


planck[f77bug]> g77 -o testbug.g77 testbug.f
planck[f77bug]> ./testbug.g77
   Should be 99: 98

This can be seen as computing (x/n+1/2)*n and comparing it to x+n/2.  Greg is 
using this in a number of places inside a library, which had never given 
trouble before when built with other compilers, like the sun, IBM, Intel and 
Lahey ones.  Now with g77 it gives the result above.

Questions:

1. Have you seen similar behavior in the past?

2. If we switch away from g77, what do you suggest moving towards?  We ran 
paranoia on ifort, lahey and g77, and lahey was the best performing of all. 
The intel one has the advantage of being free.  On the other hand, paranoia 
did complain about arithmetic issues with it (though the above code works fine 
with intel).

Any ideas you can give us would be very appreciated.

Cheers,

Fernando and Greg.

ps. Apparently g77 v 3.3.2 does NOT have this problem.