Floating point limitations can wreak havoc with loops if you're not careful. What would you expect to be the output of the following loop?

#include <stdio.h>
#include <sysexits.h>

int     main(int argc,char *argv[])

{
    double  x;
    
    x = 0.0;
    while ( x != 1.0 )
    {
        printf("%0.16f\n", x);
        x += 0.1;
    }
    return EX_OK;
}

! Main program body
program zero2one
    implicit none
    real(8) :: x
    
    x = 0.0d0
    do while ( x /= 1.0d0 )
        print *, x
        x = x + 0.1d0
    enddo
end program

It would seem reasonable to expect that it prints the numbers from 0.0 through 0.9 in increments of 0.1.

0.0000000000000000
0.1000000000000000
0.2000000000000000
0.3000000000000000
0.4000000000000000
0.5000000000000000
0.6000000000000000
0.7000000000000000
0.8000000000000000
0.9000000000000000
       

The actual results are shown below:

0.0000000000000000
0.1000000000000000
0.2000000000000000
0.3000000000000000
0.4000000000000000
0.5000000000000000
0.6000000000000000
0.7000000000000000
0.7999999999999999
0.8999999999999999
0.9999999999999999
1.0999999999999999
1.2000000000000000
1.3000000000000000
1.4000000000000001

and on and on...
        

What causes this problem? How can it be avoided? In this case, the problem is caused by the fact that 0.1₁₀ cannot be represented in binary with limited digits. It's much like trying to represent 1/3 in decimal: it requires an infinite number of digits.

The ==, !=, and /= operators should never be used with floating point. Because floating point values are prone to round-off error, floating point comparisons must always include a tolerance. We might be tempted to think the code below avoids the problem, but it does not.

    double  x;
    
    x = 0.0;
    while ( x < 1.0 )
    {
        printf("%0.16f\n", x);
        x += 0.1;
    }
        

real(8) :: x

x = 0.0d0
do while ( x < 1.0d0 )
    print *, x
    x = x + 0.1
enddo
        

The problem with this code is that different floating point hardware may produce different round off errors. For example, if round-off results in a value less than 1.0, then this loop will perform one extra iteration. The last value printed will be approximately 1.0 instead of approximately 0.9 as we intended and the final value of x will be near 1.1 instead of 1.0.

The real solution is to use a tolerance. In this case, we want to stop when x is approximately 1.0. We should always choose the largest tolerance that will work, in case round-off accumulates to a high level. Using half the increment will work for this loop. It will stop the loop when x is between 0.95 and 1.05.

    double  x, increment, tolerance;
    
    increment = 0.1;
    tolerance = increment / 2.0;

    x = 0.0;
    while ( ABS(x - 1.0) > tolerance )
    {
        printf("%0.16f\n", x);
        x += 0.1;
    }
        

real(8) :: x, increment, tolerance

increment = 0.1
tolerance = increment / 2.0d0

x = 0.0d0
do while ( abs(x - 1.0) > tolerance )
    print *, x
    x = x + increment
enddo
        

Using a tolerance this way will slow down the code somewhat, since it has to calculate the absolute value of the difference between x and the ending value each iteration. However, this is sometimes the nature of programming, especially when dealing with an imperfect system such as floating point.