Variable precision in modern floating-point computing
David H. Bailey
Lawrence Berkeley Natlional Laboratory (retired)
University of California, Davis, Department of Computer Science

1 / 33

Questions to examine in this talk
I How can we ensure accuracy and reproducibility in floating-point computing?
I What types of applications require more than 64-bit precision?
I What types of applications require less than 32-bit precision?
I What software support is required for variable precision?
I How can one move efficiently between precision levels at the user level?

2 / 33

Commonly used formats for floating-point computing
Formalname
IEEE 16-bit(none)(none)IEEE 32-bit
IEEE 64-bit
IEEE 80-bit
IEEE 128-bit(none)(none)(none)

Nickname“IEEE half”“ARM half”“bfloat16”“IEEE single”“IEEE double”“IEEE extended”“IEEE quad”“double double”“quad double”“multiple”

Sign
1
1
1
1
1
1
1
1
1
1

Number of bits
Exponent Mantissa
5
10
5
10
8
7
7
24
11
52
15
64
15
112
11
104
11
208
variesvaries

Hidden
1
1
1
1
1
0
1
2
4
variesDigits
3
3
2
7
15
19
34
31
62
varies3 / 33

Numerical reproducibility in scientific computing
A December 2012 workshop on reproducibility in scientific computing, held at Brown
University, USA, noted that
Science is built upon the foundations of theory and experiment validated andimproved through open, transparent communication. ...
Numerical round-off error and numerical differences are greatly magnified ascomputational simulations are scaled up to run on highly parallel systems. As
a result, it is increasingly difficult to determine whether a code has beencorrectly ported to a new system, because computational results quickly
diverge from standard benchmark cases. And it is doubly difficult for otherresearchers, using independently written codes and distinct computer
systems, to reproduce published results.
I V. Stodden, D. H. Bailey, J. M. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default toreproducible: Reproducibility in computational and experimental mathematics,”
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.
4 / 33

Reproducibility problems in a Large Hadron Collider code
I The 2012 discovery of the Higgs boson at the ATLAS experiment in the LHC

relied crucially on the ability to track charged particles with exquisite precision (10microns over a 10m length) and high reliability (over 99% of roughly 1000
charged particles per collision correctly identified).
I Software: five million line of C++ and Python code, developed by roughly 2000
physicists and engineers over 15 years.
Recently, in an attempt to speed up the calculation, researchers found that merely
changing the underlying math library resulted in some collisions being missed or
misidentified.
Questions:
I How extensive are these numerical difficulties?
I How can numerically sensitive code be tracked down?
I How can the library be maintained, producing numerically reliable andreproducible results?
5 / 33Other numerical reproducibility issues and remedies
I Berkeley Lab climate atmospheric modeling code: Researchers found that when

their code was merely ported from one system to another, or run with a different
number of processors, computed results quickly diverged from a benchmark run.
Had they ported the code correctly? Which run was “correct”?
I Solution: Researchers found that almost all of the numerical variability
disappeared by employing double-double arithmetic (roughly 31-digit arithmetic)in two key summation loops. Total run time was unchanged.
I Livermore Lab computational physics code: Researchers working with some large
computational physics applications reported similar difficulties with reproducibility.
I Solution: Researchers solved the problem by employing a “pseudo-double-double”
software scheme in several critical summation loops.
More examples of high precision applications will be presented later in the talk.
I Y. He and C. Ding, “Using accurate arithmetics to improve numerical reproducibility and stability inparallel applications,” J. Supercomputing, vol. 18 (2001), pg. 259–277.
I R. W. Robey, J. M. Robey, and R. Aulwes, “In search of numerical consistency in parallel programming,”
Parallel Computing, vol. 37 (2011), pg. 217–219.
6 / 33

How to enhance numerical reproducibility and solve accuracy problems
1. Employ an expert numerical analyst to examine every algorithm employed in thecode, to ensure that only the most stable and efficient schemes are being used.
2. Employ an expert numerical analyst to analyze every section of code for numericalsensitivity.
3. Employ interval arithmetic for large portions of the application (which will greatly
increases run time and code complexity).
4. Employ higher-precision arithmetic, assisted with some “smart” tools to help
determine where extra precision is needed and where it is not.
#4 is the only practical solution for real-world computing.

7 / 33

Numerical analysis expertise among U.C. Berkeley graduates
Of the 2010 U.C. Berkeley graduating class, 870 were indisciplines requiring technical computing (count by DHB):
I Division of Mathematical and Physical Sciences (Math,

Physics, Statistics).
I College of Chemistry.
I College of Engineering (including Computer Science).

Other fields whose graduates will likely do significant computing:
I Finance, biology, geology, medicine, economics, psychology and sociology.
The total count is likely well over 1000; maybe closer to 2000.
Enrollment in numerical analysis courses:
I Math 128A (Introductory numerical analysis required of math majors): 219.
I Math 128B (A more advanced course, required to do serious work): 24.
Conclusion: At most, only 2.4% of U.C. Berkeley graduates who will do technicalcomputing in their careers have had rigorous training in numerical analysis.

8 / 33

The U.C. Berkeley/U.C. Davis “Precimonious” toolObjective: Develop software facilities to find and ameliorate numerical anomalies inlarge-scale computations:
I Facilities to test the level of numerical accuracy required for an application.
I Facilities to delimit the portions of code that are inaccurate.
I Facilities to search the space of possible code modifications.
I Facilities to repair numerical difficulties, including usage of high-precision
arithmetic.
I Facilities to navigate through a hierarchy of precision levels (32-bit, 64-bit, 80-bit
or higher as needed).
The current version of this tool is known as “Precimonious.”
I C. Rubio-Gonzalez, C. Nguyen, H. D. Nguyen, J. Demmel, W. Kahan, K. Sen, D. H. Bailey and C. Iancu,“Precimonious: Tuning assistant for floating-point precision,” Proceedings of SC13, November 2013,https://www.davidhbailey.com/dhbpapers/precimonious.pdf.
I C. Nguyen, C. Rubio-Gonzalez, B. Mehne, K. Sen, C. Iancu, J. Demmel, W. Kahan, W. Lavrijsen,
D. H. Bailey and D. Hough, “Floating-point precision tuning using blame analysis,” Proceedings of the
38th International Conference on Software Engineering (ICSE 2016), 14–22 May 2016,
https://www.davidhbailey.com/dhbpapers/blame-analysis.pdf.

9 / 33

Can we achieve bit-for-bit floating-point replicability between systems?
I Although the IEEE floating-point standards guarantee bit-for-bit results (provided

the same rounding mode is used) for individual operations, different systems may
produce different higher-level results, because of non-associativity:(A ⊕ B) ⊕ C 6= A ⊕ (B ⊕ C ).
I Recently some researchers have proposed to establish a standard that would

guarantee bit-for-bit replicability, first on a single system and subsequently on aparallel system. For details, contact Michael Mascagni of Florida State University.
Dangers: If this standard becomes the default on computer systems:
I Will it lock into place computations that are seriously inaccurate?
I Will it make it more difficult for a user to even be aware that he/she has aproblem with numerical sensitivity and/or inaccuracy?
Be careful what you wish for. You may get it...
10 / 33

Why are code developers reluctant to employ high-precision arithmetic?
I Many are persuaded that physical reality fundamentally does not require highprecision, beyond, say, the limits of 64-bit IEEE arithmetic. Invalid.
I Many believe that numerical problems are always due to the usage of inferioralgorithms. Invalid.
I Many complain that easy-to-use high-precision arithmetic software facilities are

not widely available. Invalid.
I Many regard the usage of high-precision arithmetic as “cheating” or “sinful.”

Why? Because it is too easy?
I Many are concerned about the 10–50X increase in run time (for IEEE quad). Thisis a legitimate concern. Can this be improved?

11 / 33

Innocuous example where standard 64-bit precision is inadequate
Problem: Find a polynomial to fit the data (1, 1048579, 16777489,
84941299, 268501249, 655751251, 1360635409, 2523398179, 4311748609) forarguments 0, 1, · · · , 8. The usual approach is to solve the linear system:






Pn n

k=1 xk..
Pn .

n
k=1 xkPn
Pnk=1 xk2
k=1 xk..
P . n+1
k=1 xk

···
···
..
.
···

Pn

  Pnna0k
Pn k=1 xn+1
Pn k=1 yk  a1  
k=1 xk yk 
k=1 xk   ..  = .
..
..
 .  Pn . nPn . 2nan
k=1 xk yk
k=1 xk

64-bit arithmetic fails to find the correct polynomial using Matlab, Linpack or
LAPACK, even if one rounds results to nearest integer.
IEEE quad or double-double arithmetic quickly produces the correct polynomial:f (x ) = 1 + 1048577x 4 + x 8 = 1 + (220 + 1)x 4 + x 8 .

12 / 33

Innocuous example, cont.
The result on the previous page can be obtained with double precision using Lagrange
interpolation or the Demmel-Koev algorithm. But few outside of expert numericalanalysts are aware of these schemes.
Besides, even these schemes fail for larger problems, such as:(1, 134217731, 8589938753, 97845255883, 549772595201,
2097396156251, 6264239146561, 15804422886323, 35253091827713,
71611233653971, 135217729000001, 240913322581691, 409688091758593).
These values are generated by the simple polynomialf (x ) = 1 + 134217729x 6 + x 12 = 1 + (227 + 1)x 6 + x 12 ,which is easily found using quad or double-double arithmetic.

13 / 33

Types of applications that may require more than 64-bit precision
1. Large-scale, highly parallel simulations, running on systems with hundreds ofthousands or millions of processors — numerical sensitivities are greatly magnified.
2. Certain large linear programming calculations.
3. Certain large applications with highly ill-conditioned linear systems.
4. Large summations, especially those involving +/- terms and cancellations.
5. Long-time, iterative simulations.
6. Studies in computational physics and experimental mathematics — tens, hundreds
or thousands of digits may be required.

As applications continue to scale up in size, numerical difficulties will grow.

14 / 33

Some specific applications of more than 64-bit precision
1. Linear programming applications in biology and other fields (34 digits).
2. Planetary orbit calculations (32 digits).
3. Supernova simulations (32–64 digits).
4. Climate modeling (32 digits).
5. Coulomb n-body atomic system simulations (32–120 digits).
6. Schrodinger solutions for lithium and helium atoms (32 digits).
7. Electromagnetic scattering theory (32–100 digits).
8. Scattering amplitudes of fundamental particles (32 digits).
9. Discrete dynamical systems (32 digits).
10. The Taylor algorithm for ODEs (100–600 digits).
11. Problems in experimental mathematics (100–100,000 digits).
I D. H. Bailey, R. Barrio, and J. M. Borwein, “High precision computation: Mathematical physics anddynamics,” Applied Mathematics and Computation, vol. 218 (2012), pg. 10106–10121.
15 / 33

Linear programming applications in biology
I Stanford researchers are using the Constraint-Based Reconstruction and Analysis(COBRA) approach, a linear programming based scheme, to study biologicalprocesses such as metabolism and macromolecular synthesis.
I Researchers found that 32-bit arithmetic is “not useful,” and that 64-bit
arithmetic “may not ensure adequate accuracy,” but IEEE 128-bit (“quad”)
produces solutions of “exceptional accuracy.”
The authors conclude: We believe that quadruple-precision solutions are now practicalfor multiscale [linear programming] applications ..., and that they justify increased
confidence as systems biologists build ever-larger models.
I Ding Ma and Michael A. Saunders, “Solving multiscale linear programs using the simplex method in
quadruple precision,” in M. Al-Baali et al. (eds.), Numerical Analysis and Optimization, Springer
Proceedings in Mathematics and Statistics, vol. 134, DOI 10.1007/978-3-319-17689-5 9.
I D. Ma, L. Yang, R. M. T. Fleming, I. Thiele, B. O. Palsson and M. A. Saunders, “Reliable and efficientsolution of genome-scale models of metabolism and macromolecular expression,” Nature Scientific
Reports, vol. 7, 40863 (2017), https://www.nature.com/articles/srep40863.
16 / 33

Coulomb n-body atomic system simulations
I Alexei Frolov of Queen’s University in Canada has used high-precision arithmetic

to solve a generalized eigenvalue problem that arises in Coulomb n-bodyinteractions.
I Matrices are typically 5,000 x 5,000 or more in size, and are very nearly singular.
I Computations typically involve massive cancellation, and high-precision arithmetic

must be employed to obtain numerically reproducible results.
I These computations typically require 120-digit arithmetic.

Frolov: “We can consider and solve the bound state few-body problems ... beyond ourimagination even four years ago.”

I A. M. Frolov and D. H. Bailey, “Highly accurate evaluation of the few-body auxiliary functions andfour-body integrals,” Journal of Physics B, vol. 36, no. 9 (14 May 2003), pg. 1857–1867.

17 / 33

Numerical solutions in chaotic systems
1 period - TIDES (16 digits)
16 periods -TIDES (300 digits)

45
40
35

z

30
25
20
15
10
5
0
−10

x

0

10

−20

0

y

20

First point TIDES (16 digits)First-Last point TIDES (300 digits)
Last point TIDES (16 digits)

I R. Barrio, M. Rodrı́guez, A. Abad, and F. Blesa, “Breaking the limits: The Taylor series method,” Applied
Mathematics and Computation, vol. 217 (2011), 7940–7954.
18 / 33

Discovering new results in math and physics with very high precision
Methodology:
1. Compute various mathematical entities (limits, infinite series sums, definite
integrals, etc.) to high precision, typically 100–10,000 digits.
2. Use algorithms such as PSLQ to recognize these numerical values in terms ofwell-known mathematical constants.
3. When results are found experimentally, seek formal mathematical proofs of the
discovered relations.
Thousands of results have recently been found using this methodology, both in puremathematics and in mathematical physics.
If mathematics describes an objective world just like physics, there is no reason why
inductive methods should not be applied in mathematics just the same as in physics.
[Kurt Godel]

19 / 33

Examples of experimental mathematics in action
I Discovery of a previously unknown formula for Pi (200 digits):

π=


∞
X
1
n=0

16n

4
2
1
1
−
−
−
.
8n + 1 8n + 4 8n + 5 8n + 6
This formula permits one to calculate binary digits of π beginning at an arbitrary
starting position, without computing any earlier digits. It has deep implications forthe randomness of π, and for the usage of π as a pseudorandom generator.
I Evaluation of Euler sums, used in quantum physics (250–1000 digits).
I Evaluation of Ising integrals, used in quantum physics and string theory(250–2000 digits).
I Evaluation of Poisson polynomials, used in mathematical physics and image
processing (1000–64000 digits).
I D. H. Bailey and J. M. Borwein, “High-precision arithmetic in mathematical physics,” Mathematics, vol.
3 (2015), pg. 337–367.
I D. H. Bailey, J. M. Borwein, J. Kimberley and W. Ladd, “Computer discovery and analysis of large
Poisson polynomials,” Experimental Mathematics, 27 Aug 2016, vol. 26, pg. 349–363.
20 / 33

Hardware support for IEEE 128-bit (quad) floating-point arithmetic?
I Due to the growing need for higher-than-64-bit computing, particularly given therise of very highly parallel systems (thousands or even millions of processors),
hardware support for IEEE 128-bit arithmetic seems inevitable.
I However, at the present time, it is hard to make a business case for providing sucha feature, given its significant hardware cost.
I Thus, for the foreseeable future, those who do such computations will need to rely

on software implementations.
I In the meantime, can anything be done to lower the performance penalty forhigher precision (currently 10–50X performance)?

21 / 33

Quad, double-double and quad-double software
GFORTRAN:I Now includes full support for real(16), i.e., IEEE 128-bit “quad” (34 digits).
I Sadly, IEEE Quad is not yet supported in the GCC and G++ compilers.
I Available at https://gcc.gnu.org/wiki/GFortranBinaries
Binary128 mode in MPFR:I Software implementation of IEEE 128-bit “quad” (34 digits).
I Operates within the MPFR arbitrary precision library.
I Low-level library only — no high-level interface yet.
I Available in MPFR development version 4.0 at http://www.mpfr.org
QD:
I Double-double (31 digits) and quad-double (62 digits) arithmetic.
I Includes high-level interfaces for C++ and Fortran.
I Available at http://crd.lbl.gov/˜dhbailey/mpdist

22 / 33

Arbitrary precision packages
I MPFR: Low-level C library. Very high performance, thread-safe, exact rounding.

Based on GMP. http://www.mpfr.org
I MPFR++: High-level C++ interface to MPFR.

http://www.holoborodko.com/pavel/mpfr/
I MPFUN2015: High-level Fortran interface to MPFR. Thread-safe. Stand-alone

version also available. http://www.davidhbailey.com/mpdist
I ARPREC: C++ package, with high-level interface for C++ and Fortran.

http://crd.lbl.gov/˜dhbailey/mpdist
I Python bigfloat: High-level Python interface to MPFR.

https://pythonhosted.org/bigfloat/
I GPUMP: Low-level CUDA package for GPUs.

https://www.researchgate.net/publication/224175411_GPUMP_A_
Multiple-Precision_Integer_Library_for_GPUs
23 / 33

The attack from the bottom: New applications only require 16 bits
I The applications we have discussed so far all require more than 64-bit arithmetic,

most often quad (128-bit) precision.
I But other applications have recently arisen that are calling for less than IEEE

32-bit arithmetic, in return for greater performance and larger data arrays.
I These applications typically run with IEEE half format (or the similar ARM halfformat), or the new bfloat16 format.
Many of the emerging low-precision applications come from the world of artificial
intelligence and machine learning.

24 / 33

DeepMind’s AlphaGo Zero teaches itself to play Go
I In October 2017, DeepMind researchers programmed an AI

computer with the rules of Go, then had the program playitself — literally teaching itself with no human input.
I After just three days of training, the resulting program“AlphaGo Zero” defeated their earlier program (which had
defeated Ke Jie, the highest rated human) 100 games to 0.
I The Elo rating of Ke Jie, the world’s highest rated Go

player, is 3661. After 40 days of training, AlphaGo Zero’s
Elo rating was over 5000.
I AlphaGo Zero is as far ahead of the world champion as theworld champion is ahead of a good amateur.
I “The latest AI can work things out without being taught,” The Economist, 21 Oct 2017,https://www.economist.com/news/science-and-technology/
21730391-learning-play-go-only-start-latest-ai-can-work-things-out-without.
25 / 33

Artificial intelligence and big data in financeAccording to a Bloomberg report, the finance world isentering a new era:
The fraternity of bond jockeys, derivatives mavens
and stock pickers who’ve long personified theindustry are giving way to algorithms, and soon,
artificial intelligence.
Banks and investment funds have been tinkering for years, prompting anxiety
for employees. Now, firms are rolling out machine-learning software tosuggest bets, set prices and craft hedges. The tools will relieve staff ofroutine tasks and offer an edge to those who stay. But one day, machines
may not need much help.
I S. Kishan, H. Son and M. Rojanasakul, “Robots are coming for these Wall Street jobs,” Bloomberg, 18
Oct 2017, https://www.bloomberg.com/graphics/2017-wall-street-robots/.
26 / 33

Artificial intelligence and big data in finance, continued
I Credit markets: Automating subjective human decisions.
I Foreign exchange: Anticipating variations in demandand prices.
I Equity (stock) markets: Optimizing order execution,

timing purchases and assessing risk.
I Credit (bond) markets: Understanding bond deals, legal

documents and court rulings.
I Macroeconomics: Analyzing satellite images of Chinese industrial sites, farming

regions, Walmart parking lots, etc., to spot trends in the economy. Several hedge
funds and quant firms are already using such AI-based facilities to substantial
advantage.
The message to the finance world is clear: Adapt to AI and machine learning, or die!
27 / 33

Intel’s Nervana neural net software processor uses bfloat16 format
I Intel has announced a new line of processors, known as “Nervana

Neural Net L-1000,” targeted directly at AI applications.
I This processor supports the bfloat16 numerical format: one sign bit,

8 exponent bits and only 7 mantissa bits plus one hidden bit(roughly 2-3 digits).
Intel’s Naveen Rao, head of AI research, says:
Over time, Intel will be extending bfloat16 support across our AI productlines, including Intel Xeon processors and Intel FPGAs. This is part of acohesive and comprehensive strategy to bring leading AI training capabilities
to our silicon portfolio.
I Naveen Rao, “Behyond the CPU or GPU: Why enterprise-scale artificial intelligence requires a more
holistic approach,” Intel Newsroom, 23 May 2018, https://newsroom.intel.com/editorials/artificial-intelligence-requires-holistic-approach/.
28 / 33

Google has also adopted bfloat16 for AI applications
I Google researchers have adopted bfloat16 for theirTensorFlow Distributions library, which is a software
framework to support AI applications.
I Their library is designed to provide fast, numericallystable methods for generating samples and computing
statistics.

I J. V. Dillon, I. Langmore, D. Tran, E. Brevdo, S. Vasudevan, D. Moore, B. Patton, A. Alemi, M. Hoffman
and R. A. Saurous, “TensorFlow Distributions,” https://arxiv.org/pdf/1711.10604.pdf.

29 / 33

Challenges for low-precision computing
I Reproducibility: Only eight mantissa bits severely limits any measure of accuracy

and reliability. At the least, this means that software must provide a simple means
to re-run an application using 32-bit or 64-bit precision to verify some level ofreliability and reproducibility.
I Software availability: So far there are only a handful of specialized packages

available for 16-bit computing. Well-supported high-level language interfaces areconspicuously absent.
I Hardware compatibility: The proliferation of hardware formats (IEEE half,bfloat16) complicates the development of software.

30 / 33

Software support needed for variable precision
1. High-level language support. Manually calling a subroutine library is much too
tedious and error-prone. To be usable for real-world applications, precision
software must include fully integrated language support, so one can write, forinstance,d = a + b * sqrt (c + 3.)
where a, b, c and d are of possibly different precision levels, and all conversions
and computations are automatically performed. Most precision software does notinclude high-level language support.
2. Thread-safe: Variable precision software must be thread-safe, so it can support
shared-memory parallel applications. Most precision software is not thread-safe.
3. Multiple-level support: Multiple precision levels are often required in applications,since even if some sections of code can use limited precision, other sections
require more. Most precision software does not provide more than one level.

31 / 33

Software support needed for variable precision, continued
4. A way to specify that mixed-precision computations in an expression are to bedone to the maximum level of precision.
In this line of code (C or Fortran):b = a + 1./3.
the division 1./3. is done with only single or double precision, even if a and b are
declared to be of a higher-precision datatype (a similar issue with half precision).
Most languages, compilers and precision software packages do not deal with thisissue, resulting in hard-to-find errors.
Partial solution: Provide a full set of conversion routines, so that one can writeb = a + real (1.,16) / real (3.,16)
where “16” (bytes) is the type specifier for quad precision, for instance.
But it is also important that precision errors such as “b = a + 1./3.” beautomatically detected at run time. (This is done in the MPFUN2015 software.)32 / 33

SummaryI A growing number of applications (mostly in the high-performance scientificarena) require more than the IEEE 64-bit arithmetic.
I A growing number of applications (mostly in the artificial intelligence / machinelearning arena) require less than IEEE 32-bit arithmetic.
The future trend is clear: Variable precision calculations, ranging from very low
precision (e.g., bfloat16) to very high precision (128-bit or more) are emerging as animportant fixture of both business and scientific/engineering applications. Ignore it atyour peril!

This talk is available at
http://www.davidhbailey.com/dhbtalks/dhb-camb-2018.pdf
33 / 33