High-Precision Computation:
Applications and Challenges
David H. Bailey
http://www.davidhbailey.com
Computational Research Dept., Lawrence Berkeley Natl. Lab.
Computer Science Dept., University of California, Davis

7 April 2013

How much precision is needed?
For many applications, 32-bit IEEE floating-point arithmetic is adequate:
I Most graphics and gaming applications.
I Most applications involving experimental data.
I Most signal processing applications.

However, some definitely require 64-bit IEEE arithmetic:
I Certain graphical calculations (e.g., to determine whether two polygons or curved

figures intersect).
I Certain structural problems (elasticity calculations, etc.).

For additional discussion and details:
1. W. Kahan, “A floating-point trick to solve boundary-value problems faster,” May 2007,
http://www.eecs.berkeley.edu/~wkahan/Math128/FloTrik.pdf.
2. W. Kahan, “On the cost of floating-point computation without extra-precise arithmetic,” Nov
2004, http://www.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf.

Aren’t 64 bits enough?
For some other algorithms and applications, even 64-bit IEEE arithmetic isn’t sufficient:
I Certain applications with highly ill-conditioned linear systems.
I Large summations, especially those involving cancellations.
I Long-time, iterative simulations.
I Large-scale, highly parallel simulations — numerical sensitivities are magnified.
I Resolving small-scale phenomena.
I Studies in experimental mathematics.

The notion that high-precision arithmetic is not only useful, but is in fact essential fornumerous important applications has met significant resistance in the computing
community.
D.H. Bailey, R. Barrio, and J.M. Borwein, “High precision computation: Mathematical physics and dynamics,” Applied Mathematics andComputation, vol. 218 (2012), pg. 10106–10121.

Innocuous example where standard 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

Algorithm changes versus double-double?
Double-double is the pragmatic choice
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—most people use home-grown code.
Besides, even these schemes fail for higher-degree problems. For example:(1, 134217731, 8589938753, 97845255883, 549772595201,
2097396156251, 6264239146561, 15804422886323, 35253091827713,
71611233653971, 135217729000001, 240913322581691, 409688091758593) isgenerated by:
f (x) = 1 + 134217729x 6 + x 12 = 1 + (227 + 1)x 6 + x 12
In contrast, a straightforward Linpack scheme, implemented with double-double
arithmetic, works fine for this and a wide range of similar problems.

Numerical analysis expertise among U.C. Berkeley grads
Of the 2010 U.C. Berkeley graduating class, 870 were in disciplines expected to use
technical computing:
I Division of Mathematical and Physical Sciences (Mathematics, Physics,

Statistics).
I College of Chemistry.
I College of Engineering (including Computer Science).
I Other fields with significant computing: biology, geology, medicine and social

sciences.
Enrollment in numerical analysis:
I Math 128A (Introductory numerical analysis required of math majors): 219.
I Math 128B (A more advanced course): 24.

Conclusion: Fewer than 2% of Berkeley grads who will do technical computing have
had rigorous training in numerical analysis!

Reproducibility in scientific computing
The report of the recent NSF-funded Workshop on Reproducibility in Computational
and Experimental Mathematics noted:
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.
The report concluded that high-precision arithmetic was one of the best ways toameliorate these numerical difficulties.
V. Stodden, D.H. Bailey, J. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default to reproducible: Reproducibility in computationaland experimental mathematics,” http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

Enhancing reproducibility with 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, it takes

2.5

2.0

1.5

1.0

0.5

47.39 seconds seconds and yields the result
7.073157029007832.
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, π).

Techniques for high-precision computation
I Data is stored as a string of ints or floats, with initial words holding the string

length and binary exponent.
I For modest precision levels (< 1000 digits):
1. Use conventional grade-school schemes for basic arithmetic.
2. Use conventional Taylor series schemes for transcendentals.
I For extra-high precision levels (> 1000 digits):

1. Use FFT-based multiplication and Taylor series division.
2. Use quadratically convergent algorithms for transcendentals.
I Operator overloading (available in Fortran-90, C++ and other languages) can be

exploited to construct high-level translation interfaces, which greatly simplify
conversion of existing code.
1. D.H. Bailey, Y. Hida, X.S. Li and B. Thompson, “ARPREC: An arbitrary precision computation package,”
http://www.davidhbailey.com/dhbpapers/arprec.pdf.
2. Y. Hida, X.S. Li and D.H. Bailey, “Algorithms for quad-double precision floating point arithmetic,” 15th IEEE Symposium on Computer
Arithmetic, IEEE Computer Society, 2001, pg. 155–162.

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
Developing software tools to find and ameliorate numerical anomalies in large-scalecomputations:
I Tools to test the level of numerical accuracy required for an application.
I Tools to delimit the portions of code that are inaccurate.
I Tools to repair numerical difficulties, including usage of high-precision arithmetic.
I Tools to navigate through a hierarchy of precision levels (32-bit, 64-bit or higher

as needed).
ARITH21 talks on CORVETTE research:
1. Hong Diep Nguyen and Jim Demmel (UC Berkeley, USA), “Numerical Accuracy andReproducibility at ExaScale.”
2. Hong Diep Nguyen and James Demmel (UC Berkeley, USA), “Fast Reproducible Floating-Point
Summation.”

Some applications where high precision is 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. Detecting “strange” nonchaotic attractors (32 digits).
11. The Taylor algorithm for ODEs (100–600 digits).
12. Ising integrals from mathematical physics (100–1000 digits).
13. Problems in experimental mathematics (100–50,000 digits).

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 required to
avoid severe inaccuracies, even if other techniques are employed to reduce numericalerror.
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.
I The quad-double portion now dominates the entire

computation.

P.H. Hauschildt and E. Baron, “Thenumerical solution of the expanding stellar
atmosphere problem,” Journal
Computational and 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 numericalreproducibility and stability 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.”
A.M. Frolov and D.H. Bailey, “Highly accurate evaluation of the few-body auxiliary functions and four-body integrals,” Journal of Physics B, vol. 36,
no. 9 (14 May 2003), pg. 1857–1867.

Taylor’s method for the solution of ODEs
Taylor’s method is one of the oldest numerical schemes for solving ODEs, but in recent
years has re-emerged as a leading choice for problems in computational dynamics
community. Consider the initial value problem y 0 = f (t, y ). The solution at time t = tiis:y (t0 ) =: y0 ,y (ti )'

yi−1 + f (ti−1 , yi−1 ) hi + · · · +=: yi

1 d n−1 f (ti−1 , yi−1 ) nhin!
dt n−1The Taylor coefficients here may be found using automatic differentiation methods.
One significant advantage with Taylor’s method is that it can be easily implemented
using high-precision arithmetic. When this is done, Taylor’s method typically gives
superior results, compared with other available schemes.
A. Abad, R. Barrio, F. Blesa and M. Rodriguez, “TIDES: a Taylor series Integrator for Differential EquationS,” preprint, 2010.

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 periodic orbit for the Lorenz model
during 16 time periods using the TIDES code with 300 digits, and 1 time periods using

Computing the “skeleton” of periodic orbits
Jacobi constant C

−4

limit
m=1
m=2
m=3
m=4

−5

−6

A
−7

−8
−5

−4

−3

−2

−1
coordinate x

0

1

2

−4

Jacobi constant C

B
−5

−6

−7

−8
−5

−4

−3

−2

−1
coordinate x

0

1

2

Fractal properties of Lorenz attractors

On the first plot, the intersection of an arbitrary trajectory on the Lorenz attractor
with the section z = 27. The plot shows a rectangle in the x-y plane. All later plotszoom in on a tiny region (too small to be seen by the unaided eye) at the center of theLions-Mercer iterations
The Lions-Mercer iteration, also known as the Douglas-Rachford or Feinup iteration, isdefined by “reflect, reflect and average”:
x + RA (RB (x))
x 7→ T (x) :=2
In the simple 2-D case of a horizontal line of height α, we obtain the explicit iteration:
xn+1 := cos θn ,yn+1 := yn + α − sin θn ,(θn := arg zn )For 0 < α < 1, spiraling is ubiquitous. Compare α = 0.95 on left with α = 1.0 onright:

Exploring iterations using Cinderella
Iterations such as this, as well as many other graphical phenomena, may be exploredusing the Cinderella online tool: http://www.cinderella.de.
Two applets have been defined, working with Cinderella, for exploring Lions-Mercer
iterations:A1. http://users.cs.dal.ca/?jborwein/reflection.html
A2. http://users.cs.dal.ca/?jborwein/expansion.html
With Applet A1, we observe that as long as the iterate is outside the unit circle, thenext point is always closer to the origin; once inside the circle the iterates never leaveand converge to (1, 0).
Applet A2 demonstrates clear numerical degradation, compared with similar run using
high-precision arithmetic.

Iterations with Applet A1

Iterations with Applet A2: Double vs multiple precision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.
1. H.R.P. Ferguson, D.H. Bailey and S. Arno, “Analysis of PSLQ, An Integer Relation Finding Algorithm,” Mathematics of Computation, vol.
68, no. 225 (Jan 1999), pg. 351–369.
2. 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

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)

=

1X 1
8
64k



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
−
+
+
1792
212k (24k + 1)3(24k + 2)3(24k + 3)3(24k + 4)3



k=0

1536
3072
768
3072
2688
−
+
+
−
−
(24k + 5)3(24k + 6)3(24k + 7)3(24k + 9)3(24k + 10)3
192
1536
96
672
384
−
−
−
−
−
(24k + 11)3
(24k + 12)3
(24k + 13)3
(24k + 14)3
(24k + 15)3
24
48
12
120
48
+
+
−
+
+
(24k + 17)3
(24k + 18)3
(24k + 19)3
(24k + 20)3
(24k + 21)3

42
3
−
+
(24k + 22)3
(24k + 23)3

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.
1. 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.
2. H. Takahasi and M. Mori, “Double Exponential Formulas for Numerical Integration,” Publications of RIMS, Kyoto University, vol. 9 (1974),
pg. 721-741.

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
du1
dun
4 ∞
···
···
Cn :=
P
2
n! 0
u1
un
n
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 .

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 L−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 ζ is the Riemann zeta function and Lin (x) is the polylog function. D2 , D3 andD4 were originally provided to us by Craig Tracy, who hoped that our tools could help
identify D5 .

The Ising integral E5
We were able to reduce E5, which is a 5-D
integral, to an extremely complicated 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. The digits areavailable if anyone wishes to further
explore.

⇧ 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

Recursions in Ising integrals
Consider the 2-parameter class of Ising integrals (which arises in quantum field theory
for odd k):
Z ∞
Z ∞
du1
dun
4
1
Cn,k =
···
···
P
k+1
n! 0
u1
un
n
0
j=1 (uj + 1/uj )After computing 1000-digit numerical values for all n up to 36 and all k up to 75, wediscovered (using PSLQ) linear relations in the rows of this array. For example, whenn = 3:
0 = C3,0 − 84C3,2 + 216C3,4

0 = 2C3,1 − 69C3,3 + 135C3,5
0 = C3,2 − 24C3,4 + 40C3,6

0 = 32C3,3 − 630C3,5 + 945C3,7

0 = 125C3,4 − 2172C3,6 + 3024C3,8

Similar recursions have been found (and proven) for all n.

1. D.H. Bailey, D. Borwein, J.M. Borwein and R.E. Crandall, “Hypergeometric Forms for Ising-Class Integrals,” Experimental Mathematics, vol.

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

Elliptic integral functions
Research with “ramble” integrals led us to consider these integrals:
Z 1
I (n0 , n1 , n2 , n3 , n4 ) :=x n0 K n1 (x)K 0n2 (x)E n3 (x)E 0n4 (x)dx,0

where K , K 0 , E , E 0 are elliptic integral functions:
Z 1

dt
p
2
(1 − t )(1 − x 2 t 2 )0p
K 0 (x) := K ( 1 − x 2 )Z 1√
1 − x 2t 2√
E (x) :=
dt
1 − t20p
E 0 (x) := E ( 1 − x 2 )K (x) :=

1. J. Wan, “Moments of products of elliptic integrals,” Advances in Applied Mathematics, vol. 48 (2012), pg. 121–141,
http://carma.newcastle.edu.au/jamesw/mkint.pdf.
2. D.H. Bailey and J.M. Borwein, “Hand-to-hand combat with thousand-digit integrals,” Journal of Computational Science, vol. 3 (2012), pg.
77–86, http://www.davidhbailey.com/dhbpapers/combat.pdf.

Relations found among the I integrals
Thousands of relations have been found among the I integrals. For example, among
the class with n0 ≤ D1 = 4 and n1 + n2 + n3 +n4 = D2 = 3 (a set of 100 integrals),
we found that all can be expressed in terms of an integer linear combination of 8
simple integrals. Some examples:
Z 1
81

x 3 K 2 (x)E (x)dx

Z 1

?

−6

=

0

Z 1

x 3 K 3 (x)dx + 32

Z 1

0

Z 1

Z 1

0

+51
−243

K 3 (x)dx − 24

x 4 K 3 (x)dx0

x 3 K (x)E (x)K 0 (x)dx

?

=

Z 1
−59

0

K 3 (x)dx + 468

Z 1

0

Z 1
+156

x 3 K 3 (x)dx − 6240Z 1−20736Z 1?

=

x 4 K 3 (x)dx − 135

Z 1

xK (x)E (x)K 0 (x)dx

0

Z 1
3901

K 3 (x)dx − 3852

Z 1

0

0

Z 1
−1284
0
Z 1

x 3 K 3 (x)dx + 5136

x 2 K 3 (x)dx00

x 4 E 2 (x)K 0 (x)dx−972

x 2 K 3 (x)dx0Zx 2 K 3 (x)dx0

x 4 K 3 (x)dx − 2592

Z 1

x 2 K 2 (x)K 0 (x)dx0

K (x)E (x)K 0 (x)dx − 8316

Z 1

xK (x)E (x)K 0 (x)dx.

Algebraic numbers in Poisson potential functions associated with latticesums
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
1 X cos(mπx) cos(nπy )1φ2 (x, y ) =
=
log α
π2
m2 + n2πm,n odd

where α is algebraic, 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(πφ2 (x, y )). PSLQ
returns the vector of integer coefficients (a0 , a1 , a2 , · · · , an ) as output.

1. 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, http://www.davidhbailey.com/dhbpapers/PoissonLattice.pdf.
2. D.H. Bailey and J.M. Borwein, “Compressed lattice sums arising from the Poisson equation: Dedicated to Professor Hari Sirvastava,”
Boundary Value Problems, vol. 75 (2013), DOI: 10.1186/1687-2770-2013-75,
http://www.boundaryvalueproblems.com/content/2013/1/75.

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 for exp(8πφ2 (1/32, 1/32)) has degree 128, with individual
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.

Degree-128 minimal polynomial for exp(8πφ2 (1/32, 1/32))

Cautionary example
These constants agree to 42 decimal digits, but are NOT equal:
Z ∞0cos(2x)∞Ycos(x/n) dx=

n=1

0.392699081698724154807830422909937860524645434187231595926
π=
80.392699081698724154807830422909937860524646174921888227621

Richard Crandall has shown that this integral is merely the first term of a very rapidly
convergent series that converges to π/8:
∞
∞ Z ∞
X
Y
π
=
cos[2(2m + 1)x]cos(x/n) dx
8
0
m=0

n=1

1. D.H. Bailey, J.M. Borwein, V. Kapoor and E. Weisstein, “Ten Problems in Experimental Mathematics,” American Mathematical Monthly,
vol. 113, no. 6 (Jun 2006), pg. 481–409.
2. R.E. Crandall, “Theory of ROOF Walks,” 2007, available at http://people.reed.edu/~crandall/papers/ROOF.pdf.

Limitations of Maple and Mathematica
Maple and Mathematica are our first choices whenever symbolic or numericcomputations are required. However, both have limitations and bugs.
For example, in a study of Mordell-Tornheim-Witten sums (which arise inmathematical physics), we required high-precision numeric values of derivatives withrespect to the order s of polylogarithms:∂Lis (z)
,
∂s

where

Lis (z) =

∞
X
zk
k=1

ks

Maple was not able to numerically evaluate these derivatives at all. Mathematica,
when asked for 4000 digits, returned only 400 correct digits (at some arguments).
D.H. Bailey, J.M. Borwein and R.E. Crandall, “Computation and theory of extended Mordell-Tornheim-Witten sums,” Mathematics of Computation,
Ramanujan Journal, 27 Feb 2013, DOI 10.1007/s11139-012-9427-1, http://www.davidhbailey.com/dhbpapers/BBC.pdf.

What is needed for high-precision floating-point software?
I A high-performance, rock-solid-reliable arithmetic engine, with precision scalable

to 1,000,000 digits or more.
I A separate package for modest precision (32 and 64 digits)?
I FFT-based multiplication for > 1000 digits.
I A thread-safe design to facilitate multicore parallel processing, and a pathway toextend to graphics processing units (GPUs).
I A comprehensive library of tuned transcendentals: not just sin, cos, exp, etc., butalso the gamma function, polylogs (with real and complex arguments), Bessel
functions, etc.
I A robust high-level language interface for C++, Fortran-90 and possibly several

other languages as well.
I Interfaces for Maple and Mathematica.

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