In the previous sections we discussed the three most fundamental
causes of inaccuracy in floating-point computations:
Rounding in basic operations
Exceptions and library errors
This section lists the factors that can contribute to inaccuracy
or to changes in the results of an application. We also show how
each factor is derived from the previously described three fundamental
factors.
These factors are
Conversions between binary and decimal
Compiler behavior and compiler version
Hardware version of build-time system
Operating system release of build-time system
Operating system release of run-time system
Values of certain modifiable hardware status register
fields
Conversions Between Binary and
Decimal |
 |
A conversion from decimal to binary or from binary to decimal
may take place either at compile time or at run time. A compile-time
decimal-to-binary conversion takes place when a statement like the
following is compiled:
A run-time decimal-to-binary conversion takes place when the
C library routine atof
is called, or when a statement like the following is executed:
The conversion between a decimal value and a binary floating-point
value may cause a loss of accuracy for any of three reasons:
The algorithm may not be accurate,
probably because speed/accuracy tradeoffs have been made in favor
of speed.
Not all decimal values are exactly representable
in binary—for example, 0.1.
The conversion may be specified so as to deliberately
limit precision. For example, the statements
REAL X
READ (*, '(G5.2)') X |
deliver only three decimal digits of precision (5 minus 1
for the sign and 1 for the decimal point). If the user enters -1.23456,
the last three digits are lost. The same loss of precision on output
can be particularly confusing if it causes very small values to
be printed as zero, as in the following example:
REAL X
X = 1.0E-10
WRITE(*, '(F5.2)') X |
Displaying Floating-Point Values
in Binary
It is possible in each language to read or print floating-point
variables directly in binary (actually hexadecimal), where no conversion
inaccuracies occur. Floating-point programmers should familiarize
themselves with these techniques and should also examine floating-point
variables in hexadecimal in the symbolic debugger.
The following Fortran program shows how you can display floating-point
values in hexadecimal:
Example 3-2 Sample Program: flophex.f
PROGRAM FLOPHEX
DOUBLE PRECISION X, Y
X = 1.234D0
Y = DCOS(X)
WRITE(*, *) 'Y = ', Y
WRITE(*, 10) Y
10 FORMAT(' Y = ', Z16.16)
END |
The Fortran program displays results similar to the following:
Y = .3304651080717298
Y = 3FD526571FE8C7A5 |
The following C program shows how you can use a union to display
floating-point values in hexadecimal:
Example 3-3 Sample Program: flophex.c
#include <stdio.h>
#include <math.h>
union {
double y;
struct {
unsigned int ym, yl;
} i;
} di;
int main(void)
{
double x;
x = 1.234e0;
di.y = cos(x);
printf("di.y = %18.16g\n", di.y);
printf("di.y = %08x%08x\n", di.i.ym, di.i.yl);
} |
The C program displays results similar to the following:
di.y = 0.3304651080717299
di.y = 3fd526571fe8c7a5 |
You can see that, although the last digits of the floating-point
values displayed by the Fortran and C programs are not identical,
the hexadecimal values are exactly the same.
Compiler Behavior and Compiler
Version |
|
The compiler can contribute to variations in floating-point
results in several ways, including
Decisions about constant folding and algorithms
used in constant folding
Rearrangement of operations
Like the math library functions, the parsing routines (both
the run-time routines such as READ
and scanf and
the internal routines the compiler uses) attempt to come as close
as possible to the ideal accuracy of less than 1/2 ULP and to be
bit-for-bit compatible from one release to another. However, these
routines may change slightly from time to time. Because of this,
the parsing of a floating-point constant might change if the exact
value of the constant lies extremely close to the halfway point
between two representable values.
The compiler usually
performs compile-time expression evaluation, which is commonly referred
to as constant folding. A statement like
will probably be compiled as
but, because an extra division and addition took place, each
contributing some rounding error, the result may not be bit-for-bit
identical in all cases.
|
| |
|
 | NOTE: Floating-point constants without suffixes (1.1, for
example) are evaluated differently in Fortran 90 and FORTRAN 77.
FORTRAN 77 evaluates them in double precision, while Fortran 90
evaluates them in single precision. This can lead to slight differences
in the resulting constant value. |
|
| |
|
The zeal with which the compiler pursues opportunities for
constant folding may vary from one compiler release to another and
may depend on the optimization level and other compiler options.
In general, lower optimization levels produce more repeatable results.
New versions of compilers often incorporate improved constant
folders and optimizers that compile applications into more efficient
sequences of operations. However, as we have shown, different sequences
of operations produce different results, even though the new sequence
is mathematically equivalent to the old.
Compiler Options |
|
A compiler typically has several options that affect the sophistication
of the optimization algorithms used. As these compiler options change,
the final sequence of operations produced by the compiler changes.
HP-UX compilers by default
do not reorder floating-point operations, even at high optimization
levels. The order of operations is the same as the order your program
specifies. If, however, you want to allow the compiler to reorder
floating-point operations so as to improve performance, even at
the expense of possible small differences in results, specify the
+Onofltacc option.
For example, by default the compiler orders the expression
as
But if you use +Onofltacc,
it may change the ordering to
As we showed in “How Basic Operations Affect Application
Results”, this kind of reordering has an effect
on rounding errors and consequently on the final result.
The +Ovectorize
option may also reorder operations on arrays and thus cause small
differences in results. See “Optimizing Your Program” for details.
Architecture Type of
Run-Time System |
|
PA2.0 systems use two new instructions, known collectively
as FMA (fused multiply-add) instructions.
These instructions, FMPYFADD
and FMPYNFADD,
combine a multiplication and an addition (or subtraction) into a
single operation. Use the +DA2.0
compiler option (the default on PA2.0 systems) to generate code
that uses these instructions.
For example, in the statement
the multiplication of a
and c, and the
addition of the product to b,
may all be accomplished by a single FMPYFADD
instruction. When the instruction is executed, the product of a
and c is computed
to infinite precision and added to b.
The result is then rounded according to the current rounding mode.
This means, for example, that a code sequence that might be
expected to incur two rounding errors may instead incur only one.
Therefore, the use of FMA instructions may change the results of
your application slightly.
FMA
instructions are generated by default at optimization levels of
2 and higher on PA2.0 systems. If you want your optimized code to
preserve exactly the expression semantics of your source code, specify
the +Ofltacc
option to suppress the generation of FMA instructions.
Operating System Release of Build-Time
System |
|
When the compiler performs binary-decimal conversions or constant
folding, it calls system math libraries to perform certain operations.
For example, the statement
causes the constant-folding phase of the compiler to invoke
a system logarithmic function.
New operating system releases from time to time include improved
math function libraries. The changes may be for improved performance,
accuracy, or both. In any case, if the libraries on the build-time
system change, the compiler's constant folding may yield
a different result, and the application, regardless of which system
it is run on, may also yield a different result.
Operating System Release of Run-Time
System |
|
As stated in “Operating System Release of Build-Time
System”,
different operating system releases may have different libraries.
If the libraries on the run-time system change, an application run
on that system may yield different results from those yielded on
previous releases.
If your application must produce the same results on every
run-time system, you should build your application using archive
libraries instead of shared libraries. Archive libraries cause the
math routines of the build-time system to be permanently bound to
the application, making it immune to future library changes on the
run-time system. If you use archive libraries, however, you cannot
take advantage of improvements in accuracy or performance in future
library releases unless you rebuild your application. For more information
about archive libraries, see “HP-UX Library Basics” and “Shared Libraries versus Archive
Libraries”.
Values of Certain Modifiable Hardware
Status Register Fields |
|
All HP 9000 systems have a modifiable floating-point
status register. Figure 5-1 “PA-RISC
Floating-Point Status Register (fr0L)” illustrates this register. Two of
the fields in the status register can be modified in ways that may
change the results of an application:
The D bit, which affects underflow mode
The rounding-mode field is a two-bit field that specifies
which of the four rounding modes defined by the IEEE-754 standard
is in force. The default setting for this field is round to nearest.
Changing this field changes the rounding errors developed during
computations and thus changes the results yielded by an application.
Many HP 9000 systems have
a D bit in the floating-point status register. When this bit is
set to 0 (the default), the system is fully IEEE-754 conforming.
When this bit is set to 1, the system operates on denormalized operands
as if they were zeros, and results that underflow are flushed to
0 instead of denormalized. The purpose of this bit is to improve
the performance of applications that encounter denormalized values
often, since operations on denormalized values degrade system performance.
Changes to the D bit, by virtue of its effect, can obviously have
the effect of altering application results.
You may modify the floating-point status register by using
either a compiler option or any of several library functions. For
more information, see Chapter 5 “Manipulating the Floating-Point Status
Register”.
Printable version
