Using High-Precision Arithmetic to Conquer Numerical Error
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Laboratory (retired)Computer Science Department, University of California, Davis

14 May 2014

Numerical reproducibility in scientific computing
A December 2012 workshop on reproducibility in 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.

Numerical reliability
Many applications routinely use either 32-bit or 64-bit IEEE arithmetic, and employ
fairly simple algorithms, assuming that all is well. But problems can arise:
1. Large-scale, highly parallel simulations, running on systems with hundreds ofthousands or millions of processors — numerical sensitivities are greatly magnified.
2. Certain applications with highly ill-conditioned linear systems.
3. Large summations, especially those involving +/- terms and cancellations.
4. Long-time, iterative simulations (such as molecular dynamics or climate models).
5. Computations to resolve small-scale phenomena.
6. Studies in computational physics or experimental mathematics often require hugeprecision levels.
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.

Analysis of collisions at the Large Hadron Collider
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: 5 millions line of C++ and Python code, developed by roughly 2000

physicists and engineers over 15 years.
I 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 serious are these numerical difficulties?
I How can they be tracked down?
I How can the library be maintained, producing numerically reliable results?

Ways to enhance numerical reproducibility and solve accuracy problems
1. Employ an expert numerical analyst to re-examine every algorithm employed in thecomputation to ensure that only the most stable known schemes are being used.
2. Carefully analyze every step of the computation to ascertain the level of numericalsensitivity at each step.
3. Employ interval arithmetic for large portions of the application (which greatly
increases run time and code complexity).
4. Employ higher-precision arithmetic arithmetic, assisted with some “smart” tools
to help determine where extra precision is needed and where it is not.
Item #4 is often the only practical solution.

Numerical analysis expertise among U.C. Berkeley graduates
Of the 2010 U.C. Berkeley graduating class, 870 were in disciplines likely to require
technical computing:
I Division of Mathematical and Physical Sciences (Math, Physics, Statistics).
I College of Chemistry.
I College of Engineering (including Computer Science).

Other fields (not counted) that will likely involve significant computing:
I Biology, geology, medicine, economics, psychology, sociology.

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: Fewer than 2% of Berkeley graduates who will do technical computing
have had rigorous training in numerical analysis!

Enhancing reproducibility with selective usage of high-precision arithmetic
Problem: Findof the irregular function
P the arc length−kg (x) = x + 0≤k≤10 2 sin(2k x), over the interval (0, π)(using 107 abscissa points).
3.0

I If this computation is done with ordinary double

precision arithmetic, the calculation takes 2.59 seconds
and yields the result 7.073157029008510.
I If it is done using all double-double arithmetic (31-digit

2.5

2.0

1.5

1.0

0.5

accuracy), it takes 47.39 seconds seconds and yields theresult 7.073157029007832 (correct to 15 digits).
I But if only the summation is changed to double-double,

the result is identical to the double-double result (to 15
digits), yet the computation only takes 3.47 seconds.

0.5

1.0

1.5

2.0

2.5

3.0

Graph
P of g (x) =x+ 0≤k≤10 2−k sin(2k x),
over (0, π).

Aren’t 64 bits enough?
There has been considerable resistance in the scientific computing community to thenotion that more than 64-bit arithmetic is not only useful, but may even be essential insome scientific computations. Why?
I Many are persuaded that physical reality fundamentally does not require highprecision, beyond, say, the limits of 64-bit IEEE arithmetic.
I Many believe that numerical sensitivity problems are always due to the usage ofinferior algorithms.
I Some regard the usage of high-precision arithmetic as “cheating” or “sinful”(because it is too easy?).
I Until a few years ago, easy-to-use high-precision arithmetic software facilities were

not widely available to experiment with, or the computational cost was too great.
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.

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..
P .

nn
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    ..  = ..
..
 .  

.
.
P
Pn
n
n
2n
a
x
y
x
n
k=1 k kk=1 k

A 64-bit computation (e.g., using Matlab, Linpack or LAPACK) fails to find thecorrect polynomial in this instance, even if one rounds results to nearest integer.
However, if Linpack routines are converted to use double-double arithmetic (31-digit
accuracy), the above computation quickly produces the correct polynomial:f (x) = 1 + 1048577x 4 + x 8 = 1 + (220 + 1)x 4 + x 8

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 scientists, outside of expert
numerical analysts, are aware of these schemes.
Besides, even these schemes fail for higher-degree problems.
For example: (1, 134217731, 8589938753, 97845255883, 549772595201,
2097396156251, 6264239146561, 15804422886323, 35253091827713,
71611233653971, 135217729000001, 240913322581691, 409688091758593)
is generated by:
f (x) = 1 + 134217729x 6 + x 12 = 1 + (227 + 1)x 6 + x 12
Lagrange interpolation and Demmel-Koev fail for this problem, but a straightforward
Linpack scheme, implemented with double-double arithmetic, works fine.

Free software for high-precision computation
1. ARPREC. Arbitrary precision, with numerous algebraic and transcendental
functions. High-level interfaces for C++ and Fortran-90.
http://crd.lbl.gov/~dhbailey/mpdist.
2. GMP. Produced by a volunteer effort and distributed under the GNU license.
http://gmplib.org.
3. MPFR. C library for multiple-precision floating-point computations with exact
rounding, based on GMP. http://www.mpfr.org.
4. MPFR++. High-level C++ interface to MPFR.
http://perso.ens-lyon.fr/nathalie.revol/software.html.
5. MPFUN90. Similar to ARPREC, but is written entirely in Fortran-90 and provides
only a Fortran-90 interface. http://crd.lbl.gov/~dhbailey/mpdist.
6. QD. Performs “double-double” (31 digits) and “quad-double” (62 digits)
arithmetic. High-level interfaces for C++ and Fortran-90.
http://crd.lbl.gov/~dhbailey/mpdist.

U.C. Berkeley’s CORVETTE project and the “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,” manuscript, May 2013.

Precimonious in action for the arc length computation
Original code:

Automatically tuned code:

int main() {int i, n = 1000000;
long double h, t1, t2, dppi;
long double s1;...
for (i = 1; i <= n; i++) {t2 = fun (i * h);s1 = s1 + sqrt (h*h + (t1-t2)*(t1-t2));
t1 = t2;}
return 0;}

int main() {int i, n = 1000000;
double h, t1, t2, dppi;
long double s1;...
for (i = 1; i <= n; i++) {t2 = fun (i * h);s1 = s1 + sqrt (h*h + (t1-t2)*(t1-t2));
t1 = t2;}
return 0;}

The tuned code produces the same answers, to 15-digit accuracy, yet runs faster.

Some applications where high precision is useful or essential
1. Planetary orbit calculations (32 digits).
2. Supernova simulations (32–64 digits).
3. Climate modeling (32 digits).
4. Coulomb n-body atomic system simulations (32–120 digits).
5. Schrodinger solutions for lithium and helium atoms (32 digits).
6. Electromagnetic scattering theory (32–100 digits).
7. Scattering amplitudes of fundamental particles (32 digits).
8. Discrete dynamical systems (32 digits).
9. Theory of nonlinear oscillators (64 digits).
10. The Taylor algorithm for ODEs (100–600 digits).
11. Ising integrals from mathematical physics (100–1000 digits).
12. Problems in experimental mathematics (100–50,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.

Long-term planetary orbit calculations
Researchers have recognized for centuries that planetary orbits exhibit chaotic behavior:“The orbit of any one planet depends on the combined motions of all the planets,
not to mention the actions of all these on each other. To consider simultaneously allthese causes of motion and to define these motions by exact laws allowing of
convenient calculation exceeds, unless I am mistaken, the forces of the entire human
intellect.”[Isaac Newton, Principia, 1687]

Long-term simulations of planetary orbits using double precision do fairly well for long
periods, but then fail at certain key junctures.
Researchers have found that double-double or quad-double arithmetic is occasionallyrequired to avoid severe inaccuracies, even if other techniques are employed to reduce
numerical error.
I G. Lake, T. Quinn and D. C. Richardson, “From Sir Isaac to the Sloan survey: Calculating the structure
and chaos due to gravity in the universe,” Proc. of the 8th ACM-SIAM Symp. on Discrete Algorithms,
1997, pg. 1–10.

Supernova simulations
I Researchers at LBNL have used quad-double arithmetic

to solve for non-local thermodynamic equilibrium
populations of iron and other atoms in the atmospheresof supernovas.
I Iron may exist in several species, so it is necessary tosolve for all species simultaneously.
I Since the relative population of any state from thedominant state is proportional to the exponential of theionization energy, the dynamic range of these values canbe very large, and cancellations may occur.
I Quad-double arithmetic (62 digits) was employed toreduce these errors.

P. H. Hauschildt and E. Baron,“The numerical solution of the
expanding stellar atmosphereproblem,” Journal Computationaland Applied Mathematics, vol.
109 (1999), pg. 41–63.

Climate modeling: High-precision for reproducibility
I Climate and weather simulations are fundamentallychaotic: if microscopic changes are made to the current
state, soon the future state is quite different.
I In practice, computational results are altered even if

minor changes are made to the code or the system.
I This numerical variation is a major nuisance for code

maintenance.
I He and Ding found that by using double-double

arithmetic in two key inner loops, most of thisnumerical variation disappeared.

Y. He and C. Ding, “Using
accurate arithmetics to improve
numerical reproducibility andstability in parallel applications,”
Journal of Supercomputing, vol.
18, no. 3 (Mar 2001), pg.
259–277.

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 and are very nearly singular.
I Computations typically involve massive cancellation, and high-precision arithmetic

must be employed to obtain numerically reproducible results.
I Frolov has also computed elements of the Hamiltonian matrix and the overlap

matrix in four- and five-body systems.
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 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.

Taylor’s method for ODEs with high-precision arithmetic
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)

Numerical integration of the L25-R25 unstable orbit for the Lorenz model during 16
time periods using the TIDES code with 300 digits, and 1 time period using 15 digits.

Experimental mathematics:
Discovering new mathematical results by computer
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.
Many 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 inphysics.” – Kurt Godel

The PSLQ integer relation algorithm
Let (xn ) be a given vector of real numbers. An integer relation algorithm either finds
integers (an ) such thata1 x1 + a2 x2 + · · · + an xn = 0(to within the “epsilon” of the arithmetic being used), or else finds bounds within
which no relation can exist.
The “PSLQ” algorithm of mathematician-sculptor Helaman Ferguson is the mostwidely used integer relation algorithm.
Integer relation detection requires very high precision (at least n × d digits, where d isthe size in digits of the largest ak ), both in the input data and in the operation of thealgorithm.
I D. H. Bailey and D. J. Broadhurst, “Parallel integer relation detection: Techniques and applications,”
Mathematics of Computation, vol. 70, no. 236 (Oct 2000), pg. 1719–1736.

PSLQ, continued
I PSLQ constructs a sequence of integer-valued matrices Bn that reduce the vectory = x · Bn , until either the relation is found (as one of the columns of matrix Bn ),or else precision is exhausted.
I A relation is detected when the size of smallest entry of the y vector suddenly
drops to roughly “epsilon” (i.e. 10−p , where p is the number of digits ofprecision).
I The size of this drop can be viewed as a “confidence level” that the relation is nota numerical artifact: a drop of 20+ orders of magnitude almost always indicates areal relation.
Efficient variants of PSLQ:
I 2-level and 3-level PSLQ perform almost all iterations with only double precision,
updating full-precision arrays as needed. They are hundreds of times faster thanthe original PSLQ.
I Multi-pair PSLQ dramatically reduces the number of iterations required. It was
designed for parallel systems, but runs faster even on 1 CPU.

Decrease of log 10(min |yi |) in multipair PSLQ run

Simple application of PSLQ: Find minimal polynomial of algebraic number
Suppose we are given an algebraic number α, i.e., α is a root of an algebraic equation
with integer coefficients: a0 + a1 x + a2 x 2 + · · · + an x n = 0. How can we reconstruct
the minimal polynomial of α, knowing just the degree n (or some estimate of n) andthe numerical value of α?
Solution: Compute the vector x = (1, α, α2 , · · · , αn ), and then apply PSLQ to this(n + 1)-long vector.
Example: Find the minimal polynomial of the constant, assuming that it of degree nomore than 10:
87.055259441185777057425040060690481806178409055793 . . .
Answer:0 = 1 − 88α + 92α2 − 872α3 + 1990α4 − 872α5 + 92α6 − 88α7 + α8The first major PSLQ discovery: The BBP formula for π
In 1996, a PSLQ program discovered this new formula for π:


∞
X
1
4
2
1
1
π =
−
−
−
16k 8k + 1 8k + 4 8k + 5 8k + 6
n=0This formula permits one to compute binary (or hexadecimal) digits of π beginning at
an arbitrary starting position, using a very simple scheme that requires only standard
64-bit or 128-bit arithmetic.
In 2004, Borwein, Galway and Borwein proved that no base-n formulas of this type
exist for π, except when n = 2m .
BBP-type formulas (discovered with PSLQ) are now known for numerous othermathematical constants.
1. D. H. Bailey, P. B. Borwein and S. Plouffe, “On the rapid computation of various polylogarithmic
constants,” Mathematics of Computation, vol. 66, no. 218 (Apr 1997), pg. 903–913.
2. J. M. Borwein, W. F. Galway and D. Borwein, “Finding and excluding b-ary Machin-type BBP formulae,”
Canadian Journal of Mathematics, vol. 56 (2004), pg. 1339–1342.

Some other BBP-type formulas found using PSLQ
∞

π2

=

π2

=

ζ(3)

=


216
72
54
9
144
−
−
−
+
(6k + 1)2
(6k + 2)2
(6k + 3)2
(6k + 4)2
(6k + 5)2
k=0

∞
2 X 1
405
81
27
243
−
−
−
27
729k (12k + 1)2(12k + 2)2(12k + 4)2(27k + 5)2
k=0

72
9
9
5
1
−
−
−
−
+
(12k + 6)2(12k + 7)2(12k + 8)2(12k + 10)2
(12k + 11)2

∞
1 X 1
6144
43008
24576
30720
1536
−
+
+
−
1792
212k (24k + 1)3(24k + 2)3(24k + 3)3(24k + 4)3(24k + 5)3
1X 1
8
64k



k=0

3072
768
3072
2688
192
1536
+
+
−
−
−
−
(24k + 6)3(24k + 7)3(24k + 9)3(24k + 10)3
(24k + 11)3
(24k + 12)3
96
672
384
24
48
12
−
−
−
+
+
−
(24k + 13)3
(24k + 14)3
(24k + 15)3
(24k + 17)3
(24k + 18)3
(24k + 19)3

120
48
42
3
+
+
−
+
(24k + 20)3
(24k + 21)3
(24k + 22)3
(24k + 23)3
I D. H. Bailey, J. M. Borwein, A. Mattingly and G. Wightwick, “The computation of previously inaccessible
digits of Pi2 and Catalan’s constant,” Notices of the AMS, vol. 60 (2013), no. 7, pg. 844-854.

High-precision tanh-sinh numerical integration
Given f (x) defined on (−1, 1), define g (t) = tanh(π/2 sinh t). Then setting x = g (t)yieldsZ 1−1

f (x) dx

=

Z ∞

−∞

f (g (t))g 0 (t) dt ≈ hNX

wj f (xj ),j=−Nwhere xj = g (hj ) and wj = g 0 (hj ). Since g 0 (t) goes to zero very rapidly for large t, theproduct f (g (t))g 0 (t) typically is a nice bell-shaped function, so that the simple
summation above converges very rapidly. Reducing h by half typically doubles thenumber of correct digits.
We have found that tanh-sinh is the best general-purpose integration scheme forfunctions with vertical derivatives or singularities at endpoints, or for any function atvery high precision (> 1000 digits). Otherwise we use Gaussian quadrature.
I D. H. Bailey, X. S. Li and K. Jeyabalan, “A comparison of three high-precision quadrature schemes,”
Experimental Mathematics, vol. 14 (2005), no. 3, pg. 317–329.

Ising integrals from mathematical physics
We recently applied our methods to study three classes of integrals (one of which wasreferred to us by Craig Tracy of U.C. Davis) that arise in the Ising theory ofmathematical physics:
Z
Z ∞
1
dun
4 ∞
du1
···
···
Cn :=
P
2
n! 0
u
un
n
1
0
j=1 (uj + 1/uj )Q  ui −uj 2
Z ∞Z ∞i<j ui +uj
4
dun
du1
Dn :=
···
···
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )

2
Z 1
Z 1
Y uk − uj dt2 dt3 · · · dtn
En = 2···
uk + uj
0
0
1≤j<k≤nwhere in the last line uk = t1 t2 · · · tk .

I D. H. Bailey, J. M. Borwein and R. E. Crandall, “Integrals of the Ising class,” Journal of Physics A:
Mathematical and General, vol. 39 (2006), pg. 12271–12302.

Limiting value of Cn : What is this number?
Key observation: The Cn integrals can be converted to one- dimensional integrals
involving the modified Bessel function K0 (t):Z
2n ∞ nCn =tK0 (t) dt
n! 0
High-precision numerical values, computed using this formula and tanh-sinhquadrature, approach a limit. For example:
C1024 = 0.6304735033743867961220401927108789043545870787 . . .
What is this number? We copied the first 50 digits into the online Inverse Symbolic
Calculator (ISC) at http://carma-lx1.newcastle.edu.au:8087. The result was:lim Cn = 2e −2γ .

n→∞where γ denotes Euler’s constant. This is now proven.

Other Ising integral evaluations found using PSLQ
D2 = 1/3
D3 = 8 + 4π 2 /3 − 27 Li−3 (2)D4 = 4π 2 /9 − 1/6 − 7ζ(3)/2E2 = 6 − 8 log 2E3 = 10 − 2π 2 − 8 log 2 + 32 log2 2

E4 = 22 − 82ζ(3) − 24 log 2 + 176 log2 2 − 256(log3 2)/3+16π 2 log 2 − 22π 2 /3

E5 = 42 − 1984 Li4 (1/2) + 189π 4 /10 − 74ζ(3) − 1272ζ(3) log 2+40π 2 log2 2 − 62π 2 /3 + 40(π 2 log 2)/3 + 88 log4 2+464 log2 2 − 40 log 2

where ζ(x) is the Riemann zeta function and Lin (x) is the polylogarithm function.

The Ising integral E5
We were able to reduce E5, to an extremelycomplicated 3-D integral (see right).
We computed this integral to 250-digit
precision, using a highly parallel,
high-precision 3-D quadrature program.
Then we used a PSLQ program to discover
the evaluation given on the previous page.
We also computed D5 to 500 digits, butwere unable to identify it. In March 2014,
Erik Panzer proved our formula for E5 andfound a similar evaluation for D5 .
1. D. H. Bailey, J. M. Borwein and R. E. Crandall,“Integrals of the Ising class,” J. Physics A: Math.
and Gen., vol. 39 (2006), pg. 12271–12302.
2. E. Panzer, “Algorithms for the symbolic
integration of hyperlogarithms with applications
to Feynman integrals,” manuscript, 2014.

⇧ 1⇧ 1⇧ 1⇤
2(1 x)2 (1 y)2 (1 xy)2 (1 z)2 (1 yz)2 (1 xyz)2
=
0
0
0
⇤
⇥
4(x + 1)(xy + 1) log(2) y 5 z 3 x7 y 4 z 2 (4(y + 1)z + 3)x6 y 3 z y 2 + 1 z 2 + 4(y+⇥⇥1)z + 5) x5 + y 2 4y(y + 1)z 3 + 3 y 2 + 1 z 2 + 4(y + 1)z 1 x4 + y z z 2 + 4z
⇥
⇥ 3
⇥ 2
⇥ 2
2
2
2
+5) y + 4 z + 1 y + 5z + 4 x +3z
4z + 1 y4zy + 1 x(y(5z + 4)⇤
⌅ ⇤+4)x 1)] / (x 1)3 (xy 1)3 (xyz 1)3 + 3(y 1)2 y 4 (z 1)2 z 2 (yz⇥1)2 x6 + 2y 3 z 3(z 1)2 z 3 y 5 + z 2 5z 3 + 3z 2 + 3z + 5 y 4 + (z 1)2 z⇥⇥
5z 2 + 16z + 5 y 3 + 3z 5 + 3z 4 22z 3 22z 2 + 3z + 3 y 2 + 3 2z 4 + z 3 + 2⇥⇥
z 2 + z 2 y + 3z 3 + 5z 2 + 5z + 3 x5 + y 2 7(z 1)2 z 4 y 6 2z 3 z 3 + 15z 2⇥+15z + 1) y 5 + 2z 2 21z 4 + 6z 3 + 14z 2 + 6z 21 y 4 2z z 5 6z 4 27z 3⇥⇥
27z 2 6z + 1 y 3 + 7z 6 30z 5 + 28z 4 + 54z 3 + 28z 2 30z + 7 y 2 2 7z 5⇥⇥
+15z 4 6z 3 6z 2 + 15z + 7 y + 7z 4 2z 3 42z 2 2z + 7 x4 2y z 3 z 3⇥⇥
9z 2 9z + 1 y 6 + z 2 7z 4 14z 3 18z 2 14z + 7 y 5 + z 7z 5 + 14z 4 + 3⇥⇥
z 3 + 3z 2 + 14z + 7 y 4 + z 6 14z 5 + 3z 4 + 84z 3 + 3z 2 14z + 1 y 3 3 3z 5⇥⇥
+6z 4 z 3 z 2 + 6z + 3 y 2
9z 4 + 14z 3 14z 2 + 14z + 9 y + z 3 + 7z 2 + 7z⇥
+1) x3 + z 2 11z 4 + 6z 3 66z 2 + 6z + 11 y 6 + 2z 5z 5 + 13z 4 2z 3 2z 2⇥
+13z + 5) y 5 + 11z 6 + 26z 5 + 44z 4 66z 3 + 44z 2 + 26z + 11 y 4 + 6z 5 4⇥⇥
z 4 66z 3 66z 2 4z + 6 y 3 2 33z 4 + 2z 3 22z 2 + 2z + 33 y 2 + 6z 3 + 26⇥⇥⇥
z 2 + 26z + 6 y + 11z 2 + 10z + 11 x2 2 z 2 5z 3 + 3z 2 + 3z + 5 y 5 + z 22z 4⇥⇥+5z 3 22z 2 + 5z + 22 y 4 + 5z 5 + 5z 4 26z 3 26z 2 + 5z + 5 y 3 + 3z 4⇥⇥⇥
22z 3 26z 2 22z + 3 y 2 + 3z 3 + 5z 2 + 5z + 3 y + 5z 2 + 22z + 5 x + 15z 2 + 2z⇥
+2y(z 1)2 (z + 1) + 2y 3 (z 1)2 z(z + 1) + y 4 z 2 15z 2 + 2z + 15 + y 2 15z 4⇥⌅⇤
2z 3 90z 2 2z + 15 + 15 / (x 1)2 (y 1)2 (xy 1)2 (z 1)2 (yz 1)2⌅ ⇤⇥(xyz 1)24(x + 1)(y + 1)(yz + 1) z 2 y 4 + 4z(z + 1)y 3 + z 2 + 1 y 2⇥⇥⇥
4(z + 1)y + 4x y 2 1 y 2 z 2 1 + x2 z 2 y 4 4z(z + 1)y 3
z2 + 1 y2⇤⌅+4(z + 1)y + 1) 1) log(x + 1)] / (x 1)3 x(y 1)3 (yz 1)3[4(y + 1)(xy
⇥ 4
⇥ 3
⇥ 2
⇥
2
2
2
2
+1)(z + 1) x z
4z 1 y + 4x(x + 1) z
1 y
x +1 z4z 1⇥⇥⌅ ⇤y 2 4(x + 1) z 2 1 y + z 2 4z 1 log(xy + 1) / x(y 1)3 y(xy 1)3 (z⌅⇤1)34(z + 1)(yz + 1) x3 y 5 z 7 + x2 y 4 (4x(y + 1) + 5)z 6 xy 3 y 2 +⇥⇥⇥1) x2 4(y + 1)x 3 z 5 y 2 4y(y + 1)x3 + 5 y 2 + 1 x2 + 4(y + 1)x + 1 z 4 +⇥⇥y y 2 x3 4y(y + 1)x2 3 y 2 + 1 x 4(y + 1) z 3 + 5x2 y 2 + y 2 + 4x(y + 1)
⇥
⌅
⇤
⌅⇥⌅
y + 1) z 2 + ((3x + 4)y + 4)z 1 log(xyz + 1) / xy(z 1)3 z(yz 1)3 (xyz 1)3⇤⌅/ (x + 1)2 (y + 1)2 (xy + 1)2 (z + 1)2 (yz + 1)2 (xyz + 1)2 dx dy dz

E5

Box integrals
The following integrals appear in numerous applications:
Z 1
Z 1
s/2
Bn (s) :=···r12 + · · · + rn2
dR
0
0
Z 1
Z 1
s/2
···
∆n (s) :=(r1 − q1 )2 + · · · + (rn − qn )2
dRdQ
0

0

I Bn (1) is average distance of a random point from the origin.

I ∆n (1) is average distance between two random points.
I Bn (−n + 2) is average electrostatic potential in an n-cube whose origin has a unit

charge.
I ∆n (−n + 2) is average electrostatic energy between two points in a uniform

n-cube of charged “jellium.”
I Recently integrals of this type have arisen in neuroscience, e.g. the average

distance between synapses in a mouse brain.
I D. H. Bailey, J. M. Borwein and R. E. Crandall, “Box integrals,” Journal of Computational and Applied
Mathematics, vol. 206 (2007), pg. 196–208.

Sample evaluations of box integrals
n
any
1
2
2
2
2
2
2
3
3
3
3
3
3

s
even s ≥ 0
s 6= −1
-4
-3
-1
1
3
s 6= −2
-5
-4
-2
-1
1
3

Bn (s)
rational, e.g., : B2 (2) = 2/3
1
s+1
1
− π8
−4√

− 2√
2
log(1
+ 2)√
√
1
1
3√ 2 + 3 log(1 + √2)
732)
5 2 + 20 log(1 +
2
1
s 3
2F
,
−
;
;
−1
1
2+s
√2 21 2
1−√
6 3 − 12 π
− 32 2 arctan √12√√−3G + 32 π log(1 + 2) + 3 √Ti2(3 − 2 2)
13
2 + 3√ √− 4 π 1+ 2 log
1
1
3
−
π
+
log
2 + √3 
4√
24
2
2
1
7
3
5 3 − 60 π − 20 log 2 +Here F is hypergeometric function; G is Catalan; Ti is Lewin’s inverse-tan function.

Algebraic numbers in Poisson potential functions and lattice sums
Lattice sums arising from the Poisson equation have been studied widely inmathematical physics and also in image processing. We numerically discovered, andthen proved, that for rational (x, y ), the two-dimensional Poisson potential function
satisfies
φ2 (x, y ) =1 X cos(mπx) cos(nπy )
1
log α
=
π2
m2 + n2πm,n odd

where α is an algebraic number, i.e., the root of an integer polynomial
0 = a0 + a1 α + a2 α 2 + · · · + an α n
The minimal polynomials for these α were found by PSLQ calculations, with the(n + 1)-long vector (1, α, α2 , · · · , αn ) as input, where α = exp(8πφ2 (x, y )). PSLQ
returned the vector of integer coefficients (a0 , a1 , a2 , ·, an ) as output.
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), pg. 115201.

Samples of minimal polynomials found by PSLQ
k
5
6
7
8
9

10

Minimal polynomial for exp(8πφ2 (1/k, 1/k))
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

The minimal polynomial corresponding to x = y = 1/32 has degree 128, withindividual coefficients ranging from 1 to over 1056 . This PSLQ computation required
10,000-digit precision. See next page.
Other polynomials required up to 50,000-digit precision.

128-degree polynomial corresponding to x = y = 1/32
2
3
4
5
6
7
− 1 + 21888α + 5893184α + 15077928064α − 3696628330464α − 287791501240448α − 30287462976198976α + 4426867843186404992α
8
9
10
11
12
− 554156920878198587888α + 10731545733669133574528α + 120048731928709050250048α + 4376999211577765512726656α − 279045693458194222125366432α
13
14
15
16
+ 18747586287780118903854334848α − 643310226865188446831485766208α + 12047117225922787728443496655488α − 117230595100328033884939566091384α+ 667772184328316952814362214365568α17− 4130661734713288144037409932696512α− 1891420571205861612091802761809141088α

20+ 6279796382074485140847650604801614559872α

18+ 72313626239383964765274946226530432α+ 38770881730553471470590641060872686464α
23

21− 577943965397394779947709633563006963008α− 50438907678331243798448849245156136801232α− 1441007171934715336769224848138270812591296α

26

19

24

+ 5554617356232728647085822946642640269497472α

27− 20280024430170705107000630261773759070647328α

2930+ 99541720739995105011861264308551867164583808α − 754081464712315412970559119390477134883548736α
31
32+ 6271958646895434365874802435136411922022336128α − 45931349314815625339442690290912948480194150172α+ 280907040806572157908285324812126135484630889344α

33− 1427273782916972532576299009596755423149111059136α

34

35
36
+ 6055180299673737231932804443230077408291723908736α − 21609910939164553316101994301952988793013291135584α+ 65433275736596914909292838375737685959952141180288α

37+ 385709310577705218843549196766620216295554031550592α− 169928170513492897108417040254326115991438719391296α

39

38

40
− 801233230832691550861608914233661767474963249815792α

41
42+ 1706210557291030772074402183123327251333271061516160α − 4421210594351357102505784181831242174063263551938496α+ 14444199585866329915643888187597383540233619718619776α

43+ 169891313454945514927724813351516976839425267825908096α− 50968478530199956388487913417905125665738409426112032α

45

44− 506612996672385619931633440499093959534203673546181440α

46

47
48
+ 1330573388204326565144545192834096788469932897185696896α − 3069501638444045841407951432645059776135089489403138888α
49
50+ 6226636397646752257692349351542872634032398917736673152α − 11133383491631126059761752734485434504397040890449485504α+ 17601823309919260471943648355479182983209248554083752576α+ 31141043717679289808081270766611355726695735914995681664α2225+ 305806320133365055812520453224169520739712α

51− 24723027443995082126054012492323603544226813344022687712α

52

53− 35982430389670551550204799905599476866868765647852189248α

5428

128-degree polynomial corresponding to x = y = 1/32, continued
55
56+ 40292583920117898286863491450657424717015372825433076864α − 48512188214363976290470868896252008979896310883132967248α
57
58+ 69275112214095149977288310632868535966705567728055958400α − 114516830148561378617778209682642099604147034577152904128α+ 195760470467323759899736578743283333538805684128806803072α+ 468944248086031450001465269696090117959962662732817675648α+ 738516443137003178837650661261546833168555909499151978624α+ 738516443137003178837650661261546833168555909499151978624α+ 468944248086031450001465269696090117959962662732817675648α+ 195760470467323759899736578743283333538805684128806803072α5960− 317349593507106729834513764473487031789280056911012860320α

61

62− 622467103741378906100611838210632752408312516281305008960α

63

64− 781916756680856373187881889706233393197646662361906135622α

65

66− 622467103741378906100611838210632752408312516281305008960α

67

68− 317349593507106729834513764473487031789280056911012860320α

69

70− 114516830148561378617778209682642099604147034577152904128α

71
72+ 69275112214095149977288310632868535966705567728055958400α − 48512188214363976290470868896252008979896310883132967248α
73
74+ 40292583920117898286863491450657424717015372825433076864α − 35982430389670551550204799905599476866868765647852189248α
75
76+ 31141043717679289808081270766611355726695735914995681664α − 24723027443995082126054012492323603544226813344022687712α
7778+ 17601823309919260471943648355479182983209248554083752576α − 11133383491631126059761752734485434504397040890449485504α+ 6226636397646752257692349351542872634032398917736673152α+ 1330573388204326565144545192834096788469932897185696896α+ 169891313454945514927724813351516976839425267825908096α

79− 3069501638444045841407951432645059776135089489403138888α

81

82− 506612996672385619931633440499093959534203673546181440α

83− 50968478530199956388487913417905125665738409426112032α

84

85
86
+ 14444199585866329915643888187597383540233619718619776α − 4421210594351357102505784181831242174063263551938496α+ 1706210557291030772074402183123327251333271061516160α

87

88− 801233230832691550861608914233661767474963249815792α

89
90+ 385709310577705218843549196766620216295554031550592α − 169928170513492897108417040254326115991438719391296α
91
92+ 65433275736596914909292838375737685959952141180288α − 21609910939164553316101994301952988793013291135584α+ 6055180299673737231932804443230077408291723908736α

93

94− 1427273782916972532576299009596755423149111059136α

80

128-degree polynomial corresponding to x = y = 1/32, continued+ 280907040806572157908285324812126135484630889344α

95

96− 45931349314815625339442690290912948480194150172α

97
98+ 6271958646895434365874802435136411922022336128α − 754081464712315412970559119390477134883548736α+ 99541720739995105011861264308551867164583808α

99

100− 20280024430170705107000630261773759070647328α

101
102+ 5554617356232728647085822946642640269497472α− 1441007171934715336769224848138270812591296α+ 305806320133365055812520453224169520739712α

103

104
− 50438907678331243798448849245156136801232α

105
106+ 6279796382074485140847650604801614559872α− 577943965397394779947709633563006963008α
107
108
109
+ 38770881730553471470590641060872686464α− 1891420571205861612091802761809141088α+ 72313626239383964765274946226530432α− 4130661734713288144037409932696512α

110

111
112
+ 667772184328316952814362214365568α− 117230595100328033884939566091384α

113
114
115
116
+ 12047117225922787728443496655488α− 643310226865188446831485766208α+ 18747586287780118903854334848α− 279045693458194222125366432α+ 4376999211577765512726656α− 30287462976198976α

122

117

+ 120048731928709050250048α− 287791501240448α

123

118

+ 10731545733669133574528α119− 554156920878198587888α

120

124
125
126
127
128
− 3696628330464α+ 15077928064α+ 5893184α+ 21888α−α121+ 4426867843186404992α

Summary
I Numerical reproducibility is looming as a major challenge for large-scale scientificcomputations.
I One practical solution is to judiciously employ higher-precision arithmetic,

facilitated by tools to identify sensitivities and make needed code modifications.
I In most cases, only a few sensitive loops need to be computed using higher

precision (32 or 64 digits). Thus, the total run time is only slightly increased.
I Some interesting computations in experimental mathematics and mathematicalphysics require even more precision — hundreds or thousands of digits.
I Such computations are now feasible using modern high-performance computing

systems, together with relatively easy-to-use software facilities.

This talk is available at http://www.davidhbailey.com/dhbtalks/dhb-numerical-error.pdf.