While working on a simple programming exercise, I produced a while loop (DO loop in Fortran) that was meant to exit when a real variable had reached a precise value.
I noticed that due to the precision being used, the equality was never met and the loop became infinite. This is, of course, not unheard of and one is advised that, rather than comparing two numbers for equality, it is best see if the absolute difference between two numbers is less than a set threshold.
What I found disappointing was how low I had to set this threshold, even with variables at double precision, for my loop to exit properly. Furthermore, when I rewrote a "distilled" version of this loop in Perl, I had no problems with numerical accuracy and the loop exited fine.
Since the code to produce the problem is so small, in both Perl and Fortran, I'd like to reproduce it here in case I am glossing over an important detail:
Fortran Code
PROGRAM precision_test
IMPLICIT NONE
! Data Dictionary
INTEGER :: count = 0 ! Number of times the loop has iterated
REAL(KIND=8) :: velocity
REAL(KIND=8), PARAMETER :: MACH_2_METERS_PER_SEC = 340.0
velocity = 0.5 * MACH_2_METERS_PER_SEC ! Initial Velocity
DO
WRITE (*, 300) velocity
300 FORMAT (F20.8)
IF (count == 50) EXIT
IF (velocity == 5.0 * MACH_2_METERS_PER_SEC) EXIT
! IF (abs(velocity - (5.0 * MACH_2_METERS_PER_SEC)) < 1E-4) EXIT
velocity = velocity + 0.1 * MACH_2_METERS_PER_SEC
count = count + 1
END DO
END PROGRAM precision_test
Perl Code
#! /usr/bin/perl -w
use strict;
my $mach_2_meters_per_sec = 340.0;
my $velocity = 0.5 * $mach_2_meters_per_sec;
while (1) {
printf "%20.8f\n", $velocity;
exit if ($velocity == 5.0 * $mach_2_meters_per_sec);
$velocity = $velocity + 0.1 * $mach_2_meters_per_sec;
}
The commented-out line in Fortran is what I would need to use for the loop to exit normally. Notice that the threshold is set to 1E-4, which I feel is quite pathetic.
The names of the variables come from the self-study-based programming exercise I was performing and don't have any relevance.
The intent is that the loop stops when the velocity variable reaches 1700.
Here are the truncated outputs:
Perl Output
170.00000000
204.00000000
238.00000000
272.00000000
306.00000000
340.00000000
...
1564.00000000
1598.00000000
1632.00000000
1666.00000000
1700.00000000
Fortran Output
170.00000000
204.00000051
238.00000101
272.00000152
306.00000203
340.00000253
...
1564.00002077
1598.00002128
1632.00002179
1666.00002229
1700.00002280
What good is Fortran's speed and ease of parallelization if its accuracy stinks? Reminds me of the three ways to do things:
The Right Way
The Wrong Way
The Max Power Way
"Isn't that just the wrong way?"
"Yeah! But faster!"
All kidding aside, I must be doing something wrong.
Does Fortran have inherent limitations on numerical accuracy compared to other languages, or am I (quite likely) the one at fault?
My compiler is gfortran (gcc version 4.1.2), Perl v5.12.1, on a Dual Core AMD Opteron @ 1 GHZ.