Computational challenges in experimental mathematics
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Laboratory (retired)Computer Science Department, University of California, Davis

21 Jul 2014

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

I V. Stodden, D. H. Bailey, J. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default toreproducible: Reproducibility in computational and experimental mathematics,” Jan 2013, available at
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

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?

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,” Proceedings of SC13, 2013.

Innocuous example where high precision is required for reproducible results
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, continued
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
Neither the Lagrange or Demmel-Koev algorithms get the right answer with double
precision. In contrast, a straightforward Linpack scheme, implemented withdouble-double arithmetic, works fine for this and a wide range of similar problems.

Numerical analysis experts or high precision?
2010 U.C. Berkeley statistics:
I 870 seniors graduated in the Division of Mathematical and Physical Sciences(including Mathematics, Physics and Statistics), the College of Chemistry and theCollege of Engineering (including Computer Science).
I Many graduates in other fields (biology, economics, finance, medicine, psychology,etc.) will also do significant numerical computing in their professional work.
I 219 students enrolled in two sections of Math 128A, a one-semester introductory

numerical analysis course required of applied math majors.
I Only 24 enrolled in Math 128B, a more advanced course.

Conclusion: Only about 2% of Berkeley graduates that will do scientific computing intheir career work have had advanced training in numerical analysis.

Other applications where high-precision arithmetic is useful or essential
1. Planetary orbit calculations (32 digits).
2. Supernova simulations (32–64 digits).
3. Certain components of 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.

High-precision arithmetic in experimental mathematics
Typical methodology:
1. Compute various mathematical entities (limits, infinite series sums, definite
integrals, etc.) to high precision, typically 100–10,000 digits.
2. Use an integer relation algorithm such as “PSLQ” to recognize these numericalvalues in terms of well-known mathematical constants.
3. Devise formal mathematical proofs of the discovered relations.
This approach has been extremely fruitful — literally thousands of new results.
What’s more, this is a very accessible form of research — even undergraduates(if they have good programming skills) can make state-of-the-art contributions.

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 for this constant, assuming that it is of degree
no more 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. 897–925.

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 studied these three families of integrals: Cn arise in quantum field theory (see next
viewgraph); Dn arise in Ising theory of mathematical physics (referred to us by Craig
Tracy of U.C. Davis); En are the numerators of Dn :
Z
Z ∞
dun
4 ∞
1
du1
···
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 .

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)
73
2)
5 2 + 20 log(1 +

2
1
s 3
F
,
−
;
;
−1
2
1
2+s
2
2 2
√
1
1
−√
6 3 − 12 π3
− 2 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.

Lessons learned from the Ising integral, box integral and similar projects
I High-precision computation (typically 400–4000 digits), combined with the PSLQ

algorithm, is very effective in evaluating and identifying definite integrals.
I The tanh-sinh scheme is very effective in evaluating a broad range of integrals tovery high precision, even including functions with singularities at endpoints.
I For most of these problems, the numerical integration is the pacing challenge.
I Often highly parallel supercomputers are required to obtain sufficiently

high-precision results.
I By comparison, running PSLQ to identify numerical values is usually very cheap.

Thus, more efficient algorithms are needed to evaluate the following to high precision:
I Basic arithmetic operations and transcendental functions.
I A wide range of special functions.
I Ordinary one-dimensional integrals.
I Multi-dimensional integrals (this is particularly expensive).

High-precision evaluation of transcendental and special functions
1. Basic transcendentals—exp, log, sin, cos, tan, hyperbolic functions—and thecorresponding inverse functions.
2. Gamma, digamma, polygamma, incomplete gamma, beta and incomplete beta
functions.
3. Riemann zeta function, polylogarithms and Dirichlet L-functions.
4. Bessel functions (first, second and third kinds, modified, etc.).
5. Hypergeometric functions.
6. Airy functions.
7. Elliptic integral functions.
8. Jacobian elliptic functions and Weierstrass elliptic/modular functions.
9. Theta functions.
Also, derivatives are needed for all of the above.

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.

A fast series to compute π2 (x, y )The following series, due to the late Richard Crandall, can be used to rapidly compute
π2 (x, y ) to high precision:
φ2 (x, y ) =1cosh(πx) + cos(πy ) 2 X cosh(πmx) cos(πmy )
log
−
.
4π
cosh(πx) − cos(πy ) πm(1 + e πm )+m∈O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α

A conjectured formula for the degree
Jason Kimberley has observed empirically, based on our computations to date, that thedegree m(d) for the polynomial corresponding to φ2 (1/d, 1/d) appears to be given for
odd primes by m(4k + 1) = (2k) · (2k) and m(4k + 3) = (2k + 2) · (2k + 1).
If we set m(2) = 1/2, for notational convenience, then it seems that for any prime
factorization of an integer greater than 2:
!
k
Y
Q
?
2(e −1)mpiei= 4k−1 ki=1 pi im(pi ).
i=1

To determine whether or not this formula holds, even larger computer runs (and evenhigher precision) will be required.

Lessons from the Poisson lattice project
I Finding a fast series for φ2 (x, y ) was a critical step in this research.
I This problem is a good example where the PSLQ computation, not the initial

numerical evaluation, is the pacing challenge.
I This problem is also a good example of the need for extremely high precision —50,000 to 100,000 digits.
I We wanted to run larger examples, but run times were unreasonable (severaldays), and in one case the run failed due to an obscure bug in the ARPREC
package — “very reliable” isn’t good enough.
I A faster, highly parallel version of PSLQ is needed for very large problems such asthese.

Analysis of Mordell-Tornheim-Witten sums
In three recent papers, the present authors and colleagues (including the late Richard
Crandall) explored sums such asω(s1 , . . . , sK +1 )

:=

X

1

P

QM

sK
ms1 · · · mK
(m1 + · · · + mK )sK +1
m1 ,...,mK > 0 1

and, more generally,ω(s1 , . . . , sM | t1 , . . . , tN )

:=

1
= 2π

m1 ,...,mM ,n1 ,...,nN > 0
PM
P
m = N
n
j=1 j
k=1 k

R 2π QM
0

j=1 Lisj

e iθ

 QN

Here the polylogarithm of order s is denoted by Lis (z) :=We also explored character variants of these sums.

1 QN
1
j=1 mj sj
k=1 nk tk
k=1 Litk

P

n≥1 z


e −iθ dθ.

n /ns .

Computation of Mordell-Tornheim-Witten sums
To compute these sums to high precision required high-precision evaluations of:
I A new scheme for computing derivatives of polylogarithms with respect to theorder for real and complex arguments.
I The ζ function and its derivatives at integer arguments (positive and negative).
I Bernoulli numbers.
I The Γ function and derivatives.
I Definite integrals of complex functions.
and much more.
For full details, see:
I David H. Bailey, Jonathan M. Borwein and Richard E. Crandall, “Computation and theory of extended
Mordell-Tornheim-Witten sums,” Mathematics of Computation, vol. 83, no. 288 (Jul 2014), pg.
1795–1821.
I David H. Bailey and Jonathan M. Borwein, “Computation and theory of extended Mordell-TornheimWitten sums II,” 05 May 2014, available at http://www.davidhbailey.com/dhbpapers/mtw2.pdf.
I David H. Bailey and Jonathan M. Borwein, “Computation and structure of character polylogarithms withapplications to Mordell-Tornheim-Witten sums,” Mathematics of Computation, to appear, 20 Feb 2014,
available at http://www.davidhbailey.com/dhbpapers/mtw3.pdf.

Summary
I Very high-precision arithmetic is an increasingly important component of modern

experimental mathematics and mathematical physics research.
I Precision levels of several hundred digits are commonplace; some problems requireas much as 50,000 to 100,000 digits.
I Highly parallel computer implementations are often required for pacing problems.
I Faster and even easier-to-use open-source software is needed for high-precision
arithmetic, supporting computations with tens or hundreds of thousands of digits
and a wide range of functions.
I Research is needed into faster algorithms for very high-precision evaluation of
integrals (especially multi-dimensional integrals).
I Research is needed into new, faster algorithms for very high-precision evaluationof special functions.
I Research is also needed into highly parallel implementations of PSLQ.
This talk is available at http://www.davidhbailey.com/dhbtalks/dhb-icerm2014.pdf.