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

[Charles R Harris]

On 4/28/06, Fernando Perez <Fernando.Perez at colorado.edu> wrote:
>
> 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


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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20060429/c978dc43/attachment.html>