Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!usc!jarthur!uunet!van-bc!ubc-cs!alberta!mts.ucs.UAlberta.CA!userDHAL
From: userDHAL@mts.ucs.UAlberta.CA (unknown)
Newsgroups: comp.lang.fortran
Subject: Re: Explain this output to me
Message-ID: <1074@mts.ucs.UAlberta.CA>
Date: 30 May 90 20:59:16 GMT
References: <1554@ns-mx.uiowa.edu>
Organization: University of Alberta (MTS)
Lines: 69

In article <1554@ns-mx.uiowa.edu>, jlhaferman@l_eld09.icaen.uiowa.edu (Jeffrey Lawrence Haferman) writes:
>Why does this program give me -0.2775558E-16 when I expect 0.0?
>Two different versions of Fortran have produced the same result.
>I'm not a Fortran programmer, so be easy on me.  Please explain
>what I need to do to get 0.0 where I expect it.
>
>==================================================================
>C
>      IMPLICIT REAL*8 (A-H,O-Z)
>C
>      XOUT = -5.D-1
>      DX = 1.D-1
>C
>      DO 40 IOUT = 1,50
>        IF (XOUT .GT. 5.D-1) GOTO 90
>        WRITE(6,20)XOUT
>40      XOUT = XOUT + DX
>90    STOP
>20    FORMAT(D15.7)
>      END
>C
>==================================================================
>Output:
>
> -0.5000000E+00
> -0.4000000E+00
> -0.3000000E+00
> -0.2000000E+00
> -0.1000000E+00
> -0.2775558E-16   <----- ????
 
Welcome to the fun of round-off error, Jeff. The catch here is to
understand that the computer is working in base 2, and your value of
DX is a nice round number in base 10. What you see in your output is that
you have a value of zero, IF you round it to only 15 decimal places
(which is appropriate for double precision). If you printed out all
values to 16 or more decimals, you would see errors in them as well.
Making this up as we go, you might expect to see:
-0.50000000000000000
-0.40000000000000005
-0.30000000000000010
-0.20000000000000015
-0.10000000000000020
-0.00000000000000025        (which is -0.25E-16).
 
To further clarify this (or perhaps muddy the waters), think about
how these numbers are stored in base 2: 0.5 is easy because it is 2**-1
(sticking with FORTRAN notation). 0.25 is 2**-2, 0.125 is 2**-3, etc.
What does this make 0.1? It CAN'T be represented exactly in base 2!
Therefore there HAS to be rounding error in your calculations, and you
won't get zero! Ever! Even if you do it in quintiple precision. Even on
a CYBER! This is why accounting programs should not use normal floating
calculations, with cents as decimal dollars. People lose money! Use
binary coded decimal or integer cents instead.
If you do a lot of serious numerical work, rounding error is CRITICAL.
There are ways to program to avoid it, and ways to program to guarantee
it will cause problems.
 
>  0.1000000E+00
>  0.2000000E+00
>  0.3000000E+00
>  0.4000000E+00
>  0.5000000E+00
>
>
>Jeff Haferman                            internet: jlhaferman@icaen.uiowa.edu
>Department of Mechanical Engineering
>University of Iowa
>Iowa City IA  52240