A Primer on Floating-Point Arithmetic and IEEE 754

This repository addresses usnistgov/discuss#8 by exploring
floating-point rounding errors through the lens of the associative
property of addition. The programs herein demonstrate that deviations
from mathematical precision are due to the standard representation of
floating-point numbers in a binary data type, in compliance with
IEEE 754. The evidence presented here does not support any
role of the CPU, instruction sets, or dynamic execution in the
deviations.

tl;dr: If strict adherence to mathematical law is required, use a
high-precision math library. For example, instead of the
built-in data types, use

• MPFR for C/C++.
• MPMath for Python

If performance is needed, and some precision can be sacrificed in the
interests of accuracy, build your C programs with "unsafe" math
optimizations: gcc -funsafe-math-optimizations.

This primer was written by Trevor Keller trevor.keller@nist.gov.
It is intended for educational purposes only.
Please cite this work:

Trevor Keller, "A primer on floating-point arithmetic and IEEE 754."
National Institute of Standards and Technology (2018):
DOI 10.5281/zenodo.6524704.

Table of Contents

1. Three Term Addition
2. Shuffled Summation
3. Out-of-Order Execution
4. Conclusions
5. Further Reading

Dependencies

To build the C and C++ programs in this repository, you will need:

• C and C++ compilers with standard libraries, e.g. gcc
• Make build tool
• MPFR library and headers, version 4 or greater

On Debian derivatives (Ubuntu, Mint), these are all packaged:

apt install gcc libmpfr6 libmpfr-dev make

Basic Premise

Addition of real numbers is [associative][add]: inserting parentheses
must not change the computed sum. For example,

$$ (1+2)+(3+4)+5 = (1+2)+3+(4+5) = 1+(2+3)+(4+5) = 15 $$

The two's complement representation of integers naturally
supports associativity, within the representable range of $2^{n-1}-1$
for an n-bit representation. However, fractional real numbers incur
a round-off error since infinite digits cannot practically be stored.
The chief exception to this general rule is that real fractions
formed by the sum of powers of two are exactly represented, as a
natural consequence of the data format.

It should be noted that two's complement is not the internal
representation used for floating point numbers, but can be used to
visualize floating point representations.

Three Term Addition

When addition is associative, we expect the sum of three terms,
$a+b+c$, to be independent of computation as $(a+b)+c$ or $a+(b+c)$.
Due to the binary representation of floating-point numbers (specified
by IEEE 754), or due to out-of-order execution on some CPUs,
this will not be the case for numerical approximations to
floating-point summation.

Results of Three-Term Addition

The code has three variants:

• std (standard representation)
• mpf (GNU MPFR representation)
• unsafe (std with unsafe optimizations).

Use the included Makefile or the commands in the following sections
to compile the executables.

Each program will output the value of $a$ in decimal and binary form,
$a+b+c$ in binary form, and both of the associative expressions in
decimal form. When the expressions agree, a 1 is printed in the last
column; otherwise, it will be 0. The binary form of a floating-point
number is similar to an integer binary, with each bit representing a
power of two that decreases from 0 at the decimal. Therefore,

$$ 1 = 2^{ 0} = 1.00 $$ $$ \frac{1}{2} = 2^{-1} = 0.10 $$ $$ \frac{1}{4} = 2^{-2} = 0.01 $$ $$ \frac{3}{4} = \frac{1}{2} + \frac{1}{4} = 0.11 $$

Built-in floating point representation

$ make std
gcc -O3 -Wall addition.c -o std -lm && ./std

Multi-Precision Floating-point Representation

Using precision identical to the built-in float, we have the same
result:

$ make spf
gcc -O3 -Wall -D_SINGLE_ -pedantic -include "mpfr.h" addition.c -o spf -lm -lmpfr && ./spf

The MPFR library allows us to increase the precision from double
(64-bit) to quadruple (128-bit):

$ make mpf
gcc -O3 -Wall -pedantic -include "mpfr.h" addition.c -o mpf -lm -lmpfr && ./mpf

The increased precision is implemented in software, and libraries
like this allow the programmer to choose an appropriate level of
precision, the rounding scheme, and deterministic orders of
operations. However, they are necessarily slower than numerical
formats supported in hardware (like double, single, and
half-precision floats).

Unsafe floating-point representation

Per the Using the GNU Compiler Collection(GCC) §3.10,
-funsafe-math-optimizations enables optimizations that

• assume both arguments and results are valid
• violate IEEE and ANSI standards
• change the floating-point unit control word
• cause programs that rely on exact implementation of IEEE or ISO
  rules to fail
• may yield faster code for programs that do not require the
  guarantees of these specifications.

In other words, it is ill-advised for production code. However,

$ make unsafe
gcc -O3 -Wall -funsafe-math-optimizations addition.c -o unsafe -lm && ./unsafe

Discussion of Three-Term Addition

As @JRMatey-NIST noted, this is a straight-forward byproduct
of the binary representation of floating-point numbers. Exact
representation is possible for integral exponents of 2; any other
number incurs rounding error. This is good to know, especially
for codes that frequently increment large values by small amounts.

Since many scientific computing applications model diffusive
processes (heat transfer, mass diffusion, etc.), the effect is
expected to be small: the perturbations caused by round-off error
will be smoothed out by the stable numerical scheme without any
additional effort on the part of the programmer.

Shuffled Summation

Ideally, the sequence of decimals (powers of 10)

$$ \frac{1}{1000} + 9 \times \left(\frac{1}{1000} + \frac{1}{100} + \frac{1}{10} + 1 + 10 + 100 + 1000\right) = 10^4 = 10000 $$

However, due to the same floating point representation problem,
variations arise from the order of summation. As a demonstration,
this program will generate a vector of 64 numbers ($10 + 9\times 6$),
then for each of 1 million trials, the same vector gets shuffled
before summing. The histogram of values is then reported. For
additional details, see the original thread.

Similarly, the sequence of binaries (powers of 2)

$$ \frac{1}{2} + 8 \times \left(\frac{1}{16} + \frac{1}{8} + \frac{1}{4} + \frac{1}{2} + 1 + 2 + 4 + 8\right) = 2^7 = 128 $$

The program will generate a vector of 65 numbers ($9 + 8\times 7$)
and, for each of 1 million trials, shuffle the vector before
summing. Due to the exact representation of powers-of-two, only one
result (128) is expected for all million shuffles.

Results of Shuffled Summation

There are two variants, shuffle and shuffle10, which can be built
using the Makefile or the command listed below.

Floating point representation of decimal sequence $\sum 10^n$

$ make shuffle10
g++ -O3 -Wall -pedantic -std=c++11 -DDECIMAL summation.cpp -o shuffle && ./shuffle

Floating point representation of binary sequence $\sum 2^n$

Discussion of Shuffled Summation

The sequence comprised exclusively of powers of 2 is represented
exactly, i.e. 1 million repetitions produce the same result,
exactly 128, every time. The sequence of powers of 10 is
approximate, with the exact result, 10,000, computed in only 11 %
of the million repetitions.

Out-of-Order Execution

Out-of-order execution is a common method to avoid processor
pipeline stalls: instructions are fetched in the order provided by
the compiler, but executed in a way that minimizes wasted processor
cycles. While modern CPUs implement out-of-order execution with
in-order completion, a few architectures feature out-of-order
completion of the compiler-supplied instructions, e.g. Intel's
P6 (Pentium Pro/II/III), Core, and
Silvermont (Atom, Knights Landing). Since the sequence of
mathematical operations would be non-deterministic, these processors
could exhibit poor reproducibility of floating-point arithmetic.

Results of Out-of-Order Execution

Variants of the shuffled summation and
three-term addition codes can be used to
compare results between Intel Xeon (Sandy Bridge) and
Knights Landing (Silvermont) processors. Since the KNL
hardware is novel, the Intel C++ Compiler is preferred over
GCC; it sets -fp-model fast=1 by default, which is
equivalent to -Ofast in GCC and must be disabled to provide a valid
comparison using standards-compliant arithmetic. The variants, phi
and shufflePhi, can be built using the Makefile or the commands
below.

$ make phi
icc -O3 -Wall -xmic-avx512 -fp-model precise addition.c -o phi && ./phi

$ make shufflePhi
icc -O3 -Wall -pedantic -std=c++11 -DDECIMAL -xmic-avx512 -fp-model strict summation.cpp -o shufflePhi && ./shufflePhi

Discussion of Out-of-Order Execution

The results of shuffled summation, with zeros truncated, shows no difference
between the microarchitectures:

--------------------------------  --------------------------------
         Xeon E5-1650                      Xeon Phi 7210          
--------------------------------  --------------------------------
 9999.9941406250:  0.033300000 %   9999.9941406250:  0.033300001 %
 9999.9951171875:  0.624600000 %   9999.9951171875:  0.624599993 %
 9999.9960937500:  4.240700000 %   9999.9960937500:  4.240699768 %
 9999.9970703125: 15.360700000 %   9999.9970703125: 15.360699654 %
 9999.9980468750: 34.903400000 %   9999.9980468750: 34.903400421 %
 9999.9990234375: 33.413900000 %   9999.9990234375: 33.413898468 %
10000.0000000000: 11.187700000 %  10000.0000000000: 11.187700272 %
10000.0009765625:  0.235700000 %  10000.0009765625:  0.235699996 %
--------------------------------  --------------------------------

Conclusions

For both test cases — associativity of $a+b+c$ and shuffled
summations equal to 10,000 (decimal) and 128 (binary) — the programs
in this repository demonstrate that deviations from mathematically
expected results arise strictly due to the binary representation of
floating-point numbers, not computer hardware or CPU instruction
sets. While this inadequacy of IEEE 754 is nothing new, the
source code in this repository does provide simple, repeatable
examples of the phenomenon, and may be of use as a teaching aid.

Further Reading

The following resources helped me to gain understanding.

1. David Goldberg (1991):
   What Every Computer Scientist Should Know About Floating-Point Arithmetic.

   Concise summary of how the standard floating-point representation
   works, and pitfalls to avoid.

2. Donald Knuth, The Art of Computer Programming, Vol. 2, Ed. 3,
   (1997):
   §4.2.2: Accuracy of Floating Point Arithmetic

   Philosophical review of floating-point representations with
   discussion of the compromises made to arrive at a working
   standard.

3. Erik Cheever (2001):
   Representation of Numbers

   Summary, with exercises, of how the standard represents real
   numbers.

4. Rick Regan (2012):
   Why 0.1 Does Not Exist In Floating-Point

   Excellent illustration of why round-off error is a necessary evil
   in the standard representation.

5. James Demmel and Hong Diep Nguyen (2013):
   Numerical Reproducibility and Accuracy at Exascale

   Short paper summarizing reproducibility of summation algorithms.

6. James Demmel, Hong Diep Nguyen, & Peter Ahrens (2015):
   Cost of Floating-Point Reproducibility

   Detailed slides discussing sources of non-reproducibility and the
   computational cost of correcting the deficiencies.

7. Jeff Arnold (2017):
   An Introduction to Floating-Point Arithmetic and Computation

   Thorough slide deck highlighting problem areas for scientific
   computing.

8. Thomas Risse (2017):
   Better is the Enemy of Good: Unums — An Alternative to IEEE 754 Floats and Doubles

   Concise discussion of an alternative floating-point
   representation.

9. James Matey (2018):
   Re: Floating-Point Repeatability Issues on Modern Processors

   Insightful comment on the GitHub Issue with helpful
   references.