[comp.lang.fortran] precision of expression evaluation

tff@na.toronto.edu (Tom Fairgrieve) (11/29/88)

I have a question for a FORTRAN standard(s) expert.  Consider the following 
"machine epsilon" program segment:

	real    a
	a = 2.0e0
   10   a = a / 2.0e0
	if ((1.0e0 + a) .gt. 1.0e0) goto 10
	write (*,*) ' a = ', a
	stop
	end

When compiled using "f77 -fsoft" and run on a SUN 3/60 under SunOS 4.0_Export, 
the output is
 a =     5.96046E-08

When compiled using "f77 -f68881" and run on the same machine, the output is
 a =     5.42101E-20

The comparison in the "if"-statement is being executed in the precision of the 
mc68881 floating point coprocessor, not the precision which was intended.

A "fix" to the above problem is to force the sum "1.0e0 + a" to be rounded
through a store operation.  i.e., 
	real    a, b
	a = 2.0e0
   10   a = a / 2.0e0
	b = 1.0e0 + a
	if (b .gt. 1.0e0) goto 10
	write (*,*) ' a = ', a
	stop
	end

This program produces the desired results, regardless of the chosen float
option.  However the forced store operation will be undone by the optimizer 
(i.e. f77 -f68881 -O) and the "too small" results reappears. The "-fstore"
option must be used to prevent this.

Do the current or proposed FORTRAN standards say anything about what precision
expressions may be executed in?  Should they?  I know why the highest precision 
is used, but is this a safe default for the non-expert programmer?  

Tom Fairgrieve               tff@na.utoronto.ca, uunet!na.utoronto.ca!tff
Department of Computer Science, University of Toronto, Toronto, Canada M5S 1A4