[Numpy-discussion] A weekend floating point/compiler question
Charles R Harris
charlesr.harris at gmail.com
Sat Apr 29 10:43:08 EDT 2006
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>
More information about the NumPy-Discussion
mailing list