Discovery of large Poisson polynomials using a new arbitrary
precision software package
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Lab (retired)
University of California, Davis, Dept. of Computer Science

November 11, 2015

1 / 34

Some applications of high-precision arithmetic
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. Electromagnetic scattering theory (32–100 digits).
7. The Taylor algorithm for ODEs (100–600 digits).
8. Ising integrals from mathematical physics (100–1000 digits).
9. Problems in experimental mathematics (100–100,000 digits).
These and other examples are described in:
I D. H. Bailey and J. M. Borwein, “High-precision arithmetic in mathematical physics,”

Mathematics, vol. 3 (2015), pg. 337-367.
2 / 34

The PSLQ integer relation algorithm
Let X = (xk ) be an (m + 1)-long real or complex vector. An integer relation algorithm
such as PSLQ finds a nontrivial integer vector A = (ak ) such thata0 x0 + a1 x1 + · · · + am xm = 0.
I The multipair PSLQ algorithm is a more efficient and moderately parallelizablevariant of PSLQ, the most widely used integter relation algorithm (although someuse a variant of LLL).
I Integer relation detection (by any algorithm) requires very high precision: at least(m + 1) · maxk log10 |ak | digits, both in the input data and the algorithm.
I Multipair PSLQ is extremely efficient with precision — it can usually detect a
relation when the precision is only a few percent higher than this minimum level.

I D. H. Bailey and D. J. Broadhurst, “Parallel integer relation detection: Techniques andapplications,” Mathematics of Computation, vol. 70, no. 236 (Oct 2000), 1719–1736.
3 / 34

Multipair PSLQ, continued
I Iterations of the multipair PSLQ algorithm develop a sequence of invertible

integer matrices An , their inverses Bn and real matrices Hn (in lower quadratureform), so that the reduced vector Y = Bn · X has steadily smaller entries.
I A relation is detected when the size of smallest entry of the Y vector drops below
the epsilon specified for detection (typically roughly 10−p , where p is the precisionlevel 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.
I See next viewgraph for an illustration of the steadily decreasing value of min |yk |,leading up to iteration 199, when the relation was discovered.

4 / 34

Decrease of log10 (min |yi |) in a multipair PSLQ run

5 / 34

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 = 6971 − 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
6 / 34

High-precision numerical integration
Given f (x ) defined on (−1, 1), define g(t) = tanh(π/2 sinh t). By setting x = g(t),
Z ∞

Z 1

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

−∞

N
X

wj f (xj ),j=−Nwhere xj = g(jh) and wj = g 0 (jh), independent of the function f (x ), for some intervalh. The summation continues until the terms are consistently smaller than epsilon.
We have found that tanh-sinh is the best general-purpose integration scheme forfunctions with vertical derivatives or singularities at endpoints, or for any function atvery high precision (> 1000 digits). Otherwise we use Gaussian quadrature.

I D. H. Bailey, X. S. Li and K. Jeyabalan, “A comparison of three high-precision quadratureschemes,” Experimental Mathematics, vol. 14 (2005), 317–329.
7 / 34

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 ∞
4 ∞
1
dun
du1
Cn :=
···
···
P
2
n! 0
un
n
0
(u + 1/u ) u1
Dn :=

4
n!

Z ∞
0

Z 1

En

=

···

Z ∞
0

Z 1

j=1

j

j

Q



2

i<j

ui −ujui +ujPn
j=1 (uj + 1/uj )

2

dun
du1
···
u1
un

2



uk − uj 2···
dt2 dt3 · · · dtn ,u + uj0
0
1≤j<k≤n k
Ywhere in the last two lines uk = t1 t2 · · · tk .
I D. H. Bailey, J. M. Borwein and R. E. Crandall, “Integrals of the Ising class,” Journal of Physics

A: Mathematical and General, vol. 39 (2006), pg. 12271–12302.
8 / 34

Limiting value of Cn : What is this number?
The Cn integrals can be converted to 1-D integrals involving the 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.
9 / 34Other Ising integral evaluations found using multipair 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.
10 / 34

Free software for high-precision floating-point computation
High-precision is built into commerical mathematical software packages such as Maple
and Mathematica, but many researchers use one of the following:
1. ARPREC: Includes 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: Supports high-precision integer, rational and floating-point calculations.
http://gmplib.org.
4. MPFR: C library for arbitrary precision with exact rounding, based on GMP.
http://www.mpfr.org.
5. MPFR C++: A high-level C++ interface to MPFR with a thread-safe option.
http://www.holoborodko.com/pavel/mpfr.
6. Pari/GP: A computer algebra system that includes high-precision, with numerous
transcendental functions. http://pari.math.u-bordeaux.fr.
7. QD: Performs “double-double” (31 digits) and “quad-double” (62 digits) arithmetic, withtranscendental functions. High-level interfaces for C++ and Fortran-90.
http://crd.lbl.gov/˜dhbailey/mpdist.
11 / 34

Thread safety: A major challenge for high-precision software
I Because of greatly magnified run times (often 1000X or more) of high-precisioncomputations, highly parallel implementations are in order. However, mosthigh-precision packages are not thread-safe.
I Many packages employ global read/write variables, e.g., for transcendental

function evaluation, which are deadly for thread safety.
I The working precision level, a global variable that is changed frequently within the

package itself and often by users also, is particularly troublesome.
One thread-safe package: The MPFR C++ package, which calls the lower-level MPFR
package, is thread-safe (provided MPFR is compiled with the thread-safe build option).
MPFR is a very well-designed, features correct rounding and is the fastest low-levelpackage currently available.
But there is no thread-safe high-level Fortran package.

12 / 34

MPFUN2015: A new thread-safe arbitrary precision package
DHB has written a new package, 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.
The two versions are “plug-compatible” in the sense that applications written for one
version also work with the other.
A C++ version is in the works.
The software is available at http://www.davidhbailey.com/dhbsoftware.
Full details of design, algorithms, installation and usage are given in
I D. H. Bailey, “MPFUN2015: A thread-safe arbitrary precision computation package,”manuscript, 1 Oct 2015. http://www.davidhbailey.com/dhbpapers/mpfun2015.pdf.
13 / 34

Features of MPFUN-Fort and MPFUN-MPFR
I 100% thread safe (provided MPFR is compiled with the thread-safe build option).
I Both employ advanced algorithms, including FFT-based arithmetic.
I Both versions include a full-featured, high-level Fortran interface (implemented

using custom datatypes and operator overloading).
I Both provide full support for both multiprecision real (MPR) and multiprecisioncomplex (MPC) datatypes, as well as a wide variety of mixed-mode operations(e.g., double precision times multiprecision real, multiprecision complex to an
integer power, etc.).
I All common transcendentals and several special functions are included.
I Precision is scalable to over one million digits.
I Both provide an effective solution to the double precision accuracy problem that
plagues many high-precision packages.

14 / 34

Solution to the thread safety problem
I The user sets a default (maximum) working precision level in a compile-time datastatement.
I A “working precision” carried in data structure itself is associated with each

multiprecision (MP) variable and array element.
I The result of any operation or function involving MP variables or array elements

inherits the largest working precision level of the input operands.
I The built-in function mpreal allows one to convert a MP variable or arrayelement to one with a different working precision level.
I The built-in function mpwprec returns the working precision level currently

assigned to any MP variable or array element.

15 / 34

Converting existing codes to use MPFUN-Fort or MPFUN-MPFR
I Place this line in every subprogram that contains a MP variable or array:use mpmodule
I To designate a variable or array as multiprecision real (MPR) or complex (MPC):

type (mp real) a, b(10), c(10,10); type (mp complex) d, e, f(20)
I All basic and mixed-mode operations are supported, e.g.:d = a + b(i) * sqrt(3.d0 - c(i,j))
Functions provided:
I Standard Fortran functions with MPR arguments:
abs, acos, aint, anint, asin, atan, atan2, cos, cosh, dble, exp, log,
max, min, sign, sin, sqrt, tan, tanh
I Stanford Fortran functions with MPC arguments:

abs, aimag, conjg, cos, dcmplx, exp, log, sin, sqrt
I Special functions:

bessj, gamma, incgamma, mpberne, zeta, zetaem

16 / 34

A solution to the double precision accuracy loss problem
I If a user code includes a statement such as a

= b + d / 3.d0, where a and b
are multiprecision real (MPR) and d is double precision (DP), then by establishedlanguage conventions the division is only performed to DP accuracy, and this low
accuracy propagates to the result a.
I Similar difficulties afflict multiprecision programs written in other languages,
including C and C++.
I Both MPFUN-Fort and MPFUN-MPFR check every DP constant, variable andexpression value in MP statements at execution time to see if it has more than 40
significant bits. If so, it is flagged as an error. This feature catches 99.99% of the
precision loss problems.
I Special routines handle legitimate DP data with more than 40 bits.
I To convert a floating constant to full precision, enclose it in apostrophes, e.g.,
a = ’1.234’

17 / 34

Libraries and applications included with MPFUN-Fort and MPFUN-MPFR
Description
MPFUN-Fort library
MPFUN-MPFR library∗One-level standard PSLQ code
One-level multipair PSLQ code
Two-level multipair PSLQ code
Three-level multipair PSLQ code
Tanh-sinh quadrature code, with 14 integrals
Parallel tanh-sinh code, with 14 integrals
Poisson polynomial code with 3-level multipair PSLQ
Total

Lines of code
33,640
11,540
752
905
1694
2129
871
949
2217
54,697∗ To install MPFUN-MPFR, the GMP library (approximately 83,000 lines of C) and the

MPFR library (approximately 93,000 lines of C) must first be installed.

18 / 34

The Poisson potential function of mathematical physics
I Gravitation — the gravitational potential for a given mass distraibution.
I Electrostatics — the electric potential for a given charge distribution.
I Studies of cyclotomic polylogarithms and harmonic sums.
I Studies of crystal symmetry.
I Studies of the Madelung constant of mathematical physics.
I Studies of image sharpening.

I J. Borwein, L. Glasser, R. McPhedran, J. Wan and J. Zucker, Lattice Sums: Then and Now,

Encyclopedia of Mathematics and its Applications, 150, Cambridge University Press, 2013.
I S. Kanemitsu, Y. Tanigawa, H. Tsukada and M. Yoshimoto, “Crystal symmetry viewed as zeta

symmetry,” p. 91–130 in T. Aoki, S. Kanemitsu, M. Nakahara and Y. Ohno, eds., Zeta
Functions, Topology and Quantum Physics, Springer, New York, 2005 (reprinted 2010).
19 / 34

The Poisson potential function problem
In a 2013 study, we numerically discovered, and then proved the intriguing fact that for
rational (x , y ), the 2-D Poisson potential function φ2 (x , y ) 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 a some integer polynomial of degree m.
By computing high-precision numerical values of φ2 (x , y ) for various specific rational x
and y , and applying a multipair PSLQ program, we were able to produce the explicit
minimal polynomials for α in several specific cases.
This study was hampered by the huge complexity and computational costs involved.
Some runs required up to 12,000-digit precision, while other runs failed.
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), 115201.
20 / 34

Kimberley’s formula for the degree of the polynomial
Based on our preliminary results, Jason Kimberley of the University of Newcastle,
Australia observed that the degree m(s) of the minimal polynomial associated with thecase x = y = 1/s appears to be given by the following:
Set m(2) = 1/2. Otherwise for primes p congruent to 1 mod 4, set m(p) = int2 (p/2),
where int denotes greatest integer, and for primes p congruent to 3 mod 4, setm(p) = int (p/2)(int (p/2) + 1). Then for any other positive integer s whose prime
factorization is s = p1e1 p2e2 · · · prer ,?

m(s) = 4r −1

rY2(ei −1)pim(pi ).

i=1

Does Kimberley’s formula hold for larger s? Why? What is the true mathematicalconnection between the pair of rationals (x , y ) and the algebraic number α?
21 / 34

New three-level multipair PSLQ program
Employs three levels of numeric precision:
I Ordinary double precision.
I Medium precision, typically 100–2000 digits.
I Full precision, typically many thousands of digits.

When an entry of the double precision reduced vector is smaller than 10−14 , or whenan entry of one of the integer-valued double precision arrays exceeds
253 ≈ 9.007 · 1015 , the medium precision arrays are updated by matrix multiplication.
Similarly, when an entry of the medium precision reduced vector is smaller than themedium precision “epsilon,” the full-precision arrays are updated.
Substantial care must be taken to manage this three-level hierarchy, and to correctly
handle numerous atypical scenarios.

22 / 34

A fast algorithm to compute the Poisson potential function φ2 (x , y )φ2 (x , y ) =1
θ2 (z, q)θ4 (z, q)log,
2π
θ1 (z, q)θ3 (z, q)

where q = e −π and z = π2 (y + ix ). Compute the four theta functions using thefollowing rapidly convergent formulas involving complex variables:θ1 (z, q) = 2
θ2 (z, q) = 2∞
X2(−1)k−1 q (2k−1) /4 sin((2k − 1)z),

k=1
∞
X

2

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

θ3 (z, q) = 1 + 2
θ4 (z, q) = 1 + 2

∞
X
k=1
∞
X

2

q k cos(2kz),
2(−1)k q k cos(2kz).

k=1
23 / 34

Empirical formulas to estimate precision, memory, run time, etc.
m(s) is the estimated polynomial degree; D is the detection level;
P1 is the medium precision level; P2 is the full precision level;
N is the number of multipair PSLQ iterations;
T is the total one-core run time in seconds; M is the memory in Mbytes.
r −1m (s) = 4rY
2(ei −1)pim(pi )

i.e., the Kimberley formulai=1

log10 D = −0.462 m2 (s),P1 = 5 m(s),P2 = 5 m(s) + 0.462 m2 (s),log10 N = 105.278/m2 (s) − 33.1073/m(s) + 4.01963 + 0.00696322 m(s),log10 T = 374.411/m2 (s) − 79.2388/m(s) + 1.93981 + 0.0176301 m(s),M = 1.6 · 10−7 · m2 (s)(5m(s) + 0.462 m2 (s)).

24 / 34

High-level computational algorithm
1. Given rationals x = p/s and y = q/s, select a conjectured minimal polynomial
degree m(s) and other parameters for the run, using the formulas above.
2. Calculate φ2 (x , y ) to P2 -digit precision using the formulas from two viewgraphs
above. When done, calculate α = exp(8πφ2 (x , y )) and generate the (m + 1)-longvector X = (1, α, α2 , · · · , αm ), to P2 -digit precision.
3. Apply the three-level multipair PSLQ algorithm to X . For larger problems, employ
a parallel version of the three-level multipair PSLQ code, using the OpenMP
construct DO PARALLEL to perform certain time-intensive loops in parallel.
4. If no numerically significant relation is found, try again with a larger degree m or
higher precision P2 . If a relation is found, employ the polynomial factorization
facilities in Mathematica and Maple to ensure that the polynomial is irreducible.

25 / 34

Application program and libraries for the Poisson calculations
Description
Poisson polynomial program∗
MPFUN-MPFR package
MPFR package
GMP package
Total

Language
FortranFortranCC

Lines of code
2,000
12,00093,00083,000190,000∗ This includes the computation of φ (x , y ) and the 3-level multipair PSLQ program.
2

26 / 34

Samples of minimal polynomials found by multipair PSLQ
s
5
6
7
8
9

10

Minimal polynomial 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 + α8These computations are very expensive. The case x = y = 1/37 required 51,000-digit
precision and produced a polynomial of degree 324, with up to 145-digit coefficients;
the computer run consumed 7.7 million core-seconds.

27 / 34

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.

28 / 34

Timings for the case x = y = 1/35
Multiprecision software
ARPRECMPFUN-MPFR
MPFUN-MPFR

PSLQ code
2-level2-level3-level

Cores
1
1
1
2
4
8
16

Run time
1.599 · 106
5.249 · 105
1.240 · 105
7.585 · 104
4.121 · 104
2.476 · 104
1.021 · 104

Speedup
1.00
3.05
12.90
21.08
38.80
64.58
156.61

The run times are wall-clock run times (in seconds), measured on a 16-core 2.4 GHz
MacPro, in a typically busy environment with similar jobs running on other cores.

29 / 34

Parallel speedups
I While a 12X parallel speedup, and a 156X overall speedup, are certainly welcome,it does not appear possible to run this code efficiently on hundreds or thousandsof processor cores.
I Achieving high speedups is fairly straightforward in the medium precision and full
precision routimes, but is not easy for the double precision (DP) routines. Thus,
by Amdahl’s law, the DP portions limit overall parallel scalability.
I It would be easy to exhibit large parallel speedups by starting, say, with the
one-level multipair PSLQ code, which performs all iterations using full precision,
but this would violate the principle that parallel implementations and speedup
figures should be based on the most efficient single-threaded algorithm available;
otherwise speedup results can be highly misleading.

I D. H. Bailey, “Misleading performance reporting in the supercomputing field,” Scientific

Programming, vol. 1., no. 2 (Winter 1992), 141–151.
30 / 34

Selected runs (degrees, precision, timings, etc.) for x = y = 1/s
s
20
24
30
32
34
36
37
38
39
40
42
44
45
46
48
50

m
32
64
64
128
128
144
324
180
288
128
192
240
288
264
256
200

log10 (D)-463.84-1883.78-1868.01-7577.07-7574.93-9570.86-48431.32-14951.64-38330.14-7580.00-16993.99-26604.14-38315.08-32036.34-30248.55-18421.18

P1
160
320
350
650
650
750
1650
900
1450
650
1000
1200
1450
1350
1350
1000

P2
700
2200
2300
8200
8200
10300
51000
16000
40000
8200
18000
28000
40000
34000
32000
20000

N
1967
9297
9064
45893
45914
62282
93125412098466715345655
150364323762660001476902415316168947

M
0.81
11.33
11.33
168.20168.20267.106579.66642.984124.24
168.20829.41
2003.33
4124.24
2921.57
2586.39974.44

C
1
1
1
1
1
1
16
1
16
1
8
8
16
16
16
8

T (sec.)2.19 · 100
7.73 · 101
1.02 · 102
5.13 · 103
5.16 · 103
9.54 · 103
4.84 · 105
3.88 · 104
2.68 · 105
5.02 · 103
1.57 · 104
7.43 · 104
2.09 · 105
1.06 · 105
8.98 · 104
2.12 · 104

C · T (sec.)
2.19 · 100
7.73 · 101
1.02 · 102
5.13 · 103
5.16 · 103
9.54 · 103
7.74 · 106
3.88 · 104
4.29 · 106
5.02 · 103
1.26 · 105
5.94 · 105
3.35 · 106
1.70 · 106
1.44 · 106
1.69 · 105

s = denominator; m = degree; D = detection level; P1 = medium precision; P2 = full precision; N =number of iterations; M = Mbytes; C = cores; T = wall clock time; C · T = total core-seconds.
31 / 34

Selected runs (degrees, precision, timings, etc.) for x = 1/s, y = q/s
s
20
24
30
32
34
36
37
38
39
40
42
44
45
46
48
50

q
3
5
7
3
3
5
2
3
2
3
5
3
2
3
5
3

m
32
32
64
128
128
144
324
180
288
128
192
240
288
264
256
200

log10 (D)-359.23-336.46-1291.21-5578.06-5411.70-6831.83-35412.09-11011.83-27943.70-5674.61-12183.50-19581.93-27857.00-23318.37-21480.15-13409.44

P1
160
320
350
650
650
750
1650
900
1450
650
1000
1200
1450
1350
1350
1000

P2
700
2200
2300
8200
8200
10300
51000
16000
40000
8200
18000
28000
40000
34000
32000
20000

N
1637
1613
6867
33829
32842
45559
69127789722
45823834273
106770232713482959346987292974122468

M
0.81
11.33
11.33
168.20168.20267.106579.66642.984124.24
168.20829.41
2003.33
4124.24
2921.57
2586.39974.44

C
1
1
1
1
1
1
16
1
16
1
8
8
16
16
16
8

T (sec.)3.06 · 100
2.83 · 100
1.01 · 102
4.16 · 103
2.53 · 103
6.93 · 103
3.42 · 105
7.42 · 103
1.84 · 105
4.13 · 103
1.16 · 104
6.18 · 104
1.47 · 105
7.28 · 104
5.93 · 104
1.63 · 104

C · T (sec.)
3.06 · 100
2.83 · 100
1.01 · 102
4.16 · 103
2.53 · 103
6.93 · 103
5.47 · 106
7.42 · 103
2.94 · 106
4.13 · 103
9.27 · 104
4.95 · 105
2.35 · 106
1.17 · 106
9.50 · 105
1.30 · 105

s = denominator; m = degree; D = detection level; P1 = medium precision; P2 = full precision; N =number of iterations; M = Mbytes; C = cores; T = wall clock time; C · T = total core-seconds.
32 / 34

Analysis and conclusions
I As far as we are aware, these computations, which employed up to 51,000-digit

precision, producing polynomials with degrees up to 324 and integer coefficients
up to 10145 , constitute the largest successful integer relation computations
performed to date.
I Kimberley’s formula was affirmed in every case x = y = 1/s, for s up to 50(except for s = 41, 43, 47, 49, which are still too costly to test).
I Research is need to understand why Kimberley’s formula holds.
I Additional research is needed to understand many other combinations, e.g.,
x = p/s and y = q/s, for different values of p, q and s, which will require evenmore extreme computation.
I While a 12X parallel speedup is certainly welcome, a scheme to efficiently employ
hundreds or thousands of processor cores is needed. A fundamentally new integer
relation algorithm may well be required.

33 / 34

Thanks!
I This talk is available at

http://www.davidhbailey.com/dhbtalks/dhb-mpfun-pphi.pdf.
I A technical report with more details, results and analysis is available at

http://www.davidhbailey.com/dhbpapers/poisson-res.pdf.
I The MPFUN-Fort and MPFUN-MPFR software packages are available at

http://www.davidhbailey.com/dhbsoftware.
I A technical report giving full details of MPFUN-Fort and MPFUN-MPFR isavailable at http://www.davidhbailey.com/dhbpapers/mpfun2015.pdf.

34 / 34