High-precision arithmetic and PSLQDavid H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Lab (retired)
University of California, Davis

October 2, 2017

1 / 36

Numerical reproducibility in scientific computing
A December 2012 workshop on reproducibility in computing, held at Brown University,
USA, noted that
Science is built upon the foundations of theory and experiment validated andimproved through open, transparent communication. ...
Numerical round-off error and numerical differences are greatly magnified ascomputational simulations are scaled up to run on highly parallel systems. As
a result, it is increasingly difficult to determine whether a code has beencorrectly ported to a new system, because computational results quickly
diverge from standard benchmark cases. And it is doubly difficult for otherresearchers, using independently written codes and distinct computer
systems, to reproduce published results.
I V. Stodden, D. H. Bailey, J. M. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting thedefault to reproducible: Reproducibility in computational and experimental mathematics,”
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.
2 / 36

Numerical reliability
Many applications routinely use either 32-bit or 64-bit IEEE arithmetic, and employ
fairly simple algorithms, assuming that all is well. But problems can arise:
1. Large-scale, highly parallel simulations, running on systems with hundreds ofthousands or millions of processors — numerical sensitivities are greatly magnified.
2. Certain applications with highly ill-conditioned linear systems.
3. Large summations, especially those involving +/- terms and cancellations.
4. Long-time, iterative simulations (such as molecular dynamics or climate models).
5. Computations to resolve small-scale phenomena.
6. Studies in computational physics and experimental mathematics.
I D. H. Bailey, R. Barrio, and J. M. Borwein, “High precision computation: Mathematical physics

and dynamics,” Applied Mathematics and Computation, vol. 218 (2012), pg. 10106–10121.

3 / 36

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?
4 / 36

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

5 / 36

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 U.C. Berkeley graduates who will do technicalcomputing in their careers have had rigorous training in numerical analysis.
6 / 36

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

I If this computation is done with ordinary double

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

2.5

2.0

1.5

1.0

0.5

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

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

0.5

1.0

1.5

2.0

2.5

3.0

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

7 / 36

Aren’t 64 bits enough?
There has been considerable resistance in the scientific computing community to theusage of high-precision arithmetic. Why?
I Many are persuaded that physical reality fundamentally does not require highprecision, beyond, say, the limits of 64-bit IEEE arithmetic.
I Many believe that numerical sensitivity problems are always due to the usage ofinferior algorithms.
I Some regard the usage of high-precision arithmetic as “cheating” or “sinful”(because it is too easy?).
I Some complain that easy-to-use high-precision arithmetic software facilities are

not widely available, or the computational cost is too great.
But none of these excuses are valid.
I D. H. Bailey, R. Barrio, and J. M. Borwein, “High precision computation: Mathematical physics

and dynamics,” Applied Mathematics and Computation, vol. 218 (2012), pg. 10106–10121.
8 / 36

Innocuous example where standard 64-bit precision is inadequate
Problem: Find a polynomial to fit the data (1, 1048579, 16777489,
84941299, 268501249, 655751251, 1360635409, 2523398179, 4311748609) forarguments 0, 1, · · · , 8. The usual approach is to solve the linear system:






Pn n

k=1 xk..
Pn .

n
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 

 
=





..
.
..
  ..  

.
.
Pn
Pn
n
2n
an
k=1 xk yk
k=1 xk

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
9 / 36

Innocuous example, cont.
The result on the previous page can be obtained with double precision using Lagrange
interpolation or the Demmel-Koev algorithm. But few scientists, outside of expert
numerical analysts, are aware of these schemes.
Besides, even these schemes fail for higher-degree problems.
For example: (1, 134217731, 8589938753, 97845255883, 549772595201,
2097396156251, 6264239146561, 15804422886323, 35253091827713,
71611233653971, 135217729000001, 240913322581691, 409688091758593)
is generated by:f (x ) = 1 + 134217729x 6 + x 12 = 1 + (227 + 1)x 6 + x 12
Lagrange interpolation and Demmel-Koev fail for this problem, but a straightforward
Linpack scheme, implemented with double-double arithmetic, works fine.
10 / 36

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. Problems in experimental mathematics (100–50,000 digits).
11 / 36

Coulomb n-body atomic system simulations
I Alexei Frolov of Queen’s University in Canada has used high-precision arithmetic

to solve a generalized eigenvalue problem that arises in Coulomb n-bodyinteractions.
I Matrices are typically 5,000 x 5,000 and are very nearly singular.
I Computations typically involve massive cancellation, and high-precision arithmetic

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

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

Frolov: “We can consider and solve the bound state few-body problems ... beyond ourimagination even four years ago.”
I 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.
12 / 36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.

13 / 36What is needed for high-precision floating-point software?
1. A rock-solid-reliable arithmetic engine, with precision scalable to 1,000,000 digits
or more — 1,000 or even 10,000 digits is no longer enough.
2. Highly efficient algorithms, so single-processor performance is as fast as possible,
given constraints of portability and ease of installation.
3. A thread-safe design to facilitate parallel processing, both thread-based (low-level)
and message-passing (high-level) parallelism.
4. A comprehensive library of tuned transcendentals and special functions:sin, cos, exp, gamma, incomplete gamma, polylogarithms, Bessel functions, etc.
5. High-level language interfaces for C++, Fortran-90 and other languages.
6. Interoperability with Maple and Mathematica.

14 / 36

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. MPFUN90: Similar to ARPREC, but is written entirely in Fortran-90 and provides
only a Fortran-90 interface. 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.
15 / 36

How do multiprecision packages operate?
I A multiprecision floating-point number is represented as a string of DP floats orlong integers, typically 53 or 64 bits per word, with the first few words holding
sign, exponent and working precision (see below).
I Addition and subtraction are performed using relatively straightforward extensionsof familiar schemes.
I Multiplication is performed using FFTs (see below).
I Division is performed using Newton iterations, with a level of working precisionthat approximately doubles with each iteration.
I Square roots are performed using Newton iterations, similar to division.
I Transcendental functions (cos, sin, exp, log, etc.) are performed using series

expansions, after transforming the input into a small argument range. For higher
precision levels, quadratically convergent iterations are used.

16 / 36

How to perform high-precision multiplication using FFTs
Let A = (a0 , a1 , a2 , · · · , an−1 ) and B = (b0 , b1 , b2 , · · · , bn−1 ), where each element
stores successive bits, and extend A and B to length 2n by padding with zeroes. Let
C = A × B. Then except for releasing carries,
Ck =

2n−1
X

aj bk−j ,j=0where we add 2n to the b subscript if it is negative. Then
C = (C0 , C1 , C2 , · · · , C2n−2 )= (a0 b0 , a0 b1 + a1 b0 , a0 b2 + a1 b1 + a2 b0 , · · · , an−1 bn−1 ).
This is merely the convolution of A and B, which can be computed using FFTs:
C = F −1 (F (A) · F (B)),
where (F (A) · F (B)) means the element-wise product of the sequences F (A) and F (B).
For sufficiently large precision levels, FFT-based multiplication is many times fasterthan conventional grade-school multiplication — 5n log2 n instead of n2 operations.

17 / 36Thread safety: A major challenge for arbitrary precision software
Since run times are greatly magnified by high-precision arithmetic, parallel processing isincreasingly essential.
However, almost all existing arbitrary-precision packages are not thread-safe at theapplication user level, and this kills thread-based parallelization
I Almost all high-level packages use global variables, and these variables are deadly

for thread safety.
I Most global variables can be dispensed with, but not the working precision level,

which is changed frequently within the package itself and often by users also.
I The only reasonable design for user programming is via operator overloading,available in languages such as Fortran-90 and C++, but overloading does not
permit the working precision level to be passed.
What can be done?
18 / 36

MPFUN2015: A thread-safe arbitrary precision package
DHB has written a package (approx. 50,000 lines of code) for arbitrary precision
floating-point computation. It is available in two versions:I MPFUN-Fort: Completely self-contained, all-Fortran version. Compilation is a
simple one-line command, which completes in a few seconds.
I MPFUN-MPFR: Calls the MPFR package for lower-level operations. Installation
is significantly more complicated (since GMP and MPFR must first be installed),
but performance is roughly 3X faster. We used MPFUN-MPFR in in the Poisson
polynomial study.
Both versions include a high-level language interface, using Fortran custom datatypes
and operator overloading — for most applications, only a few minor changes toexisting double precision code are required.
A C++ version is planned.
Full details of design, algorithms, installation and usage are given in
I D. H. Bailey, “MPFUN2015: A thread-safe arbitrary precision computation package,”

http://www.davidhbailey.com/dhbpapers/mpfun2015.pdf.
19 / 36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 in physics.
–Ð Kurt Godel

20 / 36

The PSLQ integer relation algorithm
Given a vector (xn ) of real numbers, an integer relation algorithm finds integers (an )such thata1 x1 + a2 x2 + · · · + an xn = 0(to within the precision of the arithmetic being used), or else finds bounds within
which no relation can exist.
Helaman Ferguson’s PSLQ algorithm is the most widely used integer relation
algorithm, although variants of the LLL algorithm can also be used.
Integer relation detection requires very high numeric precision, both in the input dataand in the operation of the algorithm — at least nd digits, where d is the size in digitsof the largest |ai |.
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), 351–369.
2. D. H. Bailey and D. J. Broadhurst, “Parallel integer relation detection: Techniques andapplications,” Mathematics of Computation, vol. 70, no. 236 (Oct 2000), 1719–1736.
21 / 36

Helaman Ferguson’s “Umbilic Torus SC” sculpture at Stony Brook Univ.

22 / 36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 precision in digits).
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.
23 / 36

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

24 / 36

Application of multipair PSLQ
One simple but important application of multipair PSLQ is to recognize a computed
numerical value as the root of an integer polynomial of degree m.
Example: The following constant is suspected to be an algebraic number:
α = 1.232688913061443445331472869611255647068988824547930576057634684778 . . .
What is its minimal polynomial?
Method: Compute the vector (1, α, α2 , · · · , αm ) for m = 30, then input this vector tomultipair PSLQ.
Answer (using 250-digit arithmetic):
0 = 697 − 1440α − 20520α2 − 98280α3 − 102060α4 − 1458α5 + 80α6 − 43920α7+ 538380α8 − 336420α9 + 1215α10 − 80α12 − 56160α13 − 135540α14 − 540α15+ 40α18 − 7380α19 + 135α20 − 10α24 − 18α25 + α30
25 / 36Apéry-like summations
Apéry’s proof of the irrationality of ζ(3) has prompted considerable effort to extend
these results to larger integer arguments. The formulas in question areζ(2) = 3

∞
X

1


2 2k
k=1 k
k

,

(1)

ζ(3) =∞
5X(−1)k+1
 ,2 k=1 k 3 2kk(2)

ζ(4) =∞
36 X1
.
17 k=1 k 4 2kk

(3)

Is there an analogous formula with a rational coefficient for ζ(5) or higher? A PSLQ
computation with 10,000 digits ruled this out. What about multi-term relations?
26 / 36

Some multi-term relations found by PSLQ
ζ(5) = 2∞∞
k−1
X(−1)k+1
5 X (−1)k+1 X 1
−
,


2
j2
k 5 2kk 3 2k
k
k
j=1
k=1
k=1
∞

ζ(7) =

ζ(9) =∞

k−1

25 X (−1)k+1 X 1
5 X (−1)k+1
 +

2k
2
2
j4
k7 k
k 3 2k
k
j=1
9
4

k=1

k=1

∞
X

∞
X

k−1
k−1
∞
X
5
(−1)k+1(−1)k+1 X 1(−1)k+1 X 1
−
+
5



4
j2
j4
k 9 2kk 7 2kk 5 2k
k
k
k
j=1
j=1
k=1
k=1
k=1

∞
k−1
∞
k−1
k−1
45 X (−1)k+1 X 1
25 X (−1)k+1 X 1 X 1
−
,
+


4
j6
4
j 4 j=1 j 2
k 3 2kk 3 2k
k
k
j=1
j=1
k=1

k=1

∞
∞
k−1
5 X (−1)k+1
25 X (−1)k+1 X 1
ζ(11) =
 +

2
2
j4
k 11 2k
k 7 2k
k
k
j=1
k=1
k=1

−

∞
k−1
∞
k−1
k−1
75 X (−1)k+1 X 1
125 X (−1)k+1 X 1 X 1
+
.


4
j8
4
j 4 i=1 i 4
k 3 2kk 3 2k
k
k
j=1
j=1
k=1

k=1

27 / 36

General formulas
Using bootstrapping and an application of the “Pade” function (which in both
Mathematica and Maple produces Padé approximations to a rational function satisfied
by a truncated power series), the following general results were discovered:∞k−1
Y 1 + 4x 4 /m4
5X(−1)k+11=4
4
k 3 (1 − x 4 /k 4 )2 k=1 k 3 2k
1 − x 4 /m4
k=1
k (1 − x /k ) m=1
∞
X

∞
X

!

k−1
∞
Y 1 − 4 x 2 /m2
X
1
1
=
3

2 2k (1 − x 2 /k 2 )k2 − x2
1 − x 2 /m2
m=1
k=1
k=1 k
k

!

Power series expansions for a specific x yield all the results for ζ(n).
The second formula was proven by the Wilf-Zeilberger method, a computer-based
scheme that yields a rigorous proof. It was thus discovered and proven by computer.
I D. H. Bailey and J. M. Borwein, “Computer-assisted discovery and proof,” T. Amdeberhan and

V. H. Moll, ed., Tapas in Experimental Mathematics, AMS, May 2008, pg. 21–52.
28 / 36

The tanh-sinh numerical integration algorithm
Given f (x ) defined on (−1, 1), define g(t) = tanh(π/2 · sinh t). Then setting x = g(t)yieldsZ ∞Z 1

f (x ) dx−1=−∞f (g(t))g 0 (t) dt ≈ hNX

wj f (xj ),j=−Nwhere xj = g(hj ) and wj = g 0 (hj ) can be precomputed for a given h. Since g 0 (t) goesto zero very rapidly for large t, the integrand f (g(t))g 0 (t) typically is a nicebell-shaped function, so that the simple summation above converges very rapidly.
Reducing h by half typically doubles the number 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).
1. D. H. Bailey, X. S. Li and K. Jeyabalan, “A comparison of three high-precision quadratureschemes,” Experimental Mathematics, vol. 14 (2005), no. 3, pg. 317–329.
2. H. Takahasi and M. Mori, “Double exponential formulas for numerical integration,” Publicationsof RIMS, Kyoto University, vol. 9 (1974), pg. 721-Ð741.

29 / 36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:
Cn :=

Dn :=

4
n!

Z ∞
0

0

4
n!

Z ∞

Z ∞

=

···

Z ∞

0

Z 1

En

···

2
0

0

···

Z 1
0



1
P

n
j=1 (uj + 1/uj )

Q

i<j



ui −ujui +uj

2

dun
du1
···
u1
un

2

du1
dun
···
u1
un

2

P

n
j=1 (uj + 1/uj )2uk − uj 
dt2 dt3 · · · dtnu+ uj
1≤j<k≤n kYwhere 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.
30 / 36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):
Cn =

2n
n!

Z ∞
0

tK0n (t) dt

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 https://isc.carma.newcastle.edu.au. The result was:lim Cn = 2e −2γ .

n→∞where γ denotes Euler’s constant. This is now proven.
31 / 36Other Ising integral evaluations found using PSLQ
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.
E5 was reduced to a 3-D integral of a 60-line integrand, which was evaluated usingtanh-sinh quadrature to 250-digit arithmetic using over 1000 CPU-hours on a highly
parallel system. The PSLQ calculation required only seconds.

32 / 36

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
φ2 (x , y ) =1 X cos(mπx ) cos(nπy )
1
=
log α
2
2
2
π m,n oddm +nπ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.
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,
http://www.davidhbailey.com/dhbpapers/PoissonLattice.pdf.
33 / 36

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 (see next viewgraph). This computation
required 10,000-digit precision. Most of the time was spent on PSLQ, not quadrature.
Other polynomials required up to 50,000-digit precision, with runs extending to several
days. The ARPREC software failed on at least one of these problems.
34 / 36

192-degree minimal polynomial found by multipair PSLQ for x = y = 1/35
−1−127936α1−119729568α2+86259615424α3−63414593420752α4−23761356230294720α5−21521006149475920864α6+611939502140287079872α7−368241444434585386649704α8+34534137997535622208800832α9+190487328667825240861480928α10+70731866698900460686569086656α11−30639321248717091789168538879600α12+1922209242084456061370070376570688α13−92225937246234731111918315818791008α14+5237085972458845348795423777776516544α15−235187168740813424983476512954481139732α16+7429223399159195923652315325576416851776α17−183888805409443551177519325744878574722592α18+3271185405201792963719366162166966763587520α19−30299858585212592953207274837044089721625936α20−359992131877582010695076712447360628122541504α21+21102025339314733846086441916632636325491195040α22−476414747523697562008419279360098519444810134848α23+7169682338305684610208768903330082567075749930088α24−77147038334570163492295981291240215533152479746240α25+576164594448741453064385715717563185909190291963744α26−2736845608824484767533372475869861068967680432337984α27+14259726273826198867906144815986191950918103600120720α28−418898418536517746458185180042147268305073686192184768α29+11576727630588055309375415098495616370526328495842993184α30−215864800528660839585002162871715414654518484882221347136α31+2990075428515164844805618172618994337159659201108555443326α32−32103405822691698883721869657536188329165809065032597536320α33+266563559238616545234045835785848715794979563412588549450912α34−1588318769762074127932309766905276050093368042694667151097024α35+4203983024865006644081516487096952711041894782888827508128720α36+45741795572596990194281792252600295928845217432176402329342656α37−888295322846459251495779433538199608179522746791003115766467360α38+9318868074863097497435367536227145321575550333899160360356887616α39−75360568030730663716613169159158032198635749984420852664531432312α40+511904470165868318936871843389229757000581029769096726546021748160α41−3023546278817700158319248501909046229060364677626955289951128897760α42+15762439699463291788277556458704672972787957041647009977980115608384α43−72736254087237112123433235647299134668422519990232924303557659578512α44+294655977143150502457165877012542972911987116693855121083332478262976α45−1023418666617067140855492177628022043641503497647927192890895515496096α46+2872077837131433597858520147550093715024860106120898508958083768560704α47−5244656757384067630568130031166724732476128944651080774800327815250116α48−4375176099824059357126327245955265067024929918868045662842672259365696α49+106846263928061290720995820179791755195031418846948198250527939644422688α50−723639781157123885224054475278025279616831219096501715148796353565656000α51+3859388734169523941501510320228124865152568559427283813598867448036672464α52−18252524309522593867333335673661626604394606878814374964432136569382704704α53+77801314121940709085266334072671847312436136451041737155601244216316832608α54−295924487894698517430363594146438746974018698578472622597771256623091267264α55+990883094386701201740101120896676716687418747972363663196338999871139087736α56−2890647897108739627601089028825831836951153432467718962671967596342224331584α57+7329280206578701233923536747523603867159738577214393279268225554697992024224α58−16529926761499923756105727653427646107756611125773189972311741253023019544512α59+36621401257350141587688985844214102107297318891842289523223570545646432107248α60−97806177157646206090474078202543089426709440278475346399840243186928760712768α61+342578593578669077931375617617876157438332517966938868387799492392630391951328α62−1318494126999138496197980784859753020673610723076044612133782854663456341738176α63+4770648039063802379141989053619421884321933507524681278243277292208295535890193α64−15472388710953109453763060955215809455238864792962040591046439233238251454944640α65+44993126133327412718394996675533962755303579158603905782013384675151509701621696α66−118993678023270686331414029749780909390537125214711043931901528442137124714048640α67+291344068959030996949609862561213601556661749401236700666696596438800158224020448α68−672232544304160020140857460151229781783982323755372371953704546852887742475693952α69+1483756546102154371991621621590988762258290163435879224071206411473630010421819712α70−3160698702072288781771809020419553277725945343183442108817290538424291765729105536α71+6501409535318122589047227869235783252522832570490016457500770280864513874652673648α72−12819483008841743167047686498254296089188806357669061135254105501043843161878338944α73+23940036962388573943673383046497610566473660443396196928270109685668366770350448320α74−41783013070841149166516454269070382975961508352312427230333483156459522180289158272α75+67368007706334464267463678880196649459723255934263643619439168617608299192825564576α76−99642265418236291841000783588241107943576928164973329187609122090423042369846900608α77+135684188062925983721829145899247186881926375911313411751554577363782314702542167104α78−175027359690635875118927274192015444148378662148930449526142186378361018115616958080α79+230610470370976022179942433042240414747465002830049591030281645554186668837721744600α80−345422259636540072058441736406800364089516876786410123401855696720824391150388817792α81+607137279243399399687877325853919130478767608151821751113771985567328949196202362560α82−1146510481188778693690480353627945407015450534575302104964676472156118460868902533760α83+2105965918280150700239337266375899616661414836881685735099986426462071874731653407456α84−3578099460146102835319749769961500032509484008579871514018505538256220639466406781312α85+5536176932749319486998056607166537426841616046655786532661950492192459385736417266240α86−7794809387121374412027723614750591530039194538680900530684292991362557543590881311872α87+10030481053829943456941341746197226060710869001059273241750964485149858086288306297744α88−11855047237622313521872161223741536919023583101964392808122345150250641280790373298048α89+12893688641739677210122894074867960717441055755861579189125047363100403706385492504000α90−12809768930934167317548440121172865343713268055679115994897043173034752329114006909568α91+11262021213965304198804174698825443038955997004185306033241807189003803427064173678496α92−7849339137067651884501051001175657872891640580392350339448446525146908046759472284032α93+2136636160531235808486132100072698791511127607341553178572108601522645635713639494976α94+6182245041445709382582906443600895624739514751663138744452121338925118726793099119488α95−17027139906672640329386998091617550180769534392263392331607751724460799383582736580380α96+29712689442226074490778272413199981836245382131442863784077108073804247980419140795264α97−42913200159533980876837680464079550150161925164315385234076771249584133995790605789888α98+54875995472202404064433213065393198278520963443598737304851154803060885419026758577792α99−63825924492815000947804006660643589720936270079124038361274887833762135234355101020256α100+68414285135348430004939705729313406303651902455163077326777079330134700391931405532544α101−68052822555902228583831595467842429584162408591782046566558459654662872517242957132352α102+63030283713395534150979634221684031444907507016820847420753457009176282665669129552000α103−54395011700122510436311965134146069935113985222440557617340304628640177217303223510128α104+43661984233468332856484141573706612495184170864395287152161473759679906830742322508672α105−32445083300862179709263952078603686506241016652233926920050485448240544236200991135168α106+22121488216739201475462682897890849680076071306863200873172222890471168791244013948544α107−13610378673217079648813094175609528224090224497348403954492135388377641815307093302560α108+7303046170572758067914716975092412590455319569015290620012702916095079596345712424320α109−3131236977889234295288739605819294961076022655796155501587485430182196310350741412160α110+722365801671568238376537172560503547318051932349906390028021898906118213151314346112α111+423601672347880485375916854724736111845460547572098224965700862087464288437477713112α112−789372219071689131389195080081110795264640291084303894471284509727586944366644459136α113+754844732770124186279971574591381270567329047796107073522275693379041299799000185920α114−570721609484121135843526600643872181106278042023747642137438316501500109706499901312α115+373431689502164436293456395379326584419834315995264527851770397586202342215482658208α116−218828293376534306771048774311447660676991404850097946000689151980961553792808627328α117+116819901747981791438375631794807731777629646639313589196533853783698354327272619712α118−57453490833130202085276802289813871676642923245392648322689802787108386720310373760α119+26322973910217009898036214368526156227892184032219090112763726491148743638265345648α120−11422293344043226356444271083191575214531430352803712260273887252165192023955482240α121+4829637193522288826261580288388799608086659826584129673254995173846268064699275584α122−2078817645839949613653484866883184944970522163048142353896259321595643918566571904α123+954520484158500566039633963780071095170152665054125308901179320476061925020693472α124−476059705123351801546660809193814948416900226947841461260532426685030594421013632α125+251019218316831093023612611636151648630754206670792757985712506316495845522315200α126−133387943415348693213086926089912666295139757189579549272926455197093747072689536α127+68669026491757170599503679381653507376133436777604611890688641155842031359807249α128−33465058649746756993954154200626776375881202657693549363818927924447316311308992α129+15265533880702809827068316634690758620560027313985769068567445431127044039098336α130−6486555236218719687502565717159307954250235769707821904032108839967986971592256α131+2562655019815708913497740003785614466295686123950263933555489655843860264778480α132−940718343979230237299701844566444332546715633645270453853660861734418393188288α133+320771113247818349988983708385653420090152866910698288499686276305286454845600α134−101564452363828199943698157849015312966252040558153110340955250960209509456704α135+29840961316238615781884908158117270537546621262752881740317094715665305372024α136−8126376541879784121657187449503319934776774553449947624311536879983067197120α137+2046994405582919610833914133280166599368341899820908062960382811295175875424α138−475319952147035704834648908024038830280982173827054926876842801020099391040α139+101156139413765995787243855087479738205336119604403640417144416030882137040α140−19536944183496787665897200701914355830359619924729440672204518600726485952α141+3365653573540320310956912862417397989082549665698949749718751351520029216α142−500297582644104844305032452251740515847130656744730715102464280225674048α143+59315729762438733681847579101255419553245420619721207125629069846855484α144−4088968402271371424704252106469344820514085294011337888699549616483264α145−410467411100906673223801525015092827907382852275953565568989408591520α146+244373061239573525170364306606773136690590544856331578133540062357184α147−67137173574208581929099251980225280832815704879993996789961828208784α148+15109453726979262114404719422671376896934991181254480301068783496000α149−3041110741874930638456938010940447684323532404585513230577319188704α150+543474006539496836413675511283771054220572953580018501137140774336α151−82818210860123209370721479679382853454718857015790579069608342392α152+10208909053669116781213384172905303479740054958245095883752727616α153−950607878263673381203873433978784840581837837196127719161731872α154+61537081222427663468738382440653874375219870006734261147376320α155−3614460685557548015719619196028506524834155796764598110372912α156+735146337727640447184879802116259476500819273333226229041984α157−195982957610545601975537049907586449965158403571666750972768α158+35410924134603994553881822516846688340985279903764142254528α159−4586047122379261183418802995335871699189992574849634423682α160+461731528371721725878159689755046846967144056271095432896α161−39548982563052670609179355793571460666409257926095752160α162+3249771246267792941295436011847663684436436914874939968α163−278663520870823657308884684973631236730880421051628656α164+23298555433192575099637973614016576299882162116142016α165−1572480211244215449798675303409379604821388457521312α166+64631732188205970561818215170229051834055602691904α167−359811432336185820581657539275402497742419221400α168−81408350985197555231086181640295965142859955520α169−2670526194482509900914512598125562096574018400α170+824346510881081434287176385644420962607200832α171−40421524575366405147127382177082176049208656α172+859850260877732200189610333732782641323968α173−21746942237477507610278314875487367745056α174+1166204976072113679141136387334878097216α175−16255593447525276321168587153278001172α176−1031047179999272558707644762054830656α177+57315125130135149120405695909830560α178−2190911023714361858443237127579840α179+42156348386945234230569584900496α180−400797347138034992145950828864α181−3239264180104926121034061856α182+57518923520963556655145024α183−302326293499892000458856α184−648961315190292320832α185+4804670128512328736α186+97238301302930240α187−106616047418832α188−478181605696α189+208298592α190+25664α191−α192This distinctive semi-elliptical pattern is strong evidence that the result is correct.

35 / 36

Summary

I High-precision arithmetic can now be performed rapidly by means of widely

available, easy-to-use software.
I High-precision arithmetic is required by the PSLQ algorithm, which finds integer

relations among a set of computed floating-point real numbers.
I This combined capability (HP arithmetic and PSLQ) is a tool of enormous power

for experimental mathematics.
Thanks!
This talk is available at
http://www.davidhbailey.com/dhbtalks/dhb-hiprec-2017.pdf

36 / 36