Variable precision computing: Applications and challenges
David H. Bailey
Lawrence Berkeley Natlional Laboratory (retired)
University of California, Davis, Department of Computer Science
This talk is available at:
http://www.davidhbailey.com/dhbtalks/dhb-icerm-2020.pdf

1 / 29

Variable precision computing

I Mixed 16-bit/32-bit: Artificial intelligence, machine learning, graphics.
I Mixed 32-bit/64-bit: A broad range of current scientific applications.
I Mixed 64-bit/128-bit: Large applications with numerically sensitive portions.
I Mixed 64-bit/multiple: Applications in computational mathematics and physics.

2 / 29

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)(none)

Nickname
IEEE half
ARM half
bfloat16
IEEE single
EEE double
IEEE extended
IEEE quad
double-double
quad-double
double-quad
multiple

Sign
1
1
1
1
1
1
1
1
1
1
variesNumber of bits
Exponent Mantissa
5
10
5
10
8
7
7
24
11
52
15
64
15
112
11
104
11
208
15
224
variesvaries

Hidden
1
1
1
1
1
0
1
2
4
1
variesDigits
3
3
2
7
15
19
34
31
62
68
varies3 / 29Advantages of variable precision
Compared with using a fixed level of precision for the entire application, employing avariable precision framework features:
I Faster processing.
I Better cache utilization.
I Lower run-time memory usage.
I Lower offline data storage.
I Lower energy costs.
I Improved reproducibility and replicability.

4 / 29

Numerical reproducibility
Excerpt for 2012 ICERM workshop on reproducibility in mathematical and scientificcomputing:Numerical reproducibility has emerged as a particularly important issue, since the scaleof computations has greatly increased in recent years, particularly with computations
performed on many thousands of processors and involving similarly large datasets.
Large computations often greatly magnify the level of numeric error, so that numericaldifficulties that were once of little import now are large enough to alter the course of
the computation or to draw into question the overall validity of the 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,” manuscript, 2 Feb 2013,https://www.davidhbailey.com/dhbpapers/icerm-report.pdf.
5 / 29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.
In an attempt to speed up the calculation, researchers found that merely changing theunderlying math library (which should only affect at most the last bit) resulted in somecollisions being missed or misidentified.
Questions:
I How extensive are these numerical difficulties?
I How can numerically sensitive code be tracked down?
I How can a large code library such as this be maintained, producing numericallyreliable and reproducible results?
6 / 29

How to solve accuracy problems and ensure reproducibility
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 an expert numerical analyst to convert large portions of the code to useinterval arithmetic (greatly increasing run time and code complexity).
4. Employ significantly higher precision than is really necessary (greatly increasing
run time).
5. Employ variable-precision arithmetic, assisted with some “smart” tools to help
determine where extra precision is needed and where it is not.
Item 5 is the only practical solution for real-world scientific computing.

7 / 29

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 very likely over 1000; probably 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 about 2% of U.C. Berkeley graduates who will do technicalcomputing in their careers have had rigorous, expert-level training in numerical analysis. 8 / 29Mixed half-single applications: Machine learning
NVIDIA researchers found that a combination of IEEE 32-bit and bfloat16 arithmetic
achieved nearly the same training loss reduction curve on a English language learningmodel as with 100% IEEE 32-bit:
Maintain a master copy of training system weights in IEEE 32-bit, and select a scaling
factor S. Then for each iteration:
I Convert the array of 32-bit weights to a bfloat16 array.
I Perform forward propagation with bfloat16 weights and activations.
I Multiply the resulting loss with a scaling factor S.
I Perform backward propagation with bfloat16 weights, activations and gradients.
I Multiply the weight gradient by 1/S.
I Update the 32-bit weights using the bfloat16 data.

I “Deep learning SDK documentation,” NVIDIA, 2019,https://docs.nvidia.com/deeplearning/sdk/index.html.
9 / 29DeepMind’s AlphaGo Zero teaches itself to play Go
I In October 2017, DeepMind researchers programmed amachine learning system with the rules of Go, then had theprogram play itself — 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 was as far ahead of the world champion as

the world 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.
10 / 29

Challenges for mixed half-single computing
I Reproducibility: Only eight mantissa bits severely limits any measure of accuracy

and reliability. At the least, software must provide a simple means to re-run anapplication using 32-bit or 64-bit precision to verify some level of reliability andreproducibility.
I Software availability: So far there are relatively few specialized packages andwell-supported high-level language interfaces for 16-bit computing.
I Hardware compatibility: The proliferation of hardware formats (IEEE half,bfloat16, arm16) complicates the development of software.

11 / 29Mixed single-double application: LU decomposition
I The usual scheme is to factor a pivoted coefficient matrix PA = LU, where L islower triangular and U is upper triangular, then solve Lu = Pb and Ux = y , allusing double precision.
I In a mixed single-double scheme, the factorization and solutions are done in single
precision; then the residual and full solutions are done in double.
I Note that the factorization of A, which is the only step that is O(n3 ), is done
entirely in single.
I A similar scheme can be used for sparse systems.
I A broad range of iterative refinement methods can be performed in this way.
I Baboulin et al. (below) achieved speedups of 1.86 on a direct sparse system on asuite of test problems.

I M. Baboulin, A. Buttari, J. Dongarra, J. Kurzak, J. Langou, P. Luszczek and S. Tomov, “Accelerating
scientific computations with mixed precision algorithms,” Computational Physics Communications, vol.
180 (2009), 2256–2253.
12 / 29Mixed double-quad application: Computational biology
I Certain models of biochemical reaction networks, e.g., “metabolic expression

models,” involve solving large multiscale linear programming problems.
I Since double precision implementations of these models often produce numericallyunreliable results, researchers have resorted to exact rational arithmetic methods,
which are enormously expensive — run times are typically weeks or months.
I Ding Ma and Michael Sanders of Stanford have achieved excellent results using
double codes enhanced with quad precision for numerically sensitive operations.
I Their code was a straightforward modification of a double precision Fortran code,using the REAL(16) datatype, which is supported in the GNU gfortran compiler.

I D. Ma and M. Saunders, “Solving multiscale linear programs using the simplex method in quadruple
precision,” in M. Al-Baali, L. Grandinetti and A. Purnama, ed., Recent Developments in Numerical
Analysis and Optimization, Springer, NY, 2017,
http://web.stanford.edu/group/SOL/reports/quadLP3.pdf.
13 / 29

Innocuous example where IEEE 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.
Common approach: Use a least-squares scheme, solving the linear system:






P n

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 xkk=1 xk yk    ..  = .
..
..
 .  Pn . nPn . 2nan
k=1 xk yk
k=1 xkAn implementation using IEEE 64-bit arithmetic fails to find the correct polynomial,
even if one rounds results to nearest integer.
But an implementation using 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 .
14 / 29

Innocuous example, cont.
The result on the previous page can be obtained with 64-bit arithmetic using Lagrange
interpolation or the Demmel-Koev algorithm. But few outside of expert numericalanalysts are aware of these schemes.
Besides, even these two 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 .
This is easily found by an IEEE quad or double-double implementation of any of thethree algorithms (least-squares, Lagrange or Demmel-Koev).

15 / 29Mixed double/multiprecision applications
1. Solutions of ordinary differential equations in planetary orbit models (32 digits).
2. Supernova simulations (32–67 digits).
3. Scattering amplitudes of fundamental particles (32 digits).
4. Solutions of certain discrete dynamical systems (32 digits).
5. Coulomb n-body atomic system simulations (32–120 digits).
6. Solving ordinary differential equations using the Taylor algorithm (100–300 digits).
7. Recognizing integrals from the Ising theory of mathematical physics in terms of
simple mathematical identities (200–1000 digits).
8. Finding minimal polynomials connected to the Poisson equation of mathematicalphysics (1000–64,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), 10106–10121.
16 / 29Mixed double/multiprecision: Computational math and physics
Frequently-used “experimental math” methodology:
1. Compute a limit, infinite series sum, definite integral or other mathematicalquantity to high precision, typically 100–10,000 digits.
2. Use an integer relation algorithm such as “PSLQ” to recognize the numericalvalue in terms of a mathematical expression. In other words, if the computed
quantity is α, and if one suspects that α = x0 is given by a formula involving
known constants (x1 , x2 , · · · , xn−1 ), then use PSLQ to find integers(a0 , a1 , a2 , · · · , an−1 ) such that a0 x0 + a1 x1 + · · · + an−1 xn−1 = 0.
3. When a result is found experimentally, seek a formal mathematical proof of the
discovered relation.
4. Efficient implementations typically require three levels of precision, such as
double/250-digit/10,000-digit, varying up and down many times within the run,performing as much computation as possible using only double precision.

17 / 29

Experimental mathematics and physics results
I New formulas and identities for various mathematical constants, including

π, log(2) and others, some with remarkable properties (250–1000 digits).
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–64,000 digits).

I D. H. Bailey and J. M. Borwein, “High-precision arithmetic in mathematical physics,” Mathematics, vol.
3 (2015), 337–367.
18 / 29Mixed double/multiprecision application: Finding Poisson polynomials
In 2012 Richard Crandall of Apple Computers found that each pixel in the image
produced by a sharpening algorithm is given by the 2-D Poisson potential function ofmathematical physics:
1 X cos(mπx ) cos(nπy )φ2 (x , y ) = 2
π m,n odd
m2 + n2
Crandall and others then found that for rational x and y ,1φ2 (x , y ) = log απwhere α is algebraic, i.e., the root of a some integer polynomial of degree m.
What is the connection between x and y and the degree of α?
I D. H. Bailey, J. M. Borwein, R. E. Crandall and J. Zucker, “Lattice sums arising from the Poisson
equation,” Journal of Physics A: Mathematical and Theoretical, vol. 46 (2013), 115201.
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, 349–363.
19 / 29

Finding Poisson polynomials, cont.
By computing high-precision numerical values of φ2 (x , y ) for various specific rational x
and y , using three precision levels (e.g., double/250-digit/10,000-digit), and applying avariant of the PSLQ algorithm, we found the following results among others:
s
5
6
7
8
9

10

Minimal polynomial of α corresponding to x = y = 1/s
1 + 52α − 26α2 − 12α3 + α4
1 − 28α + 6α2 − 28α3 + α4
−1 − 196α + 1302α2 − 14756α3 + 15673α4 + 42168α5 − 111916α6 + 82264α7−35231α8 + 19852α9 − 2954α10 − 308α11 + 7α12
1 − 88α + 92α2 − 872α3 + 1990α4 − 872α5 + 92α6 − 88α7 + α8
−1 − 534α + 10923α2 − 342864α3 + 2304684α4 − 7820712α5 + 13729068α6−22321584α7 + 39775986α8 − 44431044α9 + 19899882α10 + 3546576α11−8458020α12 + 4009176α13 − 273348α14 + 121392α15−11385α16 − 342α17 + 3α18
1 − 216α + 860α2 − 744α3 + 454α4 − 744α5 + 860α6 − 216α7 + α8

20 / 29

Finding Poisson polynomials, cont.
Based on these preliminary results, Jason Kimberley conjectured that the degree m(s)of the minimal polynomial associated with the case x = y = 1/s is given by this rule:Set m(2) = 1/2. Otherwise for primes p congruent to 1 mod 4, setm(p) = int2 (p/2), where int denotes greatest integer, and for primes pcongruent to 3 mod 4, set m(p) = int (p/2)(int (p/2) + 1). Then for any
other positive integer s whose prime factorization is s = p1e1 p2e2 · · · prer ,r −1m(s) = 4rY
2(ei −1)pim(pi ).

i=1

Additional three-level precision computations, (e.g., double/2000-digit/64,000-digit),
confirmed Kimberley’s formula for all s up to 40, and for even s up to 52.
By examining the computed results, connections were found to a sequence ofpolynomials defined in a 2010 paper by Savin and Quarfoot. These observations
ultimately led to a proof of Kimberley’s formula and several other facts.
21 / 29

Software and language issues: Conversions between types
Existing programming languages have significant issues with varying precision levels,
particularly in conversions between types. For example, consider:program testprec
real(8) d1, d2, d3
d1 = sqrt (3. + 1./7.)
d2 = sqrt (3.d0 + 1./7.)
d3 = sqrt (3. + 1.d0/7.)
write (6, ’(f20.15)’) d1, d2, d3
stop
end

This code, when compiled with the gfortran compiler with default options, produces:
1.772810459136963
1.772810522656991
1.772810520855837

Only the last is a correct, full-precision value!
One could set a compiler flag to treat all single constants as double (default in somelanguages), but this defeats the goal of using varying precision levels.

22 / 29

Software and language issues: Conversions between types, cont.
Similar difficulties, greatly magnified, afflict mixed double-multiprecision applications:program testprec
use mpmodule
type (mp real) r1, r2, r3
r1 = sqrt (3.d0 + 1.d0/7.d0)
r2 = sqrt (mpreal (3.d0, 10) + 1.d0/7.d0)
r3 = sqrt (3.d0 + mpreal (1.d0, 10)/7.d0)
call mpwrite (6, 80, 60, r1, r2, r3)
stop
end

DHB’s MPFUN-MPFR package checks for conversion problems like this at runtime,
producing this error message:*** MP CHECKDP: DP value has more than 40 significant bits:
1.772810520855837D+00
and thus very likely represents an unintended loss of accuracy.

One solution: A compiler option to outlaw all mixed-mode operations and automatic
type conversions. Unpopular? Perhaps. But finding precision problems is even harder.
23 / 29

Thread safety in high precision libraries
Most multiprecision libraries in use today are not thread-safe. Thus applications
written to use them do not work in a shared-memory parallel environment.
The solution is to include the following in each multiprecision datum:
I The full length of the storage allocated to this datum.
I The current working precision associated with this datum.
I The length of nonzero data associated with this datum.

Currently the MPFR package (see below) is the only low-level multiprecision librarythat is thread-safe (provided thread-safe option is set); it is also the fastest.
DHB’s MPFUN-MPFR high-level Fortran multiprecision package is built on MPFR.

I L. Fousse, G. Hanrot, V. Lefevre, P. Pelissier and P. Zimmermann, “MPFR: A multiple precision binary
floating-point library with correct rounding,” ACM Transactions on Mathematical Software, vol. 33(2007), https://doi.org/10.1145/1236463.1236468.
24 / 29

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).
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.
25 / 29

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 an addition to the IEEE-754 standard
that would permit guaranteed bit-for-bit replicability, not only on a single systembut also on a parallel system — see references below.
I ARM (a U.K. computer processor vendor) has experimented with a high-precision
accumulator, using an adaptable vector unit.
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 users to even be aware that they have a problem
with numerical sensitivity and/or inaccuracy?
I P. Ahrens, J. Demmel and H. D. Nguyen, “Algorithms for efficient reproducible floating pointsummation,” manuscript available from authors.
I J. Riedy and J. Demmel, “Augmented arithmetic operations proposed for IEEE-754 2018,”
http://www.ecs.umass.edu/arith-2018/pdf/arith25_34.pdf.

26 / 29

New approaches to floating-point arithmetic
I Although the IEEE-754 floating-point standard has served us very well for several

decades, some wonder if it time to completely rethink floating-point arithmetic.
I Gustafson and Yonemoto have proposed “Unum,” a floating-point standard that is

flexible in size and precision — see references below.
I Lindstrom summarizes Unums and some other proposals — see below.

I J. L. Gustafson, The End of Error: Unum Computing, Chapman and Hall, New York, 2017.
I J. Gustafson and I. Yonemoto, “Beating floating point at its own game: Posit arithmetic,” 2017,
http://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf.
I P. Lindstrom, “Fixed-rate compressed floating-point arrays,” IEEE Transactions on Visualization andComputer Graphics, vol. 20 (2014), 2674–2683, https://doi.org/10.1109/TVCG.2014.2346458.
27 / 29

Research questions to be explored
1. How are various scientific disciplines currently using variable precision? What newopportunities are on the horizon?
2. What does the explosive growth in artificial intelligence and machine learning mean for
hardware, software and tools to support variable precision? What are the fundamental
accuracy requirements of such computations?
3. Are there new mathematical approaches to efficiently and accurately compute various
mathematical operations and library functions to both very modest precision (2-3 digit
accuracy) and very high precision (hundreds or thousands of digits)?
4. What new formats and representation systems to facilitate variable level precisioncomputing are available? To what extent do these new proposals overcome strengths andweaknesses in existing systems such as IEEE-754 hardware and MPFR software?
5. What new support is needed from computer hardware and software vendors for variable
precision computing? Could some modest hardware changes be made to support low,
variable precision and high precision, rather than relying on all-software implementations?
6. Almost all of published research in numerical mathematics and numerical analysis to datehas focused on either 32-bit or 64-bit results. Are there other techniques that may be
28 / 29superior in a variable precision environment?

Research questions to be explored, cont.
6. Classical error bounds used in numerical analysis, such as those developed by Wilkinson,
typically assume worst-case scenarios and thus are not very useful in a limited precision orvariable precision environment. Can more realistic error estimates be derived?
7. How can iterative refinement best be performed in a variable precision environment?
What new opportunities are there for using iterative refinement with variable precision?
8. What algorithms and software techniques can effectively compress floating-point data,on-the-fly, to save data storage and transfer? Can improved schemes be devised?
9. What software tools are available to automatically or semi-automatically analyze an
application code, identify numerically sensitive spots and propose or implement remedies?
What other tools are needed?
10. How can we ensure that computations with variable precision levels are reproducible?
What does reproducibility even mean when computing, say, with only 16-bit precision?

This talk is available at http://www.davidhbailey.com/dhbtalks/dhb-icerm-2020.pdf.
29 / 29