Petascale Computing and High-Precision Arithmetic:
Applications and Challenges
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Lab (retired) and University of California, Davis
Co-author: Jonathan M. Borwein (University of Newcastle, Australia)

Petascale supercomputers: Measured in quadrillions (1015 )The memory capacity of a supercomputer is measured in “Pbytes” or “petabytes”(one Pbyte = one quadrillion bytes).
The performance of a supercomputer is measured in “Pflop/s” or “petaflops”(one Pflop/s = one quadrillion 64-bit floating-point operations per second).
I The volume of Lake Zurich is 3.9 × 1012 litres, or roughly 1/256 quadrillion litres.
I The volume of Lake Tahoe, California is 1.51 × 1014 litres, or roughly 1/7

quadrillion litres.
I The distance to Alpha Centauri is 4.15 × 1016 metres, or 41.5 quadrillion metres.
I The time elapsed since the big bang is 7.24 × 1015 minutes, or 7.24 quadrillion

minutes.

Increasing performance of the top 500 supercomputers (1994 – present)

Orange: #1
Blue: #500
Green: Sum of #1 thru #500

I “Performance development,” Top500.org, available at http://top500.org.

Systems at the Lawrence Berkeley National Laboratory’s NERSC facility

Current system (“Edison”):

To be installed in 2016 (“Cori”):

I 5576 Intel “Ivy Bridge” nodes.

I 9300 Intel “Knight’s Landing” nodes.

I Memory per node: 64 Gbyte.

I Main memory per node: 64 Gbyte.

I Performance per node: 460.8 Gflop/s.

I Performance per node: 3 Tflop/s.

I Total main memory: 357 Tbyte.

I Total main memory: 595 Tbyte.

I Total performance: 2.57 Pflop/s.

I Total peak performance: 28 Pflop/s.

I Total disk storage: 7.56 Pbyte.

I Total disk capacity: 28 Pbyte.

The accelerating universe
In 1998, two teams of astronomers (one led
by Saul Perlmutter of LBNL, and the otherled by Brian P. Schmidt of Australian
National University), came to theparadoxical conclusion that the expansionof the universe is accelerating, not slowing
down.
Both teams based their results on carefulmeasurements of distant supernovas, which
in turn were found by sifting through reams
of digital telescope data. The U.S. team inparticular relied heavily on large computer
systems at the NERSC facility, coupledwith a worldwide network of collaborating
astronomers.

Discovery of supernova in the Pinwheel Galaxy
In 2011, Peter Nugent of LBNL, working
within a worldwide consortium of fellow
astronomers, discovered a Type Ia
supernova in the Pinwheel Galaxy, which is“only” 21 million light-years from earth. It
is the closest and brightest supernova of
this type seen in the last 30 years.
Nugent and his team utilized the Palomar
Transient Facility (PTF), a
robot-controlled telescope that produces
digital images, plus a worldwide network ofdata processing and storage facilities.

The Palomar Transient Facility (PDF) data pipeline

The Materials Project
Researchers at LBNL and Harvard
University started the “Materials Project”:
I Reduce lag time between materials

science advances and real-worldcommercialization.
I Invert the conventional materials

science paradigm: Ask “Whatproperties do I want?,” then ”Whichmaterials have them?”
I Method: Perform ab initio calculationsof tens of thousands of potentiallyinteresting compounds, then save insearchable database.
I Now in operation:

http://www.materialsproject.org.

Today’s Li ion
batteries
Interestingmaterials…

!

Voltage vs. capacity for over 20,000 Li-ion
cathode compounds using high-throughput ab
initio methods.

DANGER AHEAD: Reproducibility in scientific computing
A December 2012 workshop on reproducibility in computing, held at Brown University
in Rhode Island, USA, identified these issues:
I The need to carefully document the full context of computational
experiments—system environment, input data, code used, computed results, etc.
I The need to save the code and output data in a permanent repository.
I The need for reviewers, research institutions and funding agencies to recognize theimportance of computing and computing professionals, and to allocate funding forafter-the-grant support and repositories.
I The need to encourage publication of negative results—other researchers canoften learn from them.
I The need to ensure responsible reporting of performance.
I The increasing importance of numerical reproducibility, and the need for tools toensure and enhance numerical reliability.
I V. Stodden, D. H. Bailey, J. M. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default toreproducible: Reproducibility in computational and experimental mathematics,” manuscript, 2 Feb 2013,
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

Reproducibility in biomedicine
The biomedical field has been stung by numerous cases where pharma products lookgood based on clinical trials, but later disappoint in real-world usage, or the resultscannot be reproduced in separate studies. Examples:
I In 2004, GlaxoSmithKline acknowledged that while some trials of Paxil found it

effective for depression in children, other unpublished studies showed no benefit.
I In 2011, Bayer researchers reported that they were able to reproduce the results of

only 17 of 67 published studies they examined.
I In 2012, Amgen researchers reported that they were able to reproduce the resultsof only 6 of 53 published cancer studies.
I In 2014, a review of Tamiflu found that while it made flu symptoms disappear a

bit sooner, it did not stop serious complications or keep people out of the hospital.
These experiences have exposed a fundamental flaw in methodology:
Only publicizing the results of successful trials introduces a bias into the results.
The AllTrials movement would require all results to be public: http://www.alltrials.net

Reproducibility in physics
In March 2014, a team of researchers from Harvard
University made the dramatic announcement thatthey had discovered an interesting “twisting”
pattern in cosmic microwave background data,measured using their BICEP2 experimental system.
This pattern fit very well with the hypothesizedpattern of the most commonly assumed model of“inflation” in the first tiny fraction of a second after
the big bang, and thus has been trumpeted as thefirst experimental evidence of the inflationary
cosmology.
But other researchers had difficulty reconstructing
the claimed results. Finally, two teams challenged
the BICEP2 findings, saying that the results could
more readily be explained by dust in the Milky Way.

I Ron Cowen, “Doubt grows about
gravitational waves detection,”
Scientific American, 2 Jun 2014.

Numerical reproducibility
The ICERM reproducibility report 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,”
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

Numerical reproducibility, continued
Particularly vulnerable are:
1. Large-scale, highly parallel simulations, running on systems with hundreds ofthousands or millions of processors — numerical sensitivities are greatly magnified.
2. Certain applications with highly ill-conditioned linear systems.
3. Large summations, especially those involving cancellations.
4. Long-time, iterative simulations (such as molecular dynamics or climate models).
5. Computations to resolve small-scale phenomena.
6. Studies in computational physics or experimental mathematics often require hugeprecision levels.

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?

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

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

Enrollment in numerical analysis courses:
I Math 128A (Introductory numerical analysis required of math majors): 219.
I Math 128B (A more advanced course, required to do serious work): 24.

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

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

I If this computation is done with ordinary double

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

2.5

2.0

1.5

1.0

0.5

accuracy), it takes 47.39 seconds seconds and yields theresult 7.073157029007832.
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, π).

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.

Aren’t 64 bits enough?
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

Aren’t 64 bits enough? 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—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.

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. GFORTRAN: Now provides full REAL*16 (IEEE 128-bit, or 34-digit) support.
3. GMP: Produced by a volunteer effort and distributed under the GNU license.
http://gmplib.org.
4. MPFR: C library for multiple-precision floating-point computations with exact
rounding, based on GMP. http://www.mpfr.org.
5. MPFR++: High-level C++ interface to MPFR.
http://perso.ens-lyon.fr/nathalie.revol/software.html.
6. MPFUN2015: Similar to ARPREC, but is written entirely in Fortran-90. A 100%
thread-safe design. http://crd.lbl.gov/~dhbailey/mpdist.
7. 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).

I C. Rubio-Gonzalez, C. Nguyen, H. D. Nguyen, J. Demmel, W. Kahan, K. Sen, D. H. Bailey, C. Iancu andD. Hough, “Precimonious: Tuning assistant for floating-point precision,” Proceedings of SC13, 17 Nov
2013.

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. Optimization problems in biology and other fields (32 digits).
5. Coulomb n-body atomic system simulations (32–120 digits).
6. Schrodinger solutions for lithium and helium atoms (32 digits).
7. Electromagnetic scattering theory (32–100 digits).
8. Scattering amplitudes of fundamental particles (32 digits).
9. Discrete dynamical systems (32 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).

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.
I G. Lake, T. Quinn and D. C. Richardson, “From Sir Isaac to the Sloan survey: Calculating the structure
and chaos due to gravity in the universe,” Proc. of the 8th ACM-SIAM Symp. on Discrete Algorithms,
1997, pg. 1–10.

Supernova simulations
I Researchers at LBNL have used quad-double

arithmetic to solve for non-local thermodynamic
equilibrium populations of iron and other atoms in
the atmospheres of supernovas.
I Iron may exist in several species, so it is necessary

to solve for all species simultaneously.
I Since the relative population of any state from thedominant state is proportional to the exponential of
the ionization energy, the dynamic range of thesevalues can be very large.
I The quad-double portion now dominates the entire

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

Climate modeling: High-precision for reproducibility
I Climate and weather simulations are fundamentallychaotic: if microscopic changes are made to thecurrent 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 forcode maintenance.
I He and Ding found that by using double-double

arithmetic in two key inner loops, most of thisnumerical variation disappeared.
I 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.

Optimization solutions in biochemistry
I Michael Saunders (a well-known optimization researcher) and his graduate

students have encountered significant numerical difficulties when trying to
produce numerically reliable results using the Modular In-core Nonlinear
Optimization System (MINOS) software.
I In many applications, final results are only produced to a few significant digits,

due to numerical error.
I For one biochemistry problem in particular, Harvard researchers have asserted thatan exact rational arithmetic scheme was the only means to produce reliable results.
I When the MINOS software was converted to use gfortran’s new REAL*16 (quadprecision) facility, almost all of these troublesome applications (including thebiochemistry problem) now produce strong, numerically reliable results.
I D. Ma and M. A. Saunders, “Solving multiscale linear programs using the simplex method in quadruple
precision,” to appear in M. Al-Baali, L. Grandinetti and A. Purnama, ed., Recent Developments in
Numerical Analysis and Optimization, Springer, 2014.

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

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

matrix in four- and five-body systems.
I These computations typically require 120-digit arithmetic.

Frolov: “We can consider and solve the bound state few-body problems ... beyond ourimagination even four years ago.”
I A. M. Frolov and D. H. Bailey, “Highly accurate evaluation of the few-body auxiliary functions andfour-body integrals,” Journal of Physics B, vol. 36, no. 9 (14 May 2003), pg. 1857–1867.

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 TIDES with 300 digits, versus 1 time period using DP.

Experimental mathematics: Discovering new mathematical results
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.

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.

High-precision numerical integration: The tanh-sinh algorithm
Given f (x) defined on (−1, 1), define g (t) = tanh(π/2 sinh t). Then
Z 1
NX
f (x) dx ≈ h
wj f (xj ),
−1

j=−N

where xj = g (hj ) and wj = g 0 (hj ).
Features:
I Reducing h by half typically doubles the number of correct digits.
I Works well even for functions with singularities at the endpoints.
I Generation of abscissas and weights is much faster than Gaussian quadrature.
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).
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 applied our methods to study three classes of integrals: Cn are connected toquantum field theory, Dn arise in the Ising theory of mathematical physics, while the
En integrands are derived from Dn:
Z
Z ∞
du1
1
dun
4 ∞
···
···
Cn :=
2
P
n! 0
u
un
n
1
0
j=1 (uj + 1/uj )Q  ui −uj 2
Z ∞Z ∞i<j ui +uj
4
dun
du1
Dn :=
···
···
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )

2
Z 1
Z 1
Y uk − uj dt2 dt3 · · · dtn
En = 2···
uk + uj
0
0
1≤j<k≤nwhere in the last line uk = t1 t2 · · · tk .
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. E5 remained a“numerical conjecture” for several years, but was proven in March 2014 by Erik Panzer.

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.

Algebraic numbers in Poisson potential functions
Lattice sums arising from the Poisson equation have been studied widely inmathematical physics and image processing. We numerically discovered, and thenproved, that for rational (x, y ), the 2-D Poisson potential function satisfies
1 X cos(mπx) cos(nπy )1log αφ2 (x, y ) ==π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: Dedicatedto 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 found for k = 32MPFUN2015: A NEW thread-safe arbitrary precision package
I 100% thread safe. NO global read/write global variables.
I NO initialization required (unless over 20,000 digits required).
I NO system-dependent features or settings — the code only assumes 64-bit IEEE

format (it does not depend on rounding mode).
I NO reliance on sophisticated language features that might not be supported.
I Installation is a snap — one compile line.
I Precision is scalable to over billion digits; largest size is over 1064,000,000 .
I A full package of basic arithmetic functions. Higher-level routines forbinary-decimal conversion, I/O and transcendentals.
I A full-featured high-level language interface for Fortran-90 is provided, so thatmost users need only make minor changes to existing double precision code.
I Extensive error checking and a solution to the double precision accuracy problem.
I Currently only Fortran-90 version; C++ version is being developed in Australia.
I Available at: http://www.davidhbailey.com.

For further reading:
1. D. H. Bailey and J. M. Borwein, “High-precision arithmetic in mathematical physics,”
manuscript, 16 Jan 2015.
http://www.davidhbailey.com/dhbpapers/hp-math-physics.pdf.
2. D. H. Bailey, R. Barrio and J. 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.
3. David H. Bailey, “MPFUN2015: A thread-safe arbitrary precision package,” manuscript,
13 Mar 2015. http://www.davidhbailey.com/dhbpapers/mpfun2015.pdf.

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