High-Precision Computation and Reproducibility
David H Bailey, Lawrence Berkeley National Lab, USA
This talk is available at:
http://www.davidhbailey.com/dhbtalks/dhb-icerm.pdf

1Progress of scientific supercomputers:
Data from the Nov 2012 Top500 list

2Large-scale parallel computing andnumerical reproducibility· As numerical computations are ported from single-processor systems tolarge-scale, highly parallel supercomputers, problems are typically scaled
up by factors of millions.
· As a result, computations that previously had satisfactory numericalbehavior now may be highly ill-conditioned, and results are reproducible
to fewer digits (or none at all).
· Computational scientists who develop codes are, in most cases, notexperts in numerical analysis and may not realize the potential difficulties.
Example: Colleagues at LBNL who work with a large code for analyzing
Large Hadron Collider data recently reported that when they merely changed
the floating-point library (for transcendental functions, etc.), they observedcases where particle events no longer occurred.
§ Do their results have any numerical significance at all in such cases?

3

Large-scale parallel computing andnumerical reproducibility· Porting a code to a parallel computer inevitably destroys any specified
order of operations, particularly for global summations. As a result, digitfor-digit reproducibility cannot be guaranteed in most cases.
· Even after the port to the parallel system, the order of operation changes
when the number of processors used is changed, as is very often done byprogrammers dealing with batch queues.
One potential solution:
Perform global summations using double-double arithmetic, then reduceback to double precision for final results.

4Aren’t 64 bits enough?
Almost all scientific computers (from PCs to supercomputers) now featureIEEE-754 64-bit floating-point arithmetic. However, for a growing body ofnumerical algorithms and applications, 64 bits aren’t enough:· Ill-conditioned linear systems.
· Large summations, with cancellations.
· Long-time, iterative simulations.
· Large-scale simulations.
· Resolving small-scale phenomena.
· Studies in experimental mathematics – hundreds or thousands of digits.
Using high-precision arithmetic is often the easiest way to solve numericalproblems, even more sophisticated algorithmic solutions are possible.
The fact that high-precision arithmetic is not only useful, but is in factessential for some important applications, has encountered surprisingly
strong resistance from several quarters.
5

Innocuous-looking example where
standard precision is inadequate
Find a polynomial to fit the data (1, 1048579, 16777489, 84941299, 268501249,
655751251, 1360635409, 2523398179, 4311748609) for arguments 0, 1, …, 8.
The usual approach is to solve the linear system:

2
6
6
6
4

Pn n

k=1 xk..
Pn .

n
k=1 xk

Pn
xk
Pk=1
n
2
x
k=1 k
..
P . n+1
k=1 xk

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

Pn
32
3 2 Pn3
na0k
k=1 ykPnk=1 xPnn+1 7 6
7 6
7
x
a
x
y
1
k
k
k=1 k
k=1
76
7 6
7
7 6 .. 7 = 6
7..
..
54 . 5 4
5
.
.
P
Pn
n
n
2n
a
x
x
n
k=1 k ykk=1 kusing Matlab, Linpack or LAPACK. However, computation with 64-bit floating-point
arithmetic fails to find the correct result in this instance.

However, if the Linpack routines are converted to use double-double arithmetic (31digit accuracy), the above computation quickly produces the correct polynomial:
4
8
20
4
8

f (x) = 1 + 1048577x + x= 1 + (26+ 1)x + x

Algorithm changes versus double-double:
Double-double is the pragmatic choice· The result on the previous page can be obtained with double precisionusing Lagrange interpolation or the Demmel-Koev algorithm.
· But few computational scientists, outside of expert numerical analysts, areaware of these schemes – most people use Linpack or 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).
This is generated by:
6
12
27
6
12

f (x) = 1 + 134217729x + x= 1 + (2+ 1)x + x

In contrast, a straightforward Linpack scheme, implemented with doubledouble arithmetic, works fine for this and a wide range of similar problems.
Which is a more practical solution to a numerical anomaly?
(a) modify an existing code to use double-double, or(b) implement a new scheme from scratch?
With new easy-to-use double-double software, choice (a) is much preferable.
7Free software for high-precisioncomputation· ARPREC. Arbitrary precision, with many algebraic and transcendental functions.
High-level interfaces for C++ and Fortran-90 make code conversion easy.
Available at http://crd.lbl.gov/~dhbailey/mpdist.
· GMP. Produced by a volunteer effort and is distributed under the GNU license by
the Free Software Foundation. Available at http://gmplib.org.
· MPFR. C library for multiple-precision floating-point computations with exact
rounding, based on GMP. Available at http://www.mpfr.org.
· MPFR++. High-level C++ interface to MPFR. Available at
http://perso.ens-lyon.fr/nathalie.revol/software.html.
· GMPFRXX. Similar to MPFR++. Available at
http://math.berkeley.edu/~wilken/code/gmpfrxx.
· MPFUN90. Similar to ARPREC, but is written entirely in Fortran-90 and provides
only a Fortran-90 interface. Available at http://crd.lbl.gov/~dhbailey/mpdist.
· QD. This package perform “double-double” (approx. 31 digits) and “quaddouble” (approx. 62 digits) arithmetic. C++ and Fortran-90 high-level interfaces
makes code conversion easy. Available at http://crd.lbl.gov/~dhbailey/mpdist.
All of these packages greatly increase run time – from ~5X for double-double to~1000X for arbitrary precision with 1000 digits.
8

Berkeley’s CORVETTE project· CORVETTE: Correctness Verification and Testing of Parallel Programs· Tools to find bugs in hybrid (conventional/GPU) and large-scale systems.
· One key component: numerical reliability
§
§
§
§

Tools to easily test the level of numerical accuracy of an application.
Tools to delimit the portions of code that are the principal sources of
inaccuracy.
Tools to ameliorate numerical difficulties when they are uncovered, including
usage of double-double or higher precision arithmetic.
Tools to navigate through a hierarchy of precision levels (half, double,
double-double, higher).

9Some applications where
high-precision arithmetic is essential
·
·
·
·
·
·
·
·
·
·
·
·
·

Planetary orbit calculations (32 digits).
Supernova simulations (32-64 digits).
Climate modeling (32 digits).
Coulomb n-body atomic system simulations (32-120 digits).
Schrodinger solutions for lithium and helium atoms (32 digits).
Electromagnetic scattering theory (32-100 digits).
Scattering amplitudes of quarks, gluons and bosons (32 digits).
Discrete dynamical systems (32 digits).
Theory of nonlinear oscillators (64 digits).
Detecting “strange” nonchaotic attractors (32 digits).
The Taylor algorithm for ordinary differential equations (100-550 digits).
Ising integrals from mathematical physics (100-1000 digits).
Other examples from experimental mathematics (100-20,000 digits).

For details and references, see:David H. Bailey, Roberto Barrio, and Jonathan M. Borwein, "High precision computation: Mathematicalphysics and dynamics," Applied Mathematics and Computation, vol. 218 (2012), pg. 10106-10121.
http://www.davidhbailey.com/dhbpapers/hpmpd.pdf
10Long-term planetary orbits· Researchers have recognized for centuries thatplanetary 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 byexact 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 quaddouble arithmetic is required to avoid severe
inaccuracies, even if other techniques are employed
to reduce numerical error.
G. Lake, T. Quinn and D. C. Richardson, “From Sir Isaac to the Sloan survey: Calculating the structure andchaos due to gravity in the universe,” Proc. of the 8th ACM-SIAM Symp. on Discrete Algorithms, 1997, 1-10."11

Supernova simulations· Researchers at LBNL have used quaddouble arithmetic to solve for non-local
thermodynamic equilibrium populations ofiron and other atoms in the atmospheres ofsupernovas.
· Iron may exist in several species, so it isnecessary to solve for all species
simultaneously.
· Since the relative population of any statefrom the dominant state is proportional tothe exponential of the ionization energy,
the dynamic range of these values can bevery large.
· The quad-double portion now dominates
the entire computation.
P. H. Hauschildt and E. Baron, “The Numerical Solution of the Expanding Stellar Atmosphere Problem,”
Journal Computational and Applied Mathematics, vol. 109 (1999), pg. 41-63.
12

Climate modeling: high-precision
arithmetic for reproducibility· Climate and weather simulations arefundamentally chaotic – if microscopicchanges are made to the current state,soon the future state is quite different.
· In practice, computational results arealtered even if minor changes are
made to the code or the system.
· This numerical variation is a major
nuisance for code maintenance.
· He and Ding of LBNL found that byusing double-double arithmetic toimplement a key inner product loop,
most of this numerical variation
disappeared.
Y. He and C. Ding, “Using Accurate Arithmetics to Improve Numerical Reproducibility and Stability in Parallel
Applications,” Journal of Supercomputing, vol. 18, no. 3 (Mar 2001), pg. 259-277.
13

Coulomb n-body atomic systemsimulations· Alexei Frolov of Queen’s University in Canada has used high-precision
arithmetic to solve a generalized eigenvalue problem that arises in
Coulomb n-body interactions.
· Matrices are typically 5,000 x 5,000 and are very nearly singular.
· Computations typically involve massive cancellation, and high-precision
arithmetic must be employed to obtain numerically reproducible results.
· Frolov has also computed elements of the Hamiltonian matrix and theoverlap matrix in four- and five-body systems.
· These computations typically require 120-digit arithmetic.
“We can consider and solve the bound state few-body problems which havebeen beyond our imagination even four years ago.” – Frolov

A. M. Frolov and DHB, “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.
14

Error bounds in Chebyshev-Sobolev:
double vs high precision10

10

double precision

10

T4

0

10

0

10

T5

-10

10

-10

T5

-20

10

-30

Error

-30

10

-2

T4

10

Error

-20

10

multiple precision

10

10
-1

0

1

2

-2

-1

-0.5

0

1

2

point xpoint xBehavior of the theoretical error bounds (T4 a backward error bound andT5 for the running error bound) and the relative error in the double- andhigher-precision evaluation of the Chebyshev-Sobolev approximation ofdegree 50 of the function f(x) = (x+1)2 sin(4x), where the discrete Sobolev
measure have one mass point (c = 1) up to first derivative in the discretepart of the inner product.
R. Barrio and S. Serrano, “Generation and evaluation of orthogonal polynomials in discrete Sobolev spaces
II. Numerical stability, J. Comput. Appl. Math., vol. 181 (2005), pg. 299-320.
15Taylor’s method for the solution of ODEs
Taylor’s method is one of the oldest numerical schemes for solving ODEs, but inrecent years has re-emerged as the method of choice in the computational dynamics
community because of speed to convergence.
Consider the initial value problem y’ = f(t, y). The solution at time t = ti is:y(t0 )y(ti )

=:

y0 ,

'

1 dn 1 f (ti 1 , yi 1 ) n
yi 1 + f (ti 1 , yi 1 ) hi + · · · +
hi
n
1
n!
dt
yi

=:

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

Taylor’s method 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

20

0

y

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
double precision.
DHB, R. Barrio and J. M. Borwein, “High precision computation: Mathematical physics and dynamics,” Applied
Mathematics and Computation, vol. 218 (2012), pg. 10106-10121.
17

Computing the “skeleton” of periodic
orbitsJacobi 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

Symmetric periodic orbits (m denotes the multiplicity of the periodic orbit) in the mostchaotic zone of the (7+2) ring problem using double (A) and quadruple (B) precision.
R. Barrio and F. Blesa, “Systematic search of symmetric periodic orbits in 2DOF Hamiltonian systems,”
Chaos, Solitons and Fractals, vol. 41 (2009), 560-582.
18

Fractal properties of Lorenz attractors:
using high-precision to “zoom in”

On the first plot, the intersection of an arbitrary trajectory on the Lorenz attractor with thesection z = 27. The plot shows a rectangle in the x-y plane. All later plots zoom in on a tiny
region (too small to be seen by the unaided eye) at the center of the red rectangle of thepreceding plot to show that what appears to be a line is in fact not a line. Very high precision(hundreds of digits) arithmetic is required for these results. ""1. D. Viswanath, “The fractal property of the Lorenz attractor,” Journal of Physics D, vol. 190 (2004), 115-128."
2. D. Viswanath and S. Sahutoglu, “Complex singularities and the Lorenz attractor,” SIAM Review, to appear.
19

Lions-Mercer iterations
The Lions-Mercer iteration, also known as the Douglas-Rachford or Feinup iteration,
is defined by the procedure: reflect, reflect and average:

x 7! T (x) :=x + RA (RB (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: (α = 0.95 on left, and 1.0 on right):

20

Exploring iterations using Cinderella
Iterations such as this, as well as many other graphical phenomena, may be
explored using the Cinderella online tool: http://www.cinderella.de.
Two applets have been defined, working with Cinderella, for exploring LionsMercer iterations:A1. http://users.cs.dal.ca/∼jborwein/reflection.html
A2. http://users.cs.dal.ca/∼jborwein/expansion.html
For Applet A1, we observed that (see graphic on next slide):· As long as the iterate is outside the unit circle the next point is always
closer to the origin;· Once inside the circle the iterate never leaves;· The angle now oscillates to zero and the trajectory hence converges
to(1,0).

21

Iterations with Applet A1

22

Iterations with Applet A2:
Double vs multiple precision23

Experimental math: Discovering newmathematical results by computer· Compute various mathematical entities (limits, infinite series sums,
definite integrals) to high precision, typically 100-1000 digits.
· Use algorithms such as PSLQ to recognize these entities in terms of wellknown mathematical constants.
· When results are found experimentally, seek to find 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 justthe same as in physics.” – Kurt Godel

Mathematics

Computerscience

Scientificcomputing24

Mathematics

The PSLQ integer relation algorithm
Let (xn) be a given vector of real numbers. An integer relation algorithm
finds integers (an) such thata1 x1 + a2 x2 + · · · + an xn= 0(or within “epsilon” of zero, where epsilon = 10-p and p is the precision).
At the present time the “PSLQ” algorithm of mathematician-sculptor
Helaman Ferguson is the most widely used integer relation algorithm. It
was named one of ten “algorithms of the century” by Computing in Scienceand Engineering.
Integer relation detection requires very high precision (at least n*d digits,
where d is the size in digits of the largest ak), both in the input data and in
the operation of the algorithm.
1. H. R. P. Ferguson, DHB and S. Arno, “Analysis of PSLQ, An Integer Relation Finding Algorithm,”
Mathematics of Computation, vol. 68, no. 225 (Jan 1999), pg. 351-369.
2. DHB and D. J. Broadhurst, “Parallel Integer Relation Detection: Techniques and Applications,”
Mathematics of Computation, vol. 70, no. 236 (Oct 2000), pg. 1719-1736.
25

PSLQ, continued· PSLQ constructs a sequence of integer-valued matrices Bn that reduces
the vector y = x * Bn, until either the relation is found (as one of thecolumns of Bn), or else precision is exhausted.
· At the same time, PSLQ generates a steadily growing bound on the sizeof any possible relation.
· When a relation is found, the size of smallest entry of the y vector
suddenly drops to roughly “epsilon” (i.e. 10-p, where p is the number ofdigits of precision).
· The size of this drop can be viewed as a “confidence level” that therelation is real and not merely a numerical artifact -- a drop of 20+ ordersof magnitude almost always indicates a real relation.
Several efficient variants of PSLQ are available:· 2-level and 3-level PSLQ: performs almost all PSLQ iterations with only
double precision, updating full-precision arrays as needed. Hundreds oftimes faster than the original full-precision PSLQ algorithm.
· Multi-pair PSLQ: dramatically reduces the number of iterations required.
Designed for parallel system, but runs faster even on 1 CPU.
26

Decrease of log10(min |yi|) in multipair
PSLQ run27

PSLQ discovery:
The BBP formula for Pi
In 1996, this new formula for π was found using a PSLQ program:

=

⇤ 1
16k

k=0

4
8k + 1

2
8k + 41
8k + 51
8k + 6⇥This formula permits one to compute binary (or hexadecimal) digits of π
beginning at an arbitrary starting position, using a very simple scheme thatcan run on any system, using only standard 64-bit or 128-bit arithmetic.
Recently it was proven that no base-n formulas of this type exist for π,
except n = 2m.

1. DHB, 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.
28

High-precision tanh-sinh quadrature
Given f(x) defined on (-1,1), define g(t) = tanh (π/2 sinh t). Then setting x = g(t)yields⇥ 1

1

f (x) dx

=

⇥ ⇤

N

f (g(t))g ⇥ (t) dt⇤

wj f (xj ),hj= Nwhere xj = g(hj) and wj = g’(hj). Since g’(t) goes to zero very rapidly for large t, theproduct f(g(t)) g’(t) typically is a nice bell-shaped function for which the EulerMaclaurin formula implies that the simple summation above is remarkably accurate.
Reducing h by half typically doubles the number of correct digits.
For our applications, we have found that tanh-sinh is the best general-purpose
integration scheme for functions with vertical derivatives or singularities atendpoints, or for any function at very high precision (> 1000 digits). Otherwise weuse Gaussian quadrature.
1. DHB, 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.
29

Ising integrals from mathematicalphysics
We recently applied our methods to study three classes of integrals thatarise in the Ising theory of mathematical physics – Dn and two others:

Cn

Dn

En

:=

:=

=

4
n!

⌦ ⇤

4
n!

2

···

⌦ ⇤

⌦ ⇤

···

⌦ ⇤

⌦ 1

⌦ 1

0

0

0

···

0

0

0

⇤
⇧

du1
dun
·
·
·
⇥2
⌥n
u1
un
(u
+
1/u
)
j
j=1 j
⇥2
ui uj
i<j ui +uj
du1
dun
···
⇥2
⌥n
u1
un
(u
+
1/u
)
j
j=1 j
⌅2
uk uj ⌃
dt2 dt3 · · · dtn
uk + uj
1

1⇥j<k⇥nwhere in the last line uk = t1 t2 … tk.

DHB, J. M. Borwein and R. E. Crandall, “Integrals of the Ising class,” Journal of Physics A: Mathematicaland General, vol. 39 (2006), pg. 12271-12302.
30

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):
n Z 1
2
Cn =
tK0n (t) dt
n! 0

1000-digit numerical values, computed using this formula, approach a limit:

C1024 = 0.63047350337438679612204019271087890435458707871273234 . . .
What is this limit? We copied the first 50 digits of this numerical value into
the online Inverse Symbolic Calculator (ISC):
http://ddrive.cs.dal.ca/~isc or http://carma-lx1.newcastle.edu.au:8087/

The result was:lim Cn

n⇥⇤

=

2e 2

where γ denotes Euler’s constant.
31

Other Ising integral evaluations foundusing high-precision PSLQ
D2

=

1/3

D3

=

8 + 4⇥ 2 /3

D4=

4⇥ 2 /9E2=

68 log 2

E3

=

10

2⇥ 2

E4

=

22

82 (3)

27 L 3 (2)1/6

8 log 2 + 32 log2 2+16⇥ 2 log 2
E5

?

=

42

7 (3)/2

24 log 2 + 176 log2 2256(log3 2)/3

22⇥ 2 /3

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 mathematical physicist Craig Tracy, who hopedthat our tools could help identify D5.
32

The Ising integral E5
We were able to reduce E5, which
is a 5-D integral, to an extremelycomplicated 3-D integral.
We computed this integral to 250digit precision, using a highly
parallel, high-precision 3-D
quadrature program. Then weused a PSLQ program to discover
the evaluation given on theprevious page.
We also computed D5 to 500
digits, but were unable to identify
it. The digits are available ifanyone wishes to further explore
this question.

⇧ 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⇥⇥⇥⇥
+5) y 2 + 4 z 2 + 1 y + 5z + 4 x3 +3z 2 4z + 1 y 2 4zy + 1 x2 (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
⇥
3
2
4
3
2
+1) x + z 11z + 6z
66z + 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⇥⇥⇥⇥+1)(z + 1) x2 z 2 4z 1 y 4 + 4x(x + 1) z 2 1 y 3
x2 + 1 z 2 4z 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 +
⇥
⇥
⇥
2
1) x
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

E533Recursions in Ising integrals
Consider the 2-parameter class of Ising integrals (which arises in QFT for odd k):

Cn,k

=

4
n!

⌅

0

···

⌅

0

du1
dun
···
⇥k+1
⇤n
u1
un
(u
+
1/u
)
j
j=1 j
1

After computing 1000-digit numerical values for all n up to 36 and all k up to 75(performed on a highly parallel computer system), we discovered (using PSLQ)
linear relations in the rows of this array. For example, when n = 3:

0 = C3,084C3,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, but more complicated, recursions have been found for all n.
1. DHB, D. Borwein, J. M. Borwein and R. Crandall, “Hypergeometric Forms for Ising-Class Integrals,”
Experimental Mathematics, vol. 16 (2007), pg. 257-276.
2. J. M. Borwein and B. Salvy, “A Proof of a Recursion for Bessel Moments,” Experimental Mathematics, vol. 17(2008), pg. 223-230.
34

Box integrals
The following integrals appear in numerous arenas of math and physics:
⇤ 1
⇤ 1
⇥
2
2 s/2
Bn (s) :=···r1 + · · · + r n
dr1 · · · drn
0

n (s)

:=

⇤ 1
0

0

···

⇤ 1
0

(r1

q1 ) + · · · + (rn
2

⇥
2 s/2

qn )

dr1 · · · drn dq1 · · · dqn• Bn(1) is the expected distance of a random point from the origin of n-cube.
• Δn(1) is the expected distance between two random points in n-cube.
• Bn(-n+2) is the expected electrostatic potential in an n-cube whose origin
has a unit charge.
• Δn(-n+2) is the expected electrostatic energy between two points in a
uniform n-cube of charged “jellium.”• Recently integrals of this type have arisen in neuroscience – e.g., theaverage distance between synapses in a mouse brain. "DHB, J. M. Borwein and R. E. Crandall, “Box integrals,” Journal of Computational and Applied Mathematics,
vol. 206 (2007), pg. 196-208.

35

Evaluations of box integrals
n
any
1
2
2
2
2
2
2
3
3
3
3
3
3

s
even s ⇥ 0
s ⇤= 1-4-3-1
13s ⇤= 2
-5
-4
-2
-1
1
3

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

2⌅
2 log(1 + 2)⌅
⌅
1
1
2
+
3⌅
3 log(1 + ⌅2)
7
3
2
+
2)⇥
5
20 log(1 +
2
1
s 3
F
,
2
1
2+s
2
2; 2; 1
⌅
1
1
3
6
12
⌅
3
1
2
arctan2
2⌅⌅
3
3G + 2 log(1 + 2) + 3 ⌅Ti2⇥(3 2 2)
1
3
+
3⌅ ⇥
4
2 log 2 +
⌅
1
1
1
3
+
4⌅
24
2 log 2 + ⌅3 ⇥
1
7
2
3
3
5
60
20 log 2 +Here F is hypergeometric function; G is Catalan; Ti is Lewin’s inverse-tan function.
36

Elliptic function integrals
The research with ramble integrals led us to study integrals of the form:

I(n0 , n1 , n2 , n3 , n4 ) :=Z 1

xn0 K n1 (x)K 0n2 (x)E n3 (x)E 0n4 (x)dx,0

where K, K’, E, E’ are elliptic integral functions:

K(x) :=

Z 1
0

p

dt
(1

t2 )(1p
K (x) := K( 1 x2 )Z 1p
1 x2 t2p
E(x) :=
dt
2
1 t
0
p
E 0 (x) := E( 1 x2 )0

x2 t2 )

J. Wan, “Moments of products of elliptic integrals,” Advances in Applied Mathematics, vol. 48 (2012),
available at http://carma.newcastle.edu.au/jamesw/mkint.pdf."
DHB 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.
37

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. For example:
Z 1

x3 K 2 (x)E(x)dx =+51

x3 K 3 (x)dx + 32

81

0

243

Z 1

?

Z 1
0

6

Z 1

x3 K(x)E(x)K 0 (x)dx =+156

x3 K 3 (x)dx

624

0

20736

Z 1

59
Z 1

Z 1

1284

x3 K 3 (x)dx + 5136

0

972

Z 1

K 3 (x)dx + 468

0x4 K 3 (x)dx

135

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

Z

Z 1

Z 1

x2 K 3 (x)dx0xK(x)E(x)K 0 (x)dx0x4 E 2 (x)K 0 (x)dx = 3901
Z 1

x2 K 3 (x)dx

0

0

?

0

24

Z 1

x4 K 3 (x)dx0

Z 1

K 3 (x)dx

0

?

0

Z 1

Z 1

K 3 (x)dx

3852

0

0

Z 1
0

38

x2 K 3 (x)dx0x4 K 3 (x)dx
8316

Z 1

2592

Z 1

x2 K 2 (x)K 0 (x)dx0xK(x)E(x)K 0 (x)dx.

Algebraic numbers in Poisson potentialfunctions associated with lattice sums
Lattice sums arising from the Poisson equation have been studied widely inmathematical physics and also in image processing.
In two 2012 papers (below), we numerically discovered, and then proved, thatfor rational (x, y), the two-dimensional Poisson potential function satisfies

1
2 (x, y) = 2⇡

Xm,n odd

1
cos(m⇡x) cos(n⇡y)= log ↵
22
m +n⇡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, withthe (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. DHB, J. M. Borwein, R. E. Crandall and J. Zucker, “Lattice sums arising from the Poisson equation,”
manuscript, http://www.davidhbailey.com/dhbpapers/PoissonLattice.pdf"
2. DHB and J. M. Borwein, “Compressed lattice sums arising from the Poisson equation: Dedicated toProfessor Hari Sirvastava,” manuscript, http://www.davidhbailey.com/dhbpapers/Poissond.pdf.
39

Samples of minimal polynomials foundby 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↵735231↵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↵622321584↵7 + 39775986↵8 44431044↵9 + 19899882↵10 + 3546576↵118458020↵12 + 4009176↵13 273348↵14 + 121392↵1511385↵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, withindividual coefficients ranging from 1 to over 1056. This PSLQ computation
required 10,000-digit precision. See next slide.
40

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

41

Formal proof versus high-precisionverification· Results such as those mentioned above must still be proven rigorously.
· However, strong numerical evidence is often a good impetus to find a
proof – “discovery is 9/10 of the proof.”
What is more firmly established?
· A formally proven result, whose proof required hundreds of pages, which
crucially relies on tens of earlier results by other authors, and which has
only been read in detail by a handful of mathematicians.
· A numerically discovered experimental identity, for which no known formal
proof is available, but which has been checked to thousands of digits,
independently on separate computers.

42

Cautionary Example
These constants agree to 42 decimal digit accuracy, but are NOT equal:⇥cos(2x)0

cos(x/n) dx =n=1

0.392699081698724154807830422909937860524645434187231595926=

80.392699081698724154807830422909937860524646174921888227621
Richard Crandall has now shown that this integral is merely the first term ofa very rapidly convergent series that converges to π/8:

8

=
m=0

⇤

cos[2(2m + 1)x]0⇥

cos(x/n) dx

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

Limitations of Mathematica and Maple
for computational mathematics
Mathematica or Maple is our first choice whenever symbolic or numericcomputations are required. However, both have their limitations.
For example, in a study of Mordell-Tornheim-Witten sums (which arise inmathematical physics), we required high-precision numeric values ofderivatives with respect to the order s of polylogarithm functions:@Lis (z)
,
@s

where

Lis (z) =

1
X
zk

k=1

ks

Maple is not able to numerically evaluate these derivatives at all.
Mathematica, when asked for 4000 digits, returned only 400 correct digits (atsome arguments).
DHB, J. M. Borwein and R. E. Crandall, “Computation and theory of extended Mordell-Tornheim-Witten sums,”
Mathematics of Computation, to appear, 31 Jul 2012, http://www.davidhbailey.com/dhbpapers/BBC.pdf
44Summary· Large-scale, highly parallel computation places enormous stress on the
numerical reliability and reproducibility of scientific computations.
· Double-double or higher precision arithmetic is a practical means ofdealing with these numerical difficulties in many cases.
· Many real-world applications have now been identified that require highprecision arithmetic.
· Some research studies, particularly in experimental mathematics andmathematical physics, require hundreds or even thousands of digits.
· Software is now available, mostly for free, to facilitate conversion. In
most cases, one need only change the type of variables that are to betreated as high-precision variables.

45