Finding large Poisson polynomials using four-level variable
precision
David H. Bailey
Lawrence Berkeley National Laboratory (retired)
University of California, Davis, Department of Computer Science
This talk is available at:
http://www.davidhbailey.com/dhbtalks/dhb-sc21-correct.pdf

1 / 19ICERM report on reproducibility in mathematical and scientific computing
The field of high-performance computing has lagged efforts in other fields to ensure
reproducibility. One key issue is numerical reproducibility:
Numerical reproducibility has emerged as a particularly important issue, since
the scale of computations has greatly increased in recent years, particularly
with computations performed on many thousands of processors and involving
similarly large datasets. Large computations often greatly magnify the levelof numeric error, so that numerical difficulties that were once of little import
now are large enough to alter the course of the computation or to draw into
question the overall validity of the results.

I D. H. Bailey, J. M. Borwein, O. Caprotti, U. Martin, B. Salvy and M. Taufer, “Opportunities andchallenges in 21st century mathematical computation: ICERM workshop report,” 31 Jul 2014,
https://www.davidhbailey.com/dhbpapers/ICERM-2014.pdf.
2 / 19

Commonly used formats for floating-point computing
NameIEEE half
ARM half
bfloat16
IEEE single
IEEE double
IEEE extended
IEEE quad
double-double
quad-double
double-quad
arbitrary

Sign
1
1
1
1
1
1
1
1
1
1
variesNumber of bits
Exp. Mant. Hidden
5
10
1
5
10
1
8
7
1
8
23
1
11
52
1
15
64
0
15
112
1
11
104
2
11
208
4
15
224
1
varies varies varies

Digits
3
3
2
7
15
19
34
31
62
68
varies3 / 19Advantages of variable precision
Compared with using a high fixed level of precision for the entire application,
employing variable precision usually results in:
I Faster processing.
I Better cache utilization.
I Lower run-time memory usage.
I Lower offline data storage.
I Lower energy costs.

Challenge:
I Unless variable precision is carefully managed, there is significantly greater

potential for loss of numerical accuracy and reproducibility.

4 / 19Finding Poisson polynomials: Application of variable precision computing
Lattice sums related to the Poisson potential function, which naturally arise in studies
of gravitational and electrostatic potentials, have been studied for many years in themathematical physics community, and recently have been applied to image sharpening:
φ2 (x , y ) =1 X cos(mπx ) cos(nπy )π 2 m,n odd
m2 + n2

In 2012 researchers discovered that when x and y are rational numbers,
1
φ2 (x , y ) = log α(x , y ),πwhere α(x , y ) is an algebraic number, namely the root of a degree-m integer
polynomial. What is the connection between x and y and the degree of α?
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.
I 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.
5 / 19Finding Poisson polynomials, cont.
By computing high-precision numerical values of φ2 (x , y ) for various specific rational x
and y , using a three-level variable precision code (e.g., double/500-digit/10,000-digit),
and applying a variant of the PSLQ algorithm, we found the following results:
s
5
6
7
8
9

10

Minimal polynomial of α 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

6 / 19Kimberley’s formula
Based on these preliminary results, Jason Kimberley conjectured that the degree m(s)of the minimal polynomial associated with the case x = y = 1/s is given by this rule:Set m(2) = 1/2. Otherwise for primes p congruent to 1 mod 4, setm(p) = int2 (p/2), where int denotes greatest integer, and for primes p congruent to3 mod 4, set m(p) = int (p/2)(int (p/2) + 1). Then for any other positive integer s
whose prime factorization is s = p1e1 p2e2 · · · prer ,r −1m(s) = 4rY

2(ei −1)pim(pi ).

i=1

Additional three-level precision computations, (e.g., double/2000-digit/64,000-digit),
confirmed Kimberley’s formula for all s up to 40, and for even s up to 52.
By examining the computed results, connections were found to a sequence ofpolynomials defined in a 2010 paper by Savin and Quarfoot. These observations
ultimately led to a proof of Kimberley’s formula and several other facts.
What about the more general case x = p/s and y = q/s, for integers p, q, s?

7 / 19Thread safety in arbitrary precision libraries
Most existing arbitrary precision libraries today are not thread-safe. Thus applications
written to use them do not work in a multi-thread shared-memory parallel environment.
Currently the GNU MPFR package is the only arbitrary precision library that isthread-safe. It is also very fast (esp. for transcendentals), and correctly rounds results.
However:
I MPFR is only a low-level library; it does not include high-level language bindings
required to support full-scale application codes.
I MPFR requires a lengthy installation procedure with administrator permission (toinstall MPFR and GMP on the user’s computer), which discourages many users.
I If a thread-safe build option is not properly invoked during installation, theresulting software is not thread-safe.
I L. Fousse, G. Hanrot, V. Lefevre, P. Pelissier and P. Zimmermann, “MPFR: A multiple precision binary
floating-point library with correct rounding,” ACM Transactions on Mathematical Software, vol. 33(2007), https://doi.org/10.1145/1236463.1236468.
8 / 19

DHB’s new MPFUN2020 arbitrary precision library
I An all-Fortran design that can be compiled in a few seconds with any Fortran-2008

compliant compiler, including the gfortran, Intel and NAG Fortran compilers, on a variety
of systems, including Mac OS X and various Linux platforms.
I A 100% thread-safe design.
I A full-featured high-level Fortran language interface, so that only type statements and a

few other minor changes are required to convert most double precision codes.
I Support for common transcendental functions (sin, cos, exp, etc.) and special functions.
I Support for arbitrary precision real and complex datatypes, plus double and quad.
I Support for both a full and a medium precision datatype, which in some applications(such as Poisson problems) results in lower memory usage and memory traffic.
I FFT-based multiplication for faster performance at very high precision.
I A comprehensive test suite, together with six full-scale application test codes.

A separate version calls MPFR for all lower-level operations. It is roughly 15% fasteron most applications (several times faster for transcendental functions).
9 / 19

DHB’s new four-level multipair PSLQ code for finding integer relations
Given an (m + 1)-long input vector X = (xi ) of high-precision floating-point values, aninteger relation algorithm finds an (m + 1)-long vector of integers (ai ) such thata0 x0 + a1 x1 + a2 x2 + · · · + am xm = 0,
to within the tolerance of the numeric precision being used. As a simple illustration,
suppose one suspects that the real constant
α = 2.11959126982917513132984833493468711062807145783249543921443000325120624764609847 . . .

is an algebraic number of degree eight. After computing the vector (1, α, α2 , · · · , α8 )and applying the multipair PSLQ integer relation algorithm, the relation(1, −216, 860, −744, 454, −744, 860, −216, 1) is produced, so that α satisfies
1 − 216α + 860α2 − 744α3 + 454α4 − 744α5 + 860α6 − 216α7 + α8 = 0.
DHB’s new multipair PSLQ code employs four levels of precision: (a) IEEE double; (b)IEEE quad; (c) medium precision, typically 100-1000 digits; and (d) full precision,
typically 5,000 to 50,000 digits.
10 / 19

How can one tell if a computed polynomial is real and not an artifact?
Iterations of the multipair PSLQ algorithm develop a sequence of invertible integer
matrices An and inverses Bn so that y = Bn · X has steadily smaller entries, until one
entry of y is smaller than the epsilon specified for detection.
The size of the drop in min(|yi |) when the relation is detected can be viewed as aconfidence level that the discovered relation is a real mathematical relation.
0$

1$ 11$ 21$ 31$ 41$ 51$ 61$ 71$ 81$ 91$ 101$ 111$ 121$ 131$ 141$ 151$ 161$ 171$ 181$ 191$

!50$

Log10&[min(yi)]&

!100$

!150$

!200$

!250$

11 / 19

High-level computational algorithm
1. Given rationals x = p/s and y = q/s (typically satisfying 1 ≤ p, q < s/2 ≤ 50), select a
conjectured minimal polynomial degree m, a medium precision level P1 digits, a full
precision level P2 digits and other parameters for the run.
2. Calculate φ2 (x , y ) to P2 -digit precision using the following formula:φ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 rapidlyP∞2
convergent formulas. Example: θ1 (z, q) = 2 k=1 (−1)k−1 q (2k−1) /4 sin((2k − 1)z).
3. Calculate α = exp(8πφ2 (x , y )) and generate the (m + 1)-long vector
X = (1, α, α2 , · · · , αm ) to P2 -digit precision.
4. Apply the four-level multipair PSLQ algorithm to find an integer relation for X , if one
exists, to the precision being used (P2 digits).
5. If a numerically significant relation is not found, try again with a larger degree m or ahigher full precision level P2 . If a relation is found, employ the polynomial factorization
facilities in Mathematica or Maple to ensure that the resulting polynomial is irreducible.
12 / 19Performance: New vs old codes
Arbitrary
precision package
Old all-FortranNew all-Fortran
New MPFR-based

Multipair
PSLQ
3-level3-level4-level3-level4-levelRun time (CPU seconds)
x = y = 1/29 x = y = 1/35
247819.68
220202.68
55691.66
49221.37
52920.6747093.43
45038.87
41442.57
45411.23
39762.23

The two application cases displayed are:I x = y = 1/29: requires 19,700-digit arithmetic (full precision); produces apolynomial of degree 196 with coefficients ranging from 1 up to roughly 1087 .
I x = y = 1/35: requires 18,900-digit arithmetic (full precision), producing apolynomial of degree 192 with coefficients ranging from 1 up to roughly 1085 .
13 / 19

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
14 / 19

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
15 / 19

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
16 / 19

Research findings
Note that many of these results are not consistent with Kimberley’s formula, which has
been proven only for the cases p = q = 1 (i.e., x = 1/s; y = 1/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, when the degrees are only half
as large. For instance when s = 36, degrees alternate between 144 and 288, yet whenq = 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, 1 ≤ p, q < s/2?
I How can we greatly expand the scope of these computations? A new integer
relation algorithm?
I How can massively parallel computing be employed? Parallel implementations todate have yielded only modest speedups, compared with the best serial code.

17 / 19Degree-100 minimal polynomial found for the case x = y = 1/25
−1−33300 α1−8029750 α2+1542379500 α3−293348778825 α4−28080295402080 α5−6629781350330800 α6+58722280292148000 α7−7197243181862115100 α8+182328184490475806640 α9+1047464861734017883208 α10+1806185067754657491760 α11−2794588132471771782562700 α12+43132729144820226386860000 α13−256411580898396643830482000 α14+1396700356872030954211266400 α15+9088870493426588389401095750 α16−733599064553317639364961825000 α17+13704149574920386191809571588500 α18−157819237362090093545769457641320 α19+1321853578066306850140386747623398 α20−8347129396578224228862993897237280 α21+38212880291083166937296196441746800 α22−104022238521780343721711917072676000 α23−74032065945345957753107951351295900 α24+2695579175102470750839898367021130608 α25−17580580889577528003991569683501290200 α26+74952710236773360685253102982788566000 α27−239080432366992082679933177546734267500 α28+594175511524312220492120757330324484000 α29−1200481201725892359379030176076009068400 α30+2370791216795342465482695956290874196000 α31−6958365090143220739303965994873839612375 α32+29340961737653182271816694472006500014100 α33−118465202788489603844178657796002545834090 α34+394388553221080980775811796973709891661972 α35−1059773119301328191261276484738640092975935 α36+2274309666161593753722296548558864618267200 α37−3729748087865447793480560952748424761349600 α38+3877725636887659946158475367090512632719680 α39+846723867393661736729146827990176999234824 α40−15960428173973401824575411855716955373304480 α41+47118520726980412118527981001595834084077200 α42−96560391964393397280402148487438432608004000 α43+159442357501276027171661799158377030358819240 α44−222707090658310894211715578612740014870494016 α45+268336339651868593336977017479411555696276960 α46−280387518545606966350071856817231380200513600 α47+252610498127566038453865516310244367924136500 α48−192275255815619622137470846245628244410759600 α49+117222736648159046124836188054855685036349976 α50−47209171248166363599083489199392291307559600 α51−4942922820170238254884354265559936285783500 α52+36400498321158716446352704351511072445758400 α53−50714975271028177138621429774290438585502240 α54+52508039392118222466312798018961036691290304 α55−45426960236840080028112470330711413632425560 α56+33189561099972574095810629259221045183420000 α57−20184727758357355309386195853692979100127600 α58+10078674543222169818704320601337376637695840 α59−4216672803412456875435463721793522689401592 α60+1722566374653535529797841998864588575882560 α61−940122938034032272236892859372976755103200 α62+672500761678739947489218957049737685556800 α63−461584164001204318877039495371276336862015 α64+265582157437208924810452188771019969018068 α65−127220613278091048562581392274127224367210 α66+52420552698515377367407956646216340562900 α67−19432168360605247012293398573854487548375 α68+6715976041225158674019447203640258687520 α69−2160843064282037086551775761827735549808 α70+617282737944009515075102146542202738080 α71−145142540389934636970079472194894728300 α72+25729074826332880513711838275787958000 α73−3264527861366792560817221584368055000 α74+473260992168563618404463109451174000 α75−178858482702935977606858126791577500 α76+66028976498454634546379541782684000 α77−15761469487294858226633454090400400 α78+2480913858831352983037735837103840 α79−283939423573690705553213818042394 α80+27924016290268330781557752871960 α81−1890264935738308757140819195500 α82−197646050954205989496871249000 α83+92058025400867206020845473350 α84−15898833972285544673912261280 α85+1776153714893883182147476400 α86−134561774900933516621588000 α87+5713502521319481241590900 α88−28013036278808426259920 α89−10740617849680182084536 α90+416466536564927187120 α91−5335454032157718300 α92−61275002107276000 α93−662049081315760 α94+123800803998624 α95−330633806665 α96−8699478100 α97+13986250 α98
99

18 / 19

Suggestions for future research in variable precision computing
1. Seamlessly incorporating an arbitrary precision facility, say based on the MPFR
package, in future Fortran and C language standards.
2. Checking all available software supporting variable and arbitrary precision forthread safety, and fixing or deprecating libraries that are not strictly thread-safe.
3. Rethinking automatic type conversions and mixed-mode operations, which hinder
the development of correct and reproducible variable precision computing, in
DHB’s opinion. At the least, could some of these features be optional?
4. Investigating novel numerical algorithms, error control and error analysis
techniques, appropriate not only for single and double precision on a single
processor system, but also for the realm of very high precision and very highly
parallel computing.
5. Developing fundamentally new paradigms for variable precision computing,perhaps along the lines of Gustafson’s “unum” framework.
This talk is available at:
http://www.davidhbailey.com/dhbtalks/dhb-sc21-correct.pdf
19 / 19