New computations of Poisson polynomials:
An application of very high-precision floating-point arithmetic
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Lab (retired)

July 3, 2023

1 / 18

Polynomials arising from the Poisson potential function
The Poisson potential function appears in numerous applied math contexts, ranging
from mathematical physics to sharpening iPhone images. A simple 2-D instance is:
1 X cos(mπx ) cos(nπy )φ2 (x , y ) = 2
π m,n odd
m2 + n2
A 2013 study numerically discovered and then proved the intriguing fact that when x
and y are rational, then φ2 (x , y ) satisfies
1
φ2 (x , y ) = log β(x , y )πwhere β(x , y ) is algebraic, i.e., the root of an integer polynomial of some degree m.
By computing high-precision numerical values of φ2 (x , y ) for various specific rational x
and y , and applying variants of the PSLQ program, we were able to produce theexplicit minimal polynomials for α in several simple specific cases.
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.

2 / 18

Key breakthrough: Borwein’s fast algorithm to compute φ2 (x , y )The original formula for φ2 (x , y ) converges much too slowly for numerical evaluation.
But this formula, found by Jonathan Borwein (deceased 2016), is remarkably efficient:
φ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 ). Note that these series converge very rapidly:θ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
3 / 18

The PSLQ integer relation algorithm
Let X = (xk ) be an (m + 1)-long real or complex vector. An integer relation algorithm
finds a nontrivial integer vector A = (ak ) such thata0 x0 + a1 x1 + · · · + am xm = 0.
I The PSLQ algorithm and the multipair PSLQ algorithm are commonly used for
integer relation detection. Variants of the LLL algorithm can also be used.
I Integer relation detection requires very high precision floating-point arithmetic:
at least (m + 1) · maxk log10 |ak | digits (typically hundreds or thousands of digits),
both in the input data and the algorithm, or else the true relation will be lost in asea of numerical artifacts.
I Fast variations of PSLQ and multipair PSLQ utilize two, three or even four levelsof precision, doing as much computation as possible with only double precision,and switching to higher levels only when needed.

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.

4 / 18

Numerically finding minimal polynomials using PSLQ
Integer relation algorithms can be used to recognize a computed numerical value as theroot of an integer polynomial of degree m.
Example: α is suspected to be an algebraic number of degree 8 or less:
α = 2.1195912698291751313298483349346871106280 . . .
What is its minimal polynomial?
Method: Compute the vector (1, α, α2 , · · · , α8 ), then apply a variant of the PSLQ
algorithm. This produces the integer vector (1, −216, 860, −744, 454, −744,
860, −216, 1), so that α appears to satisfy the polynomial:
0 = 1 − 216α + 860α2 − 744α3 + 454α4 − 744α5 + 860α6 − 216α7 + α8

5 / 18

High-level computational algorithm
1. Given rationals x = p/s and y = q/s, select a conjectured minimal polynomial
degree m, precision levels P1 and P2 and other parameters for the run.
2. Calculate φ2 (x , y ) to P2 -digit precision using Borwein’s fast algorithm. Whendone, calculate α = exp(8πφ2 (x , y )) and generate the (m + 1)-long vectorx = (1, α, α2 , · · · , αm ), to P2 -digit precision. Note: the polynomials are simpler interms of α = (β(x , y ))8 = exp(8πφ2 (x , y )).
3. Apply a variant of the PSLQ algorithm to find an integer relation for x .
4. If no numerically significant relation is found, try again with a larger degree m or
higher precision. If a relation is found, employ Mathematica or Maple to ensurethat the polynomial is irreducible. Alternatively, rerun the problem with thedegree m reduced by one, to ensure that no numerically significant relation isfound with this smaller degree.

6 / 18

2013: Some initial Poisson polynomial computations
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 + α8

Questions:
I Given s, what is the degree of the corresponding minimal polynomial?
I Note that when s is even, the polynomial is palindromic, i.e., coefficientsak = am−k . Does this pattern hold for all even s?
These computations were very expensive and required very high precision.
Help! More powerful computational tools are required.

7 / 18

Kimberley’s formula for the degree of the polynomial
Based on these 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 rule: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

Questions:
I Does Kimberley’s formula hold for larger s?
I Does the palindromic property hold for larger even s?
8 / 18

2016: Improvements to the Poisson polynomial program
1. A new thread-safe arbitrary precision floating-point package.
2. Multiprecision software may optionally utilize the MPFR / GMP packages foreven faster performance.
3. A new 3-level multipair PSLQ program, up to 5X faster:a Double precision (appr. 15 digit accuracy).
b Medium multiprecision: varies from 100 digits to 1200 digits on Poisson problems.
c Full multiprecision: varies from 2,500 digits to 50,000 digits on Poisson problems.

4. Faster hardware: an 8-core MacPro system with Intel processor.

I D. H. Bailey, J. M. Borwein, J. Kimberley and W. Ladd, “Computer discovery and analysis oflarge Poisson polynomials,” Experimental Mathematics, 27 Aug 2016, vol. 26, 349-363,
https://www.davidhbailey.com/dhbpapers/poisson-res.pdf.
9 / 18

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+6700 α99

10 / 18

Palindromic polynomials
From our results, in the case x = 1/s, y = 1/s where s is even, the resulting
polynomial is always palindromic (ak = am−k ). For instance, when s = 16,
p16 (α) = 1 − 1376α1 − 12560α2 − 3550496α3 + 81241720α4 − 169589984α5+ 1334964944α6 − 24307725984α7 + 238934926108α8 − 1043027124704α9+ 2328675366384α10 − 3219896325280α11 + 4238551472456α12− 10247414430048α13 + 28552105805904α14 − 55832851687968α15+ 70020268309062α16− 55832851687968α17 + 28552105805904α18 − 10247414430048α19+ 4238551472456α20 − 3219896325280α21 + 2328675366384α22− 1043027124704α23 + 238934926108α24 − 24307725984α25 + 1334964944α26− 169589984α27 + 81241720α28 − 3550496α29 − 12560α30 − 1376α31 + α32

Nitya Mani, a student at Stanford University, pointed out that if α is a root of apalindromic polynomial such as this, then α + 1/α is a root of a transformed
polynomial of half the degree. This fact has been used to greatly accelerate thecomputation of Poisson polynomials when s is an even integer.

11 / 18

2016: A proof of Kimberley formula for the case x = 1/s, y = 1/s
Some observations from the polynomials produced by the program:
I The algebraic number αs is the largest real root of the associated polynomial.
I The polynomial has ϕ(s) real roots, where ϕ is the Euler totient function.
I Roots are connected to a sequence of polynomials defined in a 2010 paper by
Savin and Quarfoot of the University of Utah, which was found by doing an
Google search for “387221579866,” a coefficient of the polynomial for the case
(1/11, 1/11).
These observations ultimately led to a proof, by Watson Ladd of U.C. Berkeley, ofKimberley’s formula and also the palindromic property, in the specific casex = y = 1/s.
I D. H. Bailey, J. M. Borwein, J. Kimberley and W. Ladd, “Computer discovery and analysis oflarge Poisson polynomials,” Experimental Mathematics, 27 Aug 2016, vol. 26, 349–363,
https://www.davidhbailey.com/dhbpapers/poisson-res.pdf.
12 / 18

What about the many cases with x 6= y ?
All of the computations and results mentioned above are for the case x = y = 1/s forsome positive integer s.
What about rationals x = p/s, y = q/s, with p 6= q? Initial results indicated that
Kimberley’s formula does not hold for these more general rationals.
Is there a generalization of Kimberley’s formula that holds in these other cases? Does
the palindromic property hold for these other cases?
To address this more general problem requires many times more computation than forthe x = y = 1/s cases.
A familiar refrain: Help! More powerful computational tools are required.

13 / 18

2022: New computational tools for Poisson polynomials
1. A new high-level multiprecision package, based on integer arithmetic, that is
approximately 3X faster than before, nearly as fast as using MPFR / GMP.
2. A new 4-level multipair PSLQ program:a Double precision (approx. 15 digit accuracy).
b Quad precision (approx. 32 digit accuracy).
c Medium multiprecision: varies from 100 digits to 1200 digits on Poisson problems.
d Full multiprecision: varies from 2500 digits to 100,000 digits on Poisson problems.

3. The code does as much computation as possible in double precision; automaticallyshifts to higher precision when required (very challenging to reliably program).
4. Faster hardware: A 10-core MacStudio, with Apple’s low-power M1 Pro processor,
thus mercifully helping the author avoid bankruptcy from electric bills.

14 / 18

2023: New computer runs
I The new software has been used to compute the minimal polynomials for theentire set of cases (p/s, q/s), where 1 ≤ p ≤ q < s/2 and 10 ≤ s ≤ 36, and also
for s = 38, 40, 42 and s = 50 (a total of 2206 cases).
I These runs required up to 32,000-digit floating-point arithmetic and thousands ofprocessor-core hours run time.
Results: The following modification of Kimberley’s formula appears to hold in all cases(1 ≤ p ≤ q < s/2 and gcd(p, q, s) = 1):
1. For the cases x = y = p/s, Kimberley’s formula holds; further, for fixed s, allthese cases share the same minimal polynomial).
2. For the cases x = p/s, y = q/s with s odd, Kimberley’s formula holds (except fora few where the correct degree is half Kimberley’s rule).
3. For the cases x = p/s, y = q/s, with s even and both p and q odd, Kimberley’s
formula holds (except for a few where the correct degree is half Kimberley’s rule).
4. For the cases x = p/s, y = q/s, with s even and one of p or q is even, the correct
degree is twice Kimberley’s formula (except for a few where the correct degree isequal to Kimberley’s rule).
15 / 18

Sharing of minimal polynomials
One intriguing finding from the latest computations is that many cases for a given s
share the same minimal polynomial, even though the α numerical values are different.
For example, when s = 17, the (x , y ) cases(1/17, 1/17), (2/17, 2/17), (3/17, 3/17), (5/17, 5/17), (6/17, 6/17), (7/17, 7/17), (8/17, 8/17)

all satisfy the same degree-64 minimal polynomial:
1

2

3

4

5

6

7

8

+ 1 + 6912α − 1023008α + 535196800α + 7742027760α − 2451239864832α + 140264665723552α − 2494265652888704α + 18453445522215032α
910+ 21614293158955264α − 1840469978381611680α+ 21048406812137688311168α

15− 123756069205191278016740α− 14862548097474240887146768α+ 11985439092809681992002048α− 1789602095389051533149533952α21

28− 2451983939727870545406928048α− 830074840813669608610951352α− 6019517253583536219231909888α+ 4743446591055878746050587136α

62− 670240α+ 3292609205079608858656α+ 3891480748650616α

63

+ 1α

− 1408α41+ 269366792757186303037874944α

55− 89825284727190400α

29

56

50

45

18

37

42

− 96596567511508184274883040α
46

51

57

58

23

31

− 13334548483907481046238812288α

3835− 3898119322387426442055756416α+ 2101492937309911776817793664α

3943+ 46311532722057913438161792α47+ 331427117746835861477504α+ 76552613117134517712α− 3704022727520α

14

19+ 7693234584556635821267584α

2730− 3881818694457698660972764704α− 4563932639248948435424α− 4156271487999506323835036544α

34− 738833647673944491136α− 154854254425344α26+ 11966489230110362129440701856α− 2079969173966458011379616α

13− 259625664937972467300807296α+ 8780067564346216307831284640α
33+ 1143759465759937297408α− 2671051287612630032421280α+ 3055552833334608777606289376α

36+ 8153783303351403692882944α

17

22

25

− 17345855629600599241800189696α

4044− 22155672572673873192657168α

12− 219265475764921569840α+ 44351668349396581870408736α

32+ 17362857489419448630866293318α11+ 662708878348907477250816α

24+ 4903963676671959157531751248α

49

16

20

+ 803186115899676703948238664α− 7235322856083561662208α+ 26560170568288794240α59+ 404224147840α

5248− 18856552838875733014756α
53− 1424154241008650752α
60− 125943824α+ 342676113911934816α5461+ 62013440α64

16 / 18

July 2023: New results: The Poisson ψ function
The 2013 study mentioned the closely related function
ψ2 (x , y ) =1 X cos(πmx ) cos(πny ).
π 2 m,n evenm2 + n2

As with φ2 (x , y ), the authors found that when x and y are rational, thenψ2 (x , y ) = 1/π · log(β(x , y )), for algebraic β(x , y ).
As the present study was being concluded, DHB discovered formulas andcomputational techniques to find minimal polynomials for ψ2 (x , y ). Intriguing initial
results have been obtained.
However, the computations and analysis here are significantly more challenging thanwith φ2 (x , y ), requiring, among other things, much higher precision — up to
100,000-digit or more floating-point arithmetic even for modest-sized s.
Details will be provided in a separate report.
17 / 18

Conclusions
I New software has been used to compute the minimal polynomials for α
corresponding to φ2 (x , y ), for a total of 2206 cases. These runs required up to
32,000-digit floating-point and several CPU-months run time.
I A modification of Kimberley’s formula appears to hold in all cases, although noproof is known. For a given s, many cases share the same minimal polynomial.
I Exploring 3-D Poisson polynomials and polynomials associated with ψ2 (x , y ) will
require much more computation, with even higher levels of numeric precision(100,000 or more digits).
I A familiar refrain: Help! More powerful computational tools are required.
This talk is available athttps://www.davidhbailey.com/dhbtalks/dhb-fpt-2023.pdf.
Full results are available athttps://www.davidhbailey.com/dhbpapers/poisson-2023.pdf.
18 / 18