Experimental mathematics and the omega function
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Laboratory (retired)
University of California, Davis, Department of Computer Science

1 / 38

The PSLQ integer relation algorithm
Let (xn ) be a given vector of real numbers. An integer relation algorithm either finds
integers (an ) such thata1 x1 + a2 x2 + · · · + an xn = 0(to within the “epsilon” of the arithmetic being used), or else finds bounds within
which no relation can exist.
The “PSLQ” algorithm of mathematician-sculptor Helaman Ferguson is the mostwidely used integer relation algorithm.
Integer relation detection requires very high precision (at least n × d digits, where d isthe size in digits of the largest ak ), both in the input data and the algorithm itself.
I H. R. P. Ferguson, D. H. Bailey and S. Arno, “Analysis of PSLQ, An Integer Relation Finding

Algorithm,” Mathematics of Computation, vol. 68, no. 225 (Jan 1999), pg. 351–369.
I D. H. Bailey and D. J. Broadhurst, “Parallel Integer Relation Detection: Techniques andApplications,” Mathematics of Computation, vol. 70, no. 236 (Oct 2000), pg. 1719–1736.
2 / 38Ferguson’s “Umbilic Torus” sculpture at Stony Brook University

3 / 38PSLQ, continued
I PSLQ constructs a sequence of integer-valued matrices Bn that reduce the vectory = x · Bn , until either the relation is found (as one of the columns of matrix Bn ),or else precision is exhausted.
I A relation is detected when the size of smallest entry of the y vector suddenly
drops to roughly “epsilon” (i.e. 10−p , where p is the number of digits ofprecision).
I The size of this drop can be viewed as a “confidence level” that the relation is nota numerical artifact: a drop of 20+ orders of magnitude almost always indicates areal relation.
Efficient variants of PSLQ:
I 2-level and 3-level PSLQ perform almost all iterations with only double precision,
updating full-precision arrays as needed. They are hundreds of times faster thanthe original PSLQ.
I Multi-pair PSLQ dramatically reduces the number of iterations required. It was
designed for parallel systems, but runs faster even on 1 CPU.
4 / 38

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

5 / 38

The binary algorithm for exponentiation
Problem: What is 317 mod 10?
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.
Note that in Algorithm C, we never have to deal with integers larger than 9 × 9 = 81.
In general, if reducing mod n, we never have to deal with integers larger than (n − 1)2 .

6 / 38

1996 result: How to compute arbitrary binary digits of log 2
Consider this well-known formula for log 2:
log 2 =

∞
X
1
n=1

n2n

= 0.101100010111001000010111111101111101000111001111011 . . .2

Note that the binary digits of log 2 beginning after position d can be written as{2d log 2}, where {·} denotes fractional part. Thus we can write:


∞
 X
d−n 
2
+
{2d log 2} =

n
n 
n=1
n=d+1


)
( d
∞
 X
X 2d−n mod n
2d−n 
( d)
X 2d−n

=

n=1

n

+



n=d+1

n 

,

where we have inserted modn since were are only interested in the fractional partwhen divided by n. Now note that the numerator 2d−n mod n can be calculated very
rapidly using the binary algorithm for exponentiation, using only standard arithmetic.
Is there a similar formula and scheme for π?

7 / 38

The BBP formula for π
In 1996, a PSLQ program discovered this new formula for π:
π =


∞
X
1
n=0

16k

4
2
1
1
−
−
−
.
8k + 1 8k + 4 8k + 5 8k + 6
This formula permits one to compute binary (or hexadecimal) digits of π beginning at
an arbitrary starting position.
This is the first known instance of a computer program discovering a new formula for π.
BBP-type formulas (also discovered using PSLQ) are now known for numerous othermathematical constants.

I 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.
8 / 38

Some other BBP-type
formulas found using PSLQ

∞
π2 =

1X 1
8
64k
k=0
∞

2 X 1
π =
27
729k
2

144
216
72
54
9
−
−
−
+
(6k + 1)2
(6k + 2)2
(6k + 3)2
(6k + 4)2
(6k + 5)2



k=0

−

243
405
81
27
−
−
−
(12k + 1)2(12k + 2)2(12k + 4)2(27k + 5)2

72
9
9
5
1
−
−
−
+
(12k + 6)2(12k + 7)2(12k + 8)2(12k + 10)2
(12k + 11)2∞

1 X 1
ζ(3) =
1792
212k
k=0







6144
43008
24576
30720
1536
−
+
+
−
(24k + 1)3(24k + 2)3(24k + 3)3(24k + 4)3(24k + 5)3

3072
768
3072
2688
192
1536
+
+
−
−
−
−
(24k + 6)3(24k + 7)3(24k + 9)3(24k + 10)3
(24k + 11)3
(24k + 12)3
96
672
384
24
48
12
−
−
−
+
+
−
(24k + 13)3
(24k + 14)3
(24k + 15)3
(24k + 17)3
(24k + 18)3
(24k + 19)3
+

120
48
42
3
+
−
+
(24k + 20)3
(24k + 21)3
(24k + 22)3
(24k + 23)3



I D. H. Bailey, J. M. Borwein, A. Mattingly and G. Wightwick, “The computation of previously

inaccessible digits of Pi2 and Catalan’s constant,” Notices of the AMS, vol. 60 (2013), 844–854.

9 / 38

Normal numbers
Given integer b ≥ 2, a real number x is b-normal (or “normal base b”) if every m-longstring of digits appears in the base-b expansion of x with limiting frequency 1/b m .
Using measure theory, it can be shown that almost all real numbers are b-normal for agiven integer base b (in fact, for all b simultaneously).
These are widely believed to be b-normal, for all integer bases b ≥ 2 and for all m ≥ 1:
I π = 3.14159265358979323846 . . .
I e = 2.7182818284590452354 . . .
√I
2 = 1.4142135623730950488 . . .
I log(2) = 0.69314718055994530942 . . .
I Every irrational algebraic number (this conjecture is due to Borel).
But there are no proofs of normality for any of these constants. Even the question ofwhether their binary expansions have an equal number of 0s and 1s is unknown.
Until recently, normality proofs were known only for a few contrived examples such as
Champernowne’s constant = 0.12345678910111213141516 . . . (which is 10-normal).
10 / 38

BBP formulas and normality
Let {·} denote fractional part. Consider the sequence x0 = 0 and
1xn = 2xn−1 +
,
n




which is closely connected to the BBP formula for log 2. Richard Crandall and I provedthat log 2 is 2-normal if and only if this sequence is equidistributed in the unit interval.
Similarly, consider the sequence x0 = 0 and(xn =120n2 − 89n + 16
16xn−1 +512n4 − 1024n3 + 712n2 − 206n + 21),

which is closely connected to the BBP formula for π. We proved that π is 16-normal(and hence 2-normal) if and only if this sequence is equidistributed in the unit interval.
I D. H. Bailey and R. E. Crandall, “Random generators and normal numbers,” Experimental

Mathematics, vol. 11, no. 4 (2002), pg. 527–546.
I D. H. Bailey and M. Misiurewicz, “A strong hot spot theorem,” Proceedings of the American

Mathematical Society, vol. 134 (2006), no. 9, pg. 2495–2501.
11 / 38

High-precision numerical integration using the tanh-sinh scheme
Given f (x ) defined on (−1, 1), define g(t) = tanh(π/2 · sinh t). Seting x = g(t) yields
Z ∞

Z 1

f (x ) dx =−10

f (g(t))g (t) dt ≈ h

−∞

N
X

wj f (xj ),j=−Nwhere the abscissas xj and the weights wj are given by:
xj = g(hj) = tanh(π/2 · sinh(hj))
wj = g 0 (hj) = π/2 · cosh(hj)/ cosh2 (π/2 · sinh(hj)).
Features:
I Halving h (i.e., doubling N) typically doubles the number of correct digits.
I The scheme works well even for functions with blow-up singularities at endpoints.
I The cost of computing abscissas and weights increases linearly with the number Nof subdivisions, whereas it increases quadratically for most other schemes.
1. D. H. Bailey, X. S. Li and K. Jeyabalan, “A comparison of three high-precision quadratureschemes,” Experimental Mathematics, vol. 14 (2005), no. 3, pg. 317–329.
2. H. Takahasi and M. Mori, “Double exponential formulas for numerical integration,” Publicationsof RIMS, Kyoto University, vol. 9 (1974), pg. 721–741.

12 / 38

Tanh-sinh quadrature and the Euler-Maclaurin formula
The Euler-Maclaurin formula gives the error in a trapezoidal approximation:Z bf (x ) dx = h
a

−

n
X

h
f (a + jh) − (f (a) + f (b))
2
j=0

m
X
h2i B2i 
i=1

(2i)!



D 2i−1 f (b) − D 2i−1 f (a) − Em (h),where h = (b − a)/n, Bk are Bernoulli numbers, and|Em (h)| ≤ 2(b − a)(h/(2π))2m+2 max |D 2m+2 f (x )|.
a≤x ≤bNote that when f (t) and all of its derivatives are zero and a and b, then the error
Em (h) of a simple block-function approximation to the integral goes to zero more
rapidly than any power of h.
The tanh-sinh substitution x = g(t), where g(t) = tanh(π/2 · sinh t), converts almost
any reasonable integrand into a function that satisfies this condition.
13 / 38

Example of the tanh-sinh transformation
Original function:

3.0

f (x ) = − log cos

πx
2

2.5



2.0

1.5

Note singularities at x = −1 and x = 1.

1.0

0.5

-1.0

0.5

-0.5

1.0

Transformed function:
11
f (g(t))g 0 (t) = − π cosh(t)sech2
π sinh(t)
2
2




1
1
π tanhπ sinh(t)· log cos22The function and all of its derivatives go to zerovery rapidly for large negative or positive t.

0.8

0.6

0.4

0.2

-3

-2

-1

1

2

3

14 / 38

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 line 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.
15 / 38

Limiting value of Cn : What is this number?
Key observation: The Cn integrals can be converted to one- dimensional integrals
involving the modified Bessel function K0 (t):
Cn =

2n
n!

Z ∞
0

tK0n (t) dt

High-precision numerical values, computed using this formula and tanh-sinhquadrature, approach a limit. For example:
C1024 = 0.6304735033743867961220401927108789043545870787 . . .
What is this number? We copied the first 50 digits into the online Inverse Symbolic
Calculator (ISC) at https://isc.carma.newcastle.edu.au. The result was:lim Cn = 2e −2γ .

n→∞where γ denotes Euler’s constant. This is now proven.
16 / 38

Other Ising integral evaluations found using PSLQ
D3 = 8 + 4π 2 /3 − 27 L−3 (2)D4 = 4π 2 /9 − 1/6 − 7ζ(3)/2E2 = 6 − 8 log 2E3 = 10 − 2π 2 − 8 log 2 + 32 log2 2
E4 = 22 − 82ζ(3) − 24 log 2 + 176 log2 2 − 256(log3 2)/3+16π 2 log 2 − 22π 2 /3
E5 = 42 − 1984 Li4 (1/2) + 189π 4 /10 − 74ζ(3) − 1272ζ(3) log 2 + 40π 2 log2 2−62π 2 /3 + 40(π 2 log 2)/3 + 88 log4 2 + 464 log2 2 − 40 log 2
where ζ is the Riemann zeta function and Lin (x ) is the polylog function. E5 remained a“numerical conjecture” for several years, but was proven in March 2014 by Erik Panzer.

17 / 38

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 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 toa PSLQ program. The output integer vector is precisely the set of coefficients of theminimal polynomial of α.
Answer (using 250-digit arithmetic):
0 = 697 − 1440α − 20520α2 − 98280α3 − 102060α4 − 1458α5 + 80α6 − 43920α7+ 538380α8 − 336420α9 + 1215α10 − 80α12 − 56160α13 − 135540α14 − 540α15+ 40α18 − 7380α19 + 135α20 − 10α24 − 18α25 + α30
18 / 38

The Poisson potential function
In 2012, Richard Crandall, while investigating techniques to sharpen images, notedthat each pixel was given by a form of the 2-D Poisson potential function:
1 X cos(mπx ) cos(nπy )φ2 (x , y ) = 2
π m,n odd
m2 + n2
In a 2013 study, we numerically discovered, and then proved the intriguing fact that for
rational numbers x and y ,1φ2 (x , y ) = log απ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 numerous 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.
19 / 38

Samples of minimal polynomials found by 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 + α8

The case x = y = 1/32 required 10,000-digit arithmetic and ran for 45 hours.
Other runs, using even higher precision, ultimately failed. HELP!

20 / 38

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 ,r −1m(s) = 4rY
2(ei −1)pim(pi ).

i=1

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

Three improvements to the Poisson polynomial computation program

1. MPFUN-MPFR: A new thread-safe multiprecision package.
I

Speedup: 3X

2. A new 3-level multipair PSLQ program.
I

Speedup: 4.2X

3. Parallel implementation on a 16-core system.
I

Speedup: 12.2X

Overall speedup: 156X
Total size of code: 190,000 lines

22 / 38

A fast algorithm to compute the Poisson potential function φ2 (x , y )Computing φ2 (x , y ) by the definition (see earlier slide) is much too slow, but
high-precision values can be obtained by this formula:1θ2 (z, q)θ4 (z, q)log
φ2 (x , y ) =,2π
θ1 (z, q)θ3 (z, q)
where q = e −π = 0.0432139 . . ., and z = π2 (y + ix ). The four theta functions can becomputed using the following rapidly convergent formulas:
θ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
k=1

2

q k cos(2kz),
2(−1)k q k cos(2kz).
23 / 38

−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−α192

192-degree minimal polynomial found by multipair PSLQ for x = y = 1/35

This polynomial has degree 192, with coefficients as large as 1085 . This computation
required 18,000-digit arithmetic and 34 CPU-hours.
The case (1/37, 1/37) required 51,000-digit arithmetic and 90 CPU-days (5.6 days on
a 16-core parallel system).
Kimberley’s formula was upheld for (1/s, 1/s), for all s up to 52 (except fors = 41, 43, 47, 49, 51, which were too expensive), and also for s = 60 and s = 64.
24 / 38

Palindromic polynomials
From our results, in the case (1/s, 1/s) where s is even, the resulting polynomial isalways 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, an undergraduate student at Stanford University, observed that if α is aroot of a palindromic polynomial such as this, then α + 1/α is a root of a transformed
polynomial of half the degree.
25 / 38

Proofs of Kimberley’s formula and the palindromic property
I On March 16, DHB presented our results at a seminar at the University ofCalifornia, Berkeley.
I Following the presentation, Watson Ladd, a graduate student in mathematics,

brought to our attention the fact that some of our conjectures should follow fromresults in the theory of elliptic curves, Gaussian integers and ideals.
I After some effort, Ladd produced proofs of Kimberley’s formula and thepalindromic property, which proofs were then included in our paper.

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

26 / 38Mordell-Tornheim-Witten sums and the omega function
Mordell-Tornheim-Witten sums are defined as:
W (r1 , r2 , . . . , rn , t) =X1.
mr1 mr2 · · · mnrn (m1 + m2 + · · · + mn )t
m1 ,...,mn ≥1 1 2They can be continued analytically, as specified by Matsumoto.
Such sums arise in combinatorics, mathematical physics (e.g., Feynmann diagrams andstring theory), Lie algebras, number theory and numerous other fields.
We will focus on the special case r1 = r2 = · · · = rn = t = s for real s (analyticallycontinued as above), namelyX1.
ωn+1 (s) =(m1 m2 · · · mn (m1 + m2 + · · · + mn ))sm ,m ,...,m ≥1
1

2

n

I D. Romik, “On the number of n-dimensional representations of SU(3), the Bernoulli numbers,

and the Witten zeta function,” 17 Apr 2015, https://arxiv.org/pdf/1503.03776v3.pdf.
I H. Tomkins, “An exploration of multiple zeta functions,” Honours Thesis, Dalhousie University,

21 Apr 2016.

27 / 38

The degree-two omega function
The degree-two instance of the omega function isω(s) =∞
X1.
(jk(j + k))sj,k=1

For positive integer arguments, papers by Subbarao, Sitaramachandrarao, Garoufalidis,
Zagier, Huard, Williams and Nan-Yue have shown that for positive integers n,n4n − 2k − 1
4X
ζ(2k)ζ(6n − 2k)
ω(2n) =
3 k=0
2n − 1!

ω(2n + 1) = −4

n
X
k=0

!

4n − 2k − 1
ζ(2k)ζ(6n − 2k + 3).
2n

28 / 38

Specific values of omega
ω(1) = 2ζ(3)1 6
ω(2) =
π2835
ω(3) = −2π 2 ζ(7) + 20ζ(9)19ω(4) =
π 12273648375
2
70
ω(5) = − π 4 ζ(11) − π 2 ζ(13) + 252ζ(15)932062
ω(6) =
π 18116937886440375
4 656
ω(7) = −π ζ(15) − π 4 ζ(17) − 308π 2 ζ(19) + 3432ζ(21)
135
15
32293
ω(8) =
π 24
7060124893837640625
29 / 38

Romik’s formulas for omega
Dan Romik has shown that1
ω(s) =
2πΓ(s)Z

Γ(s + z)Γ(−z)ζ(2s + z)ζ(s − z) dz,(α)

where (α) means the vertical line {α + ti, −∞ < t < ∞}. This is valid for all complex
s satisfying max(−Re(s), 1 − 2Re(s)) < α < min(0, Re(s) − 1).
1
Γ(2s − 1)Γ(1 − s)ζ(3s − 1)Γ(s)ω(s) =
+
+

M−1
X

(−1)k
Γ(s + k)ζ(2s + k)ζ(s − k)
k!
k=0

1
2π

Z



Γ(s + z)Γ(−z)ζ(2s + z)ζ(s − z) dz .

(α)

This is valid for all complex s satisfying 3/4 − M/2 < Re(s) < M + 1/2.
30 / 38

The Riemann zeta function in the complex plane

31 / 38

The omega function

32 / 38

Research questions for omega
I Does omega have any zeroes in the complex plane? Where? How many?
I Does omega have nice analytic formulas for any arguments other than positive

integers? Can we discover these formulas experimentally?
I Can we do a similar analysis of the derivative of omega? For higher derivatives?
I How about for higher-degree omega?

The argument principle
Note, by the argument principle, that the number of zeroes of ω(s) in a pole-free
region bounded by a contour C in the complex plane can be computed as:
−

1
2πi

ω 0 (s)ds(C ) ω(s)Z

33 / 38

Computing the omega function
Key elements of computing ω(s) for complex s:
I Evaluating Γ(s).
I Evaluating ζ(s).
I Evaluating Romik’s integral formula.

To count zeroes in a region using the argument principle, it is also necessary tocompute the derivative ω 0 (s). This can be done byω 0 (s) ≈ω(s + h) − ω(s),hwhere h = 10−p/2 and p is the working precision level in digits (yielding the derivative
accurate to p/2 digits).
The reliability of these results strongly depends on performing the computations toconvincingly high precision — more than standard 64-bit floating-point arithmetic.
34 / 38

Computing gamma and zeta
The best way to compute Γ(s) is to use a formula originally due to Spouge:




dae−1
X ck (a)√,Γ(1 + z) ≈ (z + a)z+1/2 e −z−a 2π 1 +
z
+
k
k=1

where a = number of desired digits, andck (a) =(−1)k−1 (a − k)k−1/2 e −k+a√.
2π(k − 1)!

The best way to compute ζ(s) is to use a formula due to Richard Crandall:
∞
X
1
π s/2
Γ(s/2, πn2 ) π s−1/2 Γ((1 − s)/2, πn2 )ζ(s) =++
Γ(s/2) s(s − 1) n=1
ns
n1−s

"

!#

,

where Γ(s, x ) is the incomplete gamma function:s −xΓ(s, x ) = Γ(s) 1 − x e

∞
X

!

xk
.
Γ(s + k + 1)
k=0

35 / 38

Software for high-precision computation
I Mathematica and Maple: Both feature well-integrated arbitrary precision facilities,

but are often quite slow.
I Sage: Also includes arbitrary precision facility, but not as many functions.
I gfortran’s quad precision: Now includes full support for real(16) andcomplex(16) datatypes (128-bit IEEE arithmetic, approximately 33 digits).
I GMP: GNU-supported package for arbitrary precision integer computation.
I MPFR: Very efficient package for arbitrary precision floating-point computation,

built on top of GMP; does not include high-level language facilities.
I DHB’s software (each includes high-level language support):
QD: Double-double (31 digits) and quad-double (62 digits); C++ and Fortran.
ARPREC: Arbitrary precision floating-point package; C++ and Fortran.
I MPFUN2015: New Fortran arbitrary precision package, based on MPFR; Fortran.
II

36 / 38

What has been done

30

These have been implemented, using the MPFUN2015
multiprecision software:
25

I Routine to compute Γ(s) for any complex s.
I Routine to compute Γ(s, x ) for any complex s and all

real x .

20

I Routine to compute ζ(s) for any complex s.
I Routine to compute ω(s) for any complex s, using

15

Romik’s integral formulas and the tanh-sinh quadraturealgorithm (using MPFUN2015).
I Computation of the contour integral for the rectangular

10

contour (0, 1/8) to (2, 1/8) to (2, 32) to (0, 32), andback to (0, 1/8), all using 70-digit precision.
5

I The contour integral calculation has been running 20

days and is still not yet complete!
0.5

1.0

1.5

2.0

37 / 38

What remains to be done
I Run times are ridiculously high. Can the program be tuned to run faster? Use

different high-precision software? Run on a parallel system?
I Can a much larger range of complex arguments be explored, looking for zeroes,

perhaps not necessarily using the argument principle?
I Can we do the same for the derivative of omega or higher-order derivatives? Can

we do the same for higher-degree omega?
Any volunteers?
Contact me at: dhbailey@ucdavis.edu
This talk is available at:
http://www.davidhbailey.com/dhbtalks/dhb-ucdavis-omega.pdf
38 / 38