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

[Charles R Harris]

On 4/29/06, Charles R Harris <charlesr.harris at gmail.com> wrote:
>
>
>
> 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
>

As an example:

#include <cstdio>

int main(int argc, char** argv)
{
    int x = 100;
    long double y = .49;
    long double z = .50;
    printf("%25.22Lf\n", (y + z)*x);
    return 0;
}

prints 98.9999999999999991118216 whereas the same code with doubles instead
of long doubles prints 99.0000000000000000000000.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20060429/87ffebbd/attachment.html>