Does Fortran have inherent limitations on numerical accuracy compared to other languages?

Question

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:

1. The Right Way

2. The Wrong Way

3. 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.

Rasbora

Shouldn't #3 be "The Max Power way"? – Sippet

Whoops. You're right. Thanks! – Rasbora

Fortran is using a double here and I suspect Perl to do the same. But checking for an equal with a floating point number is asking for problems (unless you have a "safe" number like 0.0144). So I would say that it is probably your method of testing that is wrong. – My

This is a very interesting problem. I ran your Fortran code and got the same results. My first instinct was to suspect a base-2 conversion error in FORTRAN, but if so why isn't Perl susceptible to it? FORTRAN and Perl both implement IEEE-754. Also, a double precision number should only succumb to rounding errors in the 15-16 digit, but in your Fortran printout, the rounding error occurs much earlier. Very curious! – Erechtheus

@Gilead: Yes, I thought the same thing initially, regarding the rounding error being in the 15-16 digit and so I assumed that, at worst, the repeated multiplications were compounding the error dramatically by the time I reached ~1700. It wasn't until I printed out the entire sequence of values, from the initial value to what would have been the stopping point, that I saw that the rounding error was of far greater magnitude then E-15 even after the first multiplication. In any case, thank you for the comments. – Rasbora

Answer 1

Your assignment is accidentally converting the value to single precision and then back to double.

Try making your 0.1 * be 0.1D0 * and you should see your problem fixed.

Batch

Actually, he was using a REAL(kind=8), which is DOUBLE PRECISION. – Erechtheus

@Gilead: Yes, I missed that initially; my Fortran experience isn't with GNU Fortran. Removed bad answer and put in correct answer. – Batch

Ah! That's the answer. Very subtle! +1 for this. – Erechtheus

@ysth: Thank you so very very much for this. Due to the obvious fundamental nature of this question, I would have been royally screwed in programs of any reasonable size if I could not get this matter resolved. – Rasbora

So is the difference here between Fortran and Perl that Perl gives literals double precision by default while Fortran gives literals single precision by default? – Rasbora

Or if you want, say Perl has only one type of real while Fortran has several and "interesting" rules when they are mixed. – Batch

Always explicitly typing your constants in Fortran is a really good idea. – Batch

The Fortran rules for expressions with mixed types are simple: the less precise types are promoted to the higher type. The issue here is that constants such as "5.0" are default reals, likely single precision, and get converted to the internal representation using that precision. Only when they interact with a variable of higher precision are they promoted. – Tytybald

@M. S. B.: thanks, I knew it wasn't as stupid as "promoting" to the lower type, but couldn't remember exactly what the underlying cause was. – Batch

@M.S.B.: But then why was it necessary for me to write 0.1D0 rather than 0.1? When I perform 0.1 * MACH_2_METERS_PER_SEC, where MACH_2_METERS_PER_SEC is already double precision, wouldn't then the 0.1 already be promoted to double as well? – Rasbora

@CmdrGuard: yes, but the single precision value closest to .1 is a lot more inaccurate than the double precision value closest to .1. For example, the 0.1 in the following is promoted to double, but that doesn't magic back the lost precision: WRITE (*, "(F20.18)") (0.1-0.1D0) – Batch

The order is key: the decimal constant "0.1" is converted to binary representation using single precision. The conversion to single-precision binary will involve more of the base-conversion issues than if the conversion were to double-precision binary, which would be done if the constant were written "0.1D0" or "0.1_DR_K" in my example. Then when that single-precision binary value is used in an expression with a double-precision variable, it is converted to double precision and the entire calculation is done in double precision. But the base conversion already was done in single! – Tytybald

@Batch and @M.S.B.: Ahhh, yes!! Thanks very much guys. This clears it all up. In decimal 0.1 is precise but in binary, the same constant is repeating and so the truncation of single precision cause a larger difference from its exact decimal representation than the double precision binary number. – Rasbora

Answer 2

As already answered, "plain" floating point constants in Fortran will default to the default real type, which will likely be single-precision. This is an almost classic mistake.

Also, using "kind=8" is not portable -- it will give you double precision with gfortran, but not with some other compilers. The safe, portable way to specify precisions for both variables and constants in Fortran >= 90 is to use the intrinsic functions, and request the precision that you need. Then specify "kinds" on the constants where precision is important. A convenient method is to define your own symbols. For example:

integer, parameter :: DR_K = selected_real_kind (14)

REAL(DR_K), PARAMETER :: MACH_2_METERS_PER_SEC = 340.0_DR_K

real (DR_K) :: mass, velocity, energy

energy = 0.5_DR_K * mass * velocity**2

This can also be important for integers, e.g., if large values are needed. For related questions for integers, see Fortran: integer*4 vs integer(4) vs integer(kind=4) and Long ints in Fortran

Tytybald

real(kind=8) is perfectly standard Fortran 90. All Fortran 90 compliant compilers should, and do, support it. That said, the proper way to ask for a given precision in Fortran 90 is something like the following: integer, parameter :: dp = selected_real_kind(15, 307) real(kind=dp) :: a See the Fortran Wiki – Flats

@fB: @MSB is correct to assert that kind=8 is not portable, and you are right to state that it is standard Fortran 90. The non-portability of the statement arises because while kind=8 is a valid kind selector the intepretation of the 8 is implementation dependent. Certainly most current compilers will (I believe) interpret this to be an 8-byte integer, but this is not guaranteed by the standard. It's certainly not always been true that all Fortran compilers would interpret this the same way, and one of the concerns of the average Fortran programmer is portability across time. – Differentiation

Yeah, I realize I did use a less than ideal means to declare precision. I did so more for simplicity than anything else. Perhaps I should change it to make my example code portable to readers of Stack Overflow, but since this question could affect other beginners to Fortran, maybe it's best I leave the code in the current form. – Rasbora

@High Performance Mark - correct. There was some discussion about this in the standard committee, and it is general conclusion that integer kind numbers are often misleading as shown. For (i.e. kind=8) - 8 means nothing in particular; it is simply a number describing the implemented precision. It could've easily been 32 or 112 - just a label. The part which is important is verifying whether the processor supports your desired precision (cpu+compiler) and using it with a combination of sel...real_kind and a declaration of a varible with that type. – Singh