HP-UX Floating-Point Guide

86 Chapter 3
Factors that Affect the Results of Floating-Point Computations
Floating-Point Coding Practices that Affect Application Results
Consider again the example, rounderr.f, in “How Basic Operations
Affect Application Results” on page 71. A better way to code that
example is to test whether the two values are sufficiently close to each
other, rather than exactly the same:
Sample Program: roundeps.f
PROGRAM ROUNDEPS
REAL A, B, C, D, E, F, EPSILON
PRINT *, 'Enter epsilon:'
READ *, EPSILON
PRINT *, 'Enter 4 reals:'
READ *, A, B, C, D
E = (A + B) * (C + D)
F = (A * C) + (A * D) + (B * C) + (B * D)
IF (ABS(E - F) .LE. EPSILON) THEN
WRITE (*, 20) E, ' equals ', F
ELSE
WRITE (*, 20) E, ' not equal to ', F
WRITE (*, *) 'Math error!'
ENDIF
20 FORMAT(F, A, F)
END
You must choose a value for the error bound EPSILON that is appropriate
to the magnitudes of the values being computed. If computations are
yielding results with magnitudes around m, then for single-precision
computations, a reasonable value for EPSILON might be m/1.0e5; for
double-precision computations, a reasonable value might be m/1.0e14;
for quad-precision computations, a reasonable value would be m/1.0e26.
Choosing an appropriate value for the error bound, however, requires a
thorough knowledge of the mathematical nature of your application. For
example, a value of 1.0e-5 for EPSILON yields the following result,
which may or may not be acceptable in your application:
$ f90 roundeps.f
roundeps.f
program ROUNDEPS
23 Lines Compiled
$ ./a.out
Enter epsilon:
1.0e-5
Enter 4 reals:
1.1 2.2 3.3 4.4
25.4100018 equals 25.4099998
The actual difference between the values is one bit (in hexadecimal,
41CB47AF versus 41CB47AE).