g77 bug

[Toon Moene]

Matthew Bettencourt wrote:

> I have a very simple code (part of linpack) that seems to show a bug in
> g77.  Now, if I compile this just like
> 
> g77 main.f -O
> I get the right results
> 
> However, if I compile as follows
> 
> g77 main.f
> 
> I get the wrong results.

>       DOUBLE PRECISION EPS,HUGE,S,TINY
> C     ..
>       EPS = 1.0D0
>    10 EPS = EPS/2.0D0
>       S = 1.0D0 + EPS
>       IF (S.GT.1.0D0) GO TO 10
>       EPS = 2.0D0*EPS
> 
>       print*,eps
> 
>       stop
>       end

> {attica} 175 % g77 main.f -O
> {attica} 176 % a.out
>   1.08420217E-19
> {attica} 177 % g77 main.f
> {attica} 178 % a.out
>   2.22044605E-16
> {attica} 179 % pgf77 main.f
> {attica} 180 % a.out
>    1.0842021724855044E-019
> FORTRAN STOP

Q: Which of the numbers is the right one [ Note: the answer might not be
   the obvious :-) ]

A: None of the above; there is no "right" answer.  The Fortran Standard
   doesn't pin down the exact precision with which this computation has
   to be performed.

   2.22044605E-16 is consistent with the association "DOUBLE means
64-bit
   IEEE floating point arithmetic".
   1.08420217E-19 is consistent with the association "DOUBLE means
80-bit
   IEEE floating point arithmetic".

Besides, it isn't a Fortran "problem".  Try the following with gcc:

#include <stdio.h>

main()
{
        double s, eps;
        eps = 1.0;
L10:    eps = eps / 2.0;
        s = 1.0 + eps;
        if (s > 1.0)
                goto L10;
        eps = 2.0 * eps;
        printf("%e\n",eps);
}

[ Real Programmers can write Fortran in any language ]

-- 
Toon Moene - mailto:toon@moene.indiv.nluug.nl - phoneto: +31 346 214290
Saturnushof 14, 3738 XG  Maartensdijk, The Netherlands
Maintainer, GNU Fortran 77: http://gcc.gnu.org/onlinedocs/g77_news.html
Join GNU Fortran 95: http://g95.sourceforge.net/ (under construction)