Experimental mathematics: Serious research and great fun
David H. Bailey
Lawrence Berkeley National Laboratory (retired) andUniversity of California, Davis
http://www.davidhbailey.com/

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

1 / 38Experimental math makes mathematics fun again
Using computation to:
I Explore mathematical questions for insight and intuition.
I Visually represent mathematical objects.
I Test (and falsify) mathematical conjectures.
I Discover new potential identities and theorems.
I Help find rigorous proofs of discovered conjectures.
1970s: “Real mathematicians don’t compute.”
2020s: Most mathematicians now use computer tools.
All major historical mathematicians (Newton, Euler, Gauss, Hilbert, Ramanujan,
others) used experimental methods, although they lacked our modern computer tools.
Although computer-based tools are very useful in mathematical exploration, facts
discovered by experimental methods must still be proven the old-fashioned way.
2 / 38

An interdisciplinary endeavor
I Research mathematicians.
I Applied mathematicians, physicists and chemists.
I Computer scientists.
I Numerical analysts.
I Data scientists and statisticians.
I Parallel computing experts.
I Computer graphics and visualization experts.
Persons of all levels can make valuable contributions:
I Senior research mathematicians.
I Graduate students and post-docs.
I Undergraduate mathematics and computer science majors.
I High school students with good math and computer programming skills.
I have personally worked with several students (including a former Euler Circle TA)
who have made major contributions, both mathematically and computationally.

3 / 38

Tools of experimental mathematics
I Symbolic math software, such as Mathematica, Maple and Sage.
I Online facilities such as the Online Encyclopedia of Integer Sequences
http://www.oeis.org and the Inverse Symbolic Calculator
http://wayback.cecm.sfu.ca/projects/ISC/ISCmain.html.
I Graphics and visualization software, such as Cinderella and Penrose.
I Parallel computing software, such as OpenMP and MPI.
I High-precision arithmetic software packages, such as GMP and MPFR.
I Custom-written software (typically C, C++, Fortran, Java or Python) forgraphics, matrix computation, factorization, primality testing, high-precision
numerical integration, integer relation detection and other tasks.

4 / 38

A brief tour of modern experimental mathematics
In this talk, I will briefly sketch some illustrative examples of real state-of-the-art
experimental mathematics, rooted in relatively basic concepts:
I Algebra and trigonometry.
I Basic number theory — number bases, primes, modular arithmetic.
I Basic calculus — derivatives, integrals, infinite series.
I Basic computer science — algorithms, floating-point arithmetic.
However, of necessity many details are skipped. For full details on any of these topics,see online references in orange or send me an email.
Either way, please don’t freak out if you don’t immediately follow every detail.

5 / 38

Normal numbers
A real number α is said to be normal base b for integer b ≥ 2 if and only if for every
integer m ≥ 1, every m-long string of base-b digits appears with limiting frequency
1/b m (i.e., with the frequency expected “at random”).
It can be shown using measure theory that almost all real numbers are normal base b,
for any (or all) b. However, it is not easy to list any particular number that is provably
normal. Here is one that is provably normal base 2 (or base 16):
α=

∞
X
k=0

1
3k 23k

= 0.541883680831502985071252898624571682426096758 . . .

(decimal)= 0.8ab8e38f 684bda12f 684bf 35ba781948b0fcd6e9e0652 . . .

(base 16)

1. R. Stoneham, “On absolute normality in the rational fractions with applications to normal numbers,”
Acta Arithmetica, vol. 22 (1973), 277–286.
2. D. H. Bailey and R. E. Crandall, “Random generators and normal numbers,” Experimental Mathematics,
vol. 11 (2002), 527–546, https://www.davidhbailey.com/dhbpapers/bcnormal-em.pdf.
3. D. H. Bailey and M. Misiurewicz, “A strong hot spot theorem,” Proceedings of the American
Mathematical Society, vol. 134 (2006), 2495–2501,
https://www.davidhbailey.com/dhbpapers/dhb-mm-hotspot.pdf.

6 / 38Is π a normal number?
π=3.14159265358979323846264338327950288419716939937510582097494459230781
6406286208998628034825342117067982148086513282306647093844609550582231
7253594081284811174502841027019385211055596446229489549303819644288109. . .
It is not yet known whether or not π is normal base b for any b, although researcherssuspect π is normal base b for all b. Computational experimentation to date based ontrillions of decimal digits and corresponding numbers of base-16 digits have not shown
any statistical abnormalities.
Such computations cannot prove that π is normal, but if they disclosed that π is likely
not normal in some base b, that would be a very interesting result.
Similarquestions apply for most other common constants of mathematics, including√
2, log 2, e, γ (Euler’s constant), φ (golden ratio) and many others.
Conjecture: Every irrational algebraic number is normal base b for all b.
7 / 38

Visualization: A random walk on the first 100 billion base-4 digits of π

Sequence: red-orange-yellow-green-cyan-blue-purple-red. This dataset can be explored at:
http://gigapan.com/gigapans/106803
1. F. J. Aragon Artacho, D. H. Bailey, J. M. Borwein and P. B. Borwein, “Walking on real numbers,”
Mathematical Intelligencer, vol. 35 (2013), 42–60,
https://www.davidhbailey.com/dhbpapers/tools-walk.pdf.

8 / 38

The PSLQ integer relation algorithm
Given a real vector x = (xk ), an integer relation algorithm finds integers (ak ) such thata1 x1 + a2 x2 + · · · + an xn = 0
to within the “epsilon” of the numeric precision being used.
These algorithms can discover mathematical identities, based only on numerical valuesof the terms xk and conjectured forms of the terms involved, provided that very highnumerical precision (hundreds or thousands of digits) is employed.
The “PSLQ” algorithm of mathematician-sculptor Helaman Ferguson is the mostwidely used integer relation algorithm, although the “LLL” algorithm is also used.

1. H. R. P. Ferguson, D. H. Bailey and S. Arno, “Analysis of PSLQ, an integer relation finding algorithm,”
Mathematics of Computation, vol. 68 (Jan 1999), 351–369.
2. D. H. Bailey and D. J. Broadhurst, “Parallel integer relation detection: Techniques and applications,”
Mathematics of Computation, vol. 70 (Oct 2000), 1719–1736,
https://www.davidhbailey.com/dhbpapers/ppslq.pdf.
9 / 38

How PSLQ operates
I PSLQ constructs a sequence Bn of integer-valued matrices that reduce the size ofthe entries of the vector y = x · Bn , until either the relation is found (as one of
the columns of 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 numeric precision in digits).
The size of this drop is a “confidence level” that the relation is a true relation.
I The “multipair variant” of PSLQ greatly reduces the number of iterations
required. “Multilevel” implementations greatly reduce run time by using multiple
levels of precision, doing as much computation as possible in low precision.
I The state-of-the-art is “four-level multipair PSLQ.”
Relations discovered using PSLQ must still be proven the old-fashioned way, but it istypically easier to find a proof if one already “knows” that the result is likely correct.

10 / 38

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

Note the abrupt drop at iteration 199, by nearly 200 orders of magnitude.

11 / 38

The binary algorithm for exponentiation
Problem: What is 317 mod 10 (i.e., what is the last decimal digit of 317 )?
Algorithm A:
317 = 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 = 129140163,
so answer = 3.
Algorithm B (faster): 317 = ((((32 )2 )2 )2 ) × 3 = 129140163, so answer = 3.
Algorithm C (fastest):
317 = ((((32 mod 10)2 mod 10)2 mod 10)2 mod 10) × 3 mod 10 = 3.
Another example: 317,777,777 mod 100 = 63 is easy using Algorithm C.
Note that in Algorithm C mod n, we never need to do arithmetic with integers larger
than (n − 1)2 .

12 / 38

How to compute arbitrary binary digits of log 2
Consider this well-known formula for log 2 (due to Euler):log 2 =∞
X1
= 0.6931471805599 . . .10 = 0.101100010111001000010111 . . .2
k 2k
k=1

In 1996, Peter Borwein observed that the binary digits of log 2 beginning after positiond can be written as {2d log 2}, where {·} denotes fractional part, so one can write
( d
) ( ∞
)
X 2d−k
X 2d−k
{2d log 2} =
+
k
k
k=1
k=d+1
( d
) ( ∞
)
X 2d−k mod k
X 2d−k
=
+
,
k
k
k=1

k=d+1

where “modk” has been inserted since we are only interested in the fractional partwhen divided by k. The numerator 2d−k mod k can be calculated very rapidly using
the binary algorithm for exponentiation, using standard computer arithmetic.
Is there a similar formula and algorithm for π? Prior to 1996, none was known.

13 / 38

The BBP formula for π
In 1996, Simon Plouffe used a PSLQ-type program to discover this new formula for π:


∞
X
4
1
2
1
1
π=
−
−
−
16k 8k + 1 8k + 4 8k + 5 8k + 6
k=0This formula permits one to directly compute a string of binary or base-16 digits of π
beginning at an arbitrary starting position, using a simple algorithm similar to log 2.
Example computation: The two quadrillionth binary digit of π is 0.
This result also has implications for the normality of π (see reference #2 below).

1. D. H. Bailey, P. B. Borwein and S. Plouffe, “On the rapid computation of various polylogarithmic
constants,” Mathematics of Computation, vol. 66 (Apr 1997), 903–913,
https://www.davidhbailey.com/dhbpapers/digits.pdf.
2. D. H. Bailey and R. E. Crandall, “On the random character of fundamental constant expansions,”
Experimental Mathematics, vol. 10 (Jun 2001), 175–190,
https://www.davidhbailey.com/dhbpapers/baicran.pdf.

14 / 38

Some other BBP-type formulas discovered using PSLQ
∞


216
72
54
9
144
−
−
−
+
(6k + 1)2
(6k + 2)2
(6k + 3)2
(6k + 4)2
(6k + 5)2
k=0

∞
405
81
27
2 X 1
243
−
−
−
π2 =
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
4096
8192
26112
15360
1024
π log(2) =
−
−
+
−
256
4096k (24k + 1)2(24k + 2)2(24k + 3)2(24k + 4)2(24k + 5)2
π2 =

1X 1
8
64k



k=0

9984
11520
2368
512
768
+
+
+
−
+
(24k + 6)2(24k + 8)2(24k + 9)2(24k + 10)2
(24k + 12)2
64
408
720
16
196
−
+
+
+
+
(24k + 13)2
(24k + 15)2
(24k + 16)2
(24k + 17)2
(24k + 18)2

37
60
+
−
(24k + 20)2
(24k + 21)2
1. D. H. Bailey, J. M. Borwein, A. Mattingly and G. Wightwick, “The computation of previously inaccessible
digits of π 2 and Catalan’s constant,” Notices of the AMS, vol. 60 (2013), 844–854,
15 / 38https://www.davidhbailey.com/dhbpapers/bbp-bluegene.pdf.

Using PSLQ to find the minimal polynomial of an algebraic number
One simple but important application of PSLQ is to recognize a computed numericalvalue as the root of a polynomial of degree m with integer coefficients.
Example: α is suspected to be an algebraic number of degree 30 or less:
α = 1.232688913061443445331472869611255647068988824547930576057634684778 . . .
What is the minimal polynomial satisfied by α?
Method: Compute the vector (1, α, α2 , · · · , α30 ) to 250-digit precision, then input thisvector to a multipair PSLQ program.
Answer: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
Mathematica verified that this polynomial is irreducible.

16 / 38

Using PSLQ and high-precision numerical integration to identify integrals
Another common application of PSLQ is to identify the formula for a definite integral,
given only its high-precision numerical value.
Problem: Most numerical integration techniques do not work well for hundreds or
thousands of digits, particularly for functions with singularities.
For example, Gaussian quadrature (look this up if you are not familiar) is impracticaldue to the very high cost of computing abscissas and weights at high precision.
Curious but important fact: For a few bell-shaped functions,such as the Gaussian probability function
12
f (x) = √ · e −t /2 ,2π
a simple block or trapezoidal rule approximation to theintegral works remarkably well — reducing the block size hin half roughly doubles the number of correct digits. Why?
Answer: The Euler-Maclaurin summation formula.

17 / 38

The Euler-Maclaurin summation formulaZ b
f (x) dx = h
a

n
X
j=0

where

m
X

h
h2i B2if (xj ) − (f (a) + f (b)) −D 2i−1 f (b) − D 2i−1 f (a) − E (h),
2
(2i)!
i=1

2m+2

|E (h)| ≤ 2(b − 1)(h/(2π))

max |D 2m+2 f (x)|.

a≤x≤bHere h = (b − a)/n and xj = a + jh; B2i are Bernoulli numbers (look this up if you are
not familiar); and D k f (x) is the k-th derivative of f (x).
In effect, the E-M formula provides high-order correction terms to the trapezoidal rule.
Note that in the case when f (x) and all of its derivatives are zero at the endpoints (asfor some bell-shaped curves), the error E (h) of the trapezoidal approximation goes tozero more rapidly than any power of h.
Can we apply this fact to evaluate more general definite integrals?
18 / 38

The hyperbolic trigonometric functions
The hyperbolic trigonometric functions are defined as:
e x − e −x2
e x + e −xcosh(x) =2
e x − e −xsinh(x)= xtanh(x) =cosh(x)e + e −xsinh(x) =Red: cosh; Blue: sinh

Their properties are analogous to standard trigonometric functions, for example:
D[sinh(x)] = cosh(x)D[cosh(x)] = sinh(x)D[tanh(x)] = sech2 (x) = 1/ cosh2 (x)There are numerous applications in pure and applied mathematics.
19 / 38

The tanh-sinh numerical integration algorithm
Given h > 0, define x = g (t) = tanh(π/2 · sinh t) and compute, for −N ≤ j ≤ N,
xj = g (hj) = tanh(π/2 · sinh(hj)), wj = g 0 (hj) = π/2 · cosh(hj)/ cosh(π/2 · sinh(hj))2 .
Then for a given integrable C ∞ function f (x) on (−1, 1),
Z 1
Z ∞
N
X
0
f (x) dx =f (g (t))g (t) dt ≈ h
wj f (xj ).
−1

−∞

j=−N

It looks wacky, but it works great:
I Reducing h by half typically doubles the number of correct digits.
I It works well even for many functions with blow-up singularities at endpoints.
I Computing abscissas (xj ) and weights (wj ) is much faster than other methods;
the (xj ) and (wj ) values can be precomputed and used for different functions.
In effect, the substitution x = tanh(π/2 · sinh t) converts almost any C ∞ function f (x)on (−1, 1) into a bell-shaped function for which the Euler-Maclaurin trick works.
1. D. H. Bailey, X. S. Li and K. Jeyabalan, “A comparison of three high-precision quadrature schemes,”
Experimental Mathematics, vol. 14 (2005), 317–329,
https://www.davidhbailey.com/dhbpapers/quadrature.pdf.

20 / 38

Curious sinc integrals
A colleague of Jonathan Borwein (my long-time collaborator, deceased 2016)
discovered the following curious fact using computations [here sinc(x) = sin(x)/x]:Z ∞sinc(x) dx =

Q0 =
Z0 ∞

π
,
2

sinc(x) sinc(x/3) dx =

Q1 =
Z0 ∞
Q2 =

π
,
2

sinc(x) sinc(x/3) sinc(x/5) dx =
Z0 ∞

π
,
2

sinc(x) sinc(x/3) · · · sinc(x/13) dx =

Q6 =
0

···
π
,
2

Z ∞
sinc(x) sinc(x/3) · · · sinc(x/13) sinc(x/15) dx

Q7 =
0

=

467807924713440738696537864469 π= 0.499999999992647 . . . π.
935615849440640907310521750000

Why does the pattern abruptly change beginning at Q7 ?
1. D. H. Bailey, J. M. Borwein, V. Kapoor and E. Weisstein, “Ten problems in experimental mathematics,”
American Mathematical Monthly, vol. 113 (Jun 2006), 481–409,
https://www.davidhbailey.com/dhbpapers/tenproblems.pdf.

21 / 38

Mystery solved[Don’t worry if you don’t follow the details]
R Q
∞

n

Define Tn = 0
k=0 sinc(ak x) dx, for positive reals (ak ), a0 = 1. Since sinc(ak x) = F [fk ](x),where fk (x) = 1/(2ak ) · χ(−ak , ak ) and F denotes the Fourier transform, it follows by the
convolution theorem that Tn = F (f0 ? f1 ? · · · ? fn ), and by the Fourier inversion theorem,Z ∞Z 1
Tn (x) = 2π(f0 ? f1 ? · · · ? fn )(0) = π(f1 ? f2 ? · · · ? fn )(u) du.
−∞−1

SincePthe support of fk is [−ak , ak ], the support of (f1 ? f2 ? · · · ? fn ) is [−s, s], wherens = k=1 ak . Now note that if s ≤ 1, then
Z 1nY
(f1 ? f2 ? · · · ? fn )(u) du = F (f1 ? f2 ? · · · ? fn )(0) =sinc(ak 0) = 1,−1

k=1

so that Tn = π/2. On the other hand, if s > 1, then the interval [−1, 1] is strictly inside thesupport of (f1 ? f2 ? · · · ? fn ), and one can conclude thatZ 1
(f1 ? f2 ? · · · ? fn )(x) dx < 1,−1

so that Tn < π/2. This explains why Qn abruptly drops at n = 7:
1/3 + 1/5 + · · · + 1/13 < 1, but 1/3 + 1/5 + · · · + 1/15 > 1.
22 / 38

Strange cosine integral
The following integral arose in analysis of random walks:
Z ∞∞Y
cos(2x)cos(x/n) dx
0

n=1

The infinite product converges slowly, but after converting it to a sum and some additional
manipulation, the integral can be calculated using the tanh-sinh scheme. Oddly, the result
agrees with π/8 to 42 decimal digits, but is not equal to π/8. Why?
Z ∞∞Y
cos(2x)cos(x/n) dx = 0.39269908169872415480783042290993786052464543418 . . .
0

n=1

π
= 0.39269908169872415480783042290993786052464617492 . . .
8Richard Crandall showed that the integral is the first term of a rapidly convergent series for π/8:
∞ Z ∞
∞
X
Y
π
=
cos[2(2m + 1)x]cos(x/n) dx
8
m=0 0
n=1
One term is accurate to 42 digits; two terms are accurate to 350 digits; etc.
1. R. E. Crandall, “Theory of ROOF walks,” Algorithmic reflections: Selected works, 2012,https://www.reed.edu/physics/faculty/crandall/papers/ROOF12.pdf.

23 / 38

Three integrals from mathematical physics
Cn are connected to quantum field theory, Dn arise in the Ising theory of mathematicalphysics, while the En integrands are derived from Dn :
Z
Z ∞
4 ∞
1
dun
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
du1
dun
Dn =
···
···
P
2
n! 0
u
un
n
1
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 . Can we study these integrals computationally?
1. 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), 12271–12302,
https://www.davidhbailey.com/dhbpapers/ising.pdf.
24 / 38

Limiting value of Cn : What is this number?
We were able to convert the Cn multi-dimensional integrals integrals into 1-D integrals:
Z
2n ∞
Cn =
t K0n (t) dt,
n! 0
where K0 (t) is the modified Bessel function of the second kind (look it up).
When we evaluated this formula using tanh-sinh, we found that Cn approach a limit:
Cn → C∞ = 0.6304735033743867961220401927108789043545870787 . . .
What is this number? We copied the first 50 digits into the online Inverse Symbolic
Calculator (ISC) at http://wayback.cecm.sfu.ca/projects/ISC/ISCmain.html. The
result was:C∞ = 2e −2γ .
where γ = 0.5772156649 . . . denotes Euler’s constant. We subsequently were able toprove this and numerous other related identities.

25 / 38

Other Ising identities discovered using tanh-sinh and 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 Lk is a Dirichlet L series value, Lik is the polylog function, and ζ(x) denotes the
Riemann zeta function (don’t worry if you are not familiar with these).
Both the computations and the proofs were quite challenging!
1. 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), 12271–12302,
https://www.davidhbailey.com/dhbpapers/ising.pdf.
26 / 38

Polynomials related to the Poisson potential function
In 2012 Richard Crandall found that each pixel in the image produced by a sharpening
algorithm is given by the 2-D Poisson potential function of mathematical physics:
1 X cos(mπx) cos(nπy )φ2 (x, y ) = 2
π
m2 + n2m,n odd

Crandall, Jon Borwein and others subsequently proved that for rational x and y ,1φ2 (x, y ) = log α,πwhere α is algebraic, i.e., the root of a some polynomial with integer coefficients.
√
When x = y = 1/4, α = 5.8284271247 . . . = 3 + 2 2, a root of 1 − 6α + α2 .
When x = y = 1/5, α = 13.627503029 . . ., a root of 1 + 52α − 26α2 − 12α3 + α4 .
Problem: The above series converges much too slowly for high-precision computation.
1. D. H. Bailey, J. M. Borwein, J. Kimberley and W. Ladd, “Computer discovery and analysis of large
Poisson polynomials,” Experimental Mathematics, 27 Aug 2016, vol. 26, 349–363,
https://www.davidhbailey.com/dhbpapers/poisson-res.pdf.
27 / 38

A rapidly convergent formula
A key breakthrough, due to Jon Borwein, was to show that the φ2 (x, y ) function canbe written in a form much more conducive to rapid computation:
θ2 (z, q)θ4 (z, q)
1,
φ2 (x, y ) =log
2π
θ1 (z, q)θ3 (z, q)
where q = e −π = 0.043213918 . . . and z = π2 (y + ix). The four theta functions can becomputed using these rapidly convergent formulas (note the exponents of q):∞
X2
θ1 (z, q) = 2(−1)k−1 q (2k−1) /4 sin((2k − 1)z),
θ2 (z, q) = 2

k=1
∞
X

2

q (2k−1) /4 cos((2k − 1)z),k=1

θ3 (z, q) = 1 + 2∞
X2

q k cos(2kz),

k=1
∞
X
2
θ4 (z, q) = 1 + 2(−1)k q k cos(2kz).
k=1

28 / 38

Some initial results
Given x and y , we used the above formulas to compute α to high precision and thenused the multipair PSLQ algorithm to find the minimal polynomial for α:
s
5
6
7
8
9

10

Minimal polynomial for α corresponding to x = y = 1/s
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

I What is the connection between x and y and the degree of α?
I Note that for even s, the polynomials are palindromic — coefficients are the same
backwards as forwards. Does this pattern hold for larger s?
29 / 38

Kimberley’s formula for the degree
Jason Kimberley (a student at University of Newcastle, Australia) observed that thedegree m(s) of the minimal polynomial of α for x = y = 1/s is given by this formula:Set m(2) = 1/2; otherwise for primes p congruent to 1 mod 4, set m(p) = (p − 1)2 /4,
and for primes p congruent to 3 mod 4, set m(p) = (p 2 − 1)/4. Then for any otherpositive integer s whose prime factorization is s = p1e1 p2e2 · · · prer ,m(s) = 4r −1

rY2(ei −1)pim(pi ).

i=1

Subsequently, in 2016, Watson Ladd (a student at U.C. Berkeley) proved that
Kimberley’s formula is correct for the case x = y = 1/s, and that the palindromic
property holds for even s.
Does Kimberley’s formula hold for pairs of rationals other than x = y = 1/s?
1. D. H. Bailey, J. M. Borwein, J. Kimberley and W. Ladd, “Computer discovery and analysis of large
Poisson polynomials,” Experimental Mathematics, 27 Aug 2016, vol. 26, 349–363,
https://www.davidhbailey.com/dhbpapers/poisson-res.pdf.

30 / 38

High-level computational algorithm
1. Given rationals x = p/s and y = q/s, select a polynomial degree m and precisionlevels P1 and P2 (for our computer runs, P2 varies from 1,000 to 64,000 digits).
2. Calculate φ2 (x, y ) to P2 -digit precision using the fast formulas on page 28.
3. Calculate α = exp(8πφ2 (x, y )) and the (m + 1)-long vector
X = (1, α, α2 , · · · , αm ) to P2 -digit precision.
4. Apply a four-level multipair PSLQ algorithm to find an integer relation for X .
5. If a numerically significant relation is not found, try again with a larger degree mor a higher precision level P2 . If a relation is found, use Mathematica or Maple toensure that the resulting polynomial is irreducible.
Code size: 25,000 lines (using MPFUN2020) or 140,000 lines (using GMP/MPFR).

1. D. H. Bailey, “Finding large Poisson polynomials using four-level variable precision,” Correctness 2021:
5th International Workshop on Software Correctness for HPC Applications, 19 Nov 2021,
https://www.davidhbailey.com/dhbpapers/correct-2021.pdf.
31 / 38

Nitya Mani’s palindromic polynomial shortcut for even s
In all cases studied to date when s is even, the the corresponding minimal polynomial
a0 + a1 α + a2 α2 + · · · am αm is palindromic (i.e, aj = am−j for each j ≤ m).
For example, when x = y = 1/8, the minimal polynomial is1 − 88α + 92α2 − 872α3 + 1990α4 − 872α5 + 92α6 − 88α7 + α8 .
Nitya Mani, a former Euler Circle TA, pointed out that if α satisfies a palindromic
polynomial, then β = α + 1/α satisfies a polynomial of half the degree, and theoriginal polynomial can be easily reconstructed from the half-degree polynomial.
Since the cost of performing the multipair PSLQ algorithm increases very rapidly withthe degree m, using this fact dramatically reduces the cost for handling even s cases.

1. N. Mani, “Palindrome degree reduction,” 19 Dec 2015, available athttps://www.davidhbailey.com/dhbpapers/mani-palindrome.pdf.
32 / 38

Sample result: Minimal polynomial for case x = y = 1/23
−1α0
1
+276α
−48990α23+6243028α−268081353α4−7007943232α56+699168895168α−8507658502848α7−228421295042000α89+5752931816856576α−69577168674471200α1011+1342467622983507840α−20886025588167398672α12
13+349167061432874646208α−3452961624838699374016α1415+15640784440670709844544α−104217914365387271401240α16−707297164619952782619488α17−1010779866294321116930832α18−38938027282982631315388896α1920+15273847812560162975303784α−252294875204439839252885056α21−1814499584808466663818882112α2223+7122860530259481805651039552α−29751952082955733043119937840α24
25+69112629323019561967486699392α26+277026903082879050281512565408α−1750707356769870221042718523648α2728+7841746142835495601991495671440α−6369728528150841895179072801088α29−35953899240983246926617945898688α30
31+253492876456503829747447461470784α−667834023356310407949079857214780α32

33+162422492508454890608937648198512α
34+3931664673190042230498787905679544α−18133691752545610684537913799858192α3536+23789518577406441345286599055526948α37+37469246645313996740590585412048064α−281341968248747952781979853250623040α3839+685450730380390572422936104954217024α−144807877454540334657454671748536144α40−2676374507247534118825617453738202880α4142+9268558185390732238617127106511230176α−8169131934494953200462056091539832192α43−24973838730281646460618615909661528464α4445+86480873940520233375619792198402276800α−135145074117747229141553700760222342848α46−136878295961379583616822852976503569600α4748+678907009350034701692481661741570615448α−826551987822128918951623497620329820704α49−250029772766791506678263063944048615280α50
51+4203025326736104693470282522553895192928α−3878919678334634848760927707284752064872α52
53+162350544202944626514876598545210959040α
54+4823912376464138143603163479868453706944α−32447159140179954019440351777350445058496α55
56+10793210639006798965181292630236352423632α
57+47798645963669738247513432401008300599936α−53267062459588215249236047567693700327008α58
59+75013990401402134185716434599423648706560α
60+64458486492326422733823977462792116677520α−215173808790139348785993658139108405629504α6162+85620600429649929065268861189687330676800α
63+139350732807546237327599795099241600771904α−250813976164701012680304124927666981877222α6465+54309279545871700785244094943285407971768α

66+191983097542468140804868975448582309836108α−106896043545407401874530608602325810772552α67−142528636680982846937656120402002070771798α68
69+73854370040353117163017414423037785534784α
70+140156083674583453941489740970247416329792α−69130044528962858268002108127048139470912α71−71616687080402664876669283116574668250864α72
73+36784745627135753478698634095086817636864α
74+35394737651242262546479699690525861027744α−17956071634537519607161038659310053274496α75−20336190903346145735503496223269203568816α76
77+13630911262027292414127427055252676137024α
78+5572729792117138939086526506626923795136α−4089865120803933198217666901071614415168α79−2615092808097454015621055806574050460968α80
81+2209795446144060278513925847082003764064α
82+536972872828811706650471125454348811920α−838590764716382848440202241096095828000α83
84−10999416202903176320023415762860858792α
85+179376397229534413041533561964022758720α−15609601791658126093273272542469773504α86−85912879928606509932720682745289619520α87
88+35900403296202961761181993065193949424α
89+11889588476079748335308838031449347200α−14452120597436246105835100451405878560α9091−3928001780218000982842803096677131520α
92+702309189177306708681044564957381680α−900667846828975766802897215449019328α93−1070482587263099407527831014338023488α9495+347733729184112119976507013855153856α
96+91316928183854286074743342756853892α−53432493225801914733293041393576464α9798+29297268341789315496992698093351608α

99+7690044340064802181369044817359216α−3188354423894290085673630211864476α100−58750370497093670583364169980864α101102+183505586931227063763291585126208α−93640045058313908718602479122240α103−25202486755374245498242332804080α104105+7119062987888570394024596150528α
106+1030356560081795611678995024032α−393078180816678972203397758592α107−8326519137079696206002948272α108109+21913274492964339461296318272α
110+105999845651185538221855168α−863563217740448639591640128α111−48073238371689968208546776α112
113+18565893213139990900332576α
114+2105838095925307433972464α−225795463065220784914784α115−50162337393304568016088α116117+447902847210305133632α
118+773461151618304399424α
119+53770166932871224512α−4202551956967976592α120−844977474669598336α121−36912569450921760α122123+2411687976601088α
124+371577683528112α
125+16571351703360α
126+1279246016α−26017245248α127−694105569α128
129+7591380α
130+522146α131+6164α132+23αThe distinctive crescent shape of this set of coefficients is strong numerical evidence
that all software worked correctly and this is the true minimal polynomial for this case.
33 / 38

Results: Degrees of minimal polynomials for x = 1/s; y = q/s
s
10
11
12
13
14
15
16
17
18
19
20
21
22
23

1
11
8
30
16
36
24
32
32
64
36
90
32
96
60
132
132

2
12
16
30
32
36
48
32
64
64
72
90
64
96
120
132

q: first row 1–10; second row, 11–20
3
4
5
6
7
8
9
13
14
15
16
17
18
19
4
16
30
30
30
16
32
8
36
36
18
36
24
48
24
48
32
16
32
32
32
32
64
32
64
16
64
32
64
64
64
64
36
72
36
72
36
72
90
90
90
90
90
90
90
32
64
32
64
32
64
16
96
96
96
96
96
48
96
60 12060 12060 12060
132 132 132 132 132 132 132

10
20

96
120
132
34 / 38

Results: Degrees of minimal polynomials for x = 1/s; y = q/s
s
24
25
26
27
28
29
30
31

1
11
64
32
100
100
72
72
162
162
96
96
196
196
64
32
240
240

2
12
128
100
100
144
144
162
162
192
192
196
98
128
128
240
240

q: first row 1–10; second row, 11–20
3
4
5
6
7
8
13
14
15
16
17
18
64 12832 12832 128

9
19
64

10
20
128

100

100

100

100

50

100

100

100

72

144

36

144

72

144

72

144

162
162
96
48
196
196
64
64
240
240

162

162

162

162

162

162

162

192

96

192

96

192

96

192

196
196
128
128
240
240

196

196

196

196

196

196

64

128

64

128

64

128

240
240

240

240

240

240

240
35 / 38

Results: Degrees of minimal polynomials for x = 1/s; y = q/s
s
32
33
34
35
36
38
39
40

1
11
128
128
240
240
128
128
192
192
144
144
180
180
288
288
128
64

2
12
256
256
240
240
256
256
192
192
288
288
360
360
288
288
256
256

q: first row 1–10; second row, 11–20
3
4
5
6
7
8
13
14
15
16
17
18
128 256 128 256 128 256
128 25664
240 240 240 240 240 240
240 240 240 240
128 256 128 256 128 256
64 256 128 256
192 192 192
96 192 192
192 192 192 192 192
144 288 144 288 144 288
144 288 144 288
72
180 360 180 360 180 360
180 360 180 360 180 360
288 288 288 288 288 288
288 144 288 288 288 288
128 256 128 256 128 256
128 256 128 256 128 256

9
19
128

10
20
256

240

120

128

256

192

192

144

288

180

360

288
288
64
64

288
256
36 / 38

Latest research findings
Note that many of these results are not consistent with Kimberley’s formula, which has
been proven only for the cases x = y = 1/s for integer s.
For example, note that when s is even, the degrees for even q are, in most cases,
double the degrees for q = 1. But there are exceptions. For instance when s = 36,
degrees alternate between 144 and 288, yet when q = 17, the degree is 72.
There are also anomalies when s is odd. Note, for instance, that when s = 35, alldegrees are 192, as given by Kimberley’s formula, except the degree is 96 when q = 6.
Questions:
I Is there a generalization of Kimberley’s formula for the much larger classx = p/s, y = q/s, for 1 ≤ p, q < s/2 ?
I How can we expand these computations? New algorithms? Parallel computing?
I What about φ3 (x, y , z) and n-dimensional cases?
There is still much to be learned.
1. D. H. Bailey, “Finding large Poisson polynomials using four-level variable precision,” Correctness 2021:
5th International Workshop on Software Correctness for HPC Applications, 19 Nov 2021,
https://www.davidhbailey.com/dhbpapers/correct-2021.pdf.

37 / 38

Thanks!
Questions?
This talk is available at http://www.davidhbailey.com/dhbtalks/dhb-eulerc-2022.pdf

38 / 38