Computation and theoryof extended Mordell-Tornheim-Witten sums
David H. Bailey
http://www.davidhbailey.com
Computational Research Dept., Lawrence Berkeley National Lab (retired)Computer Science Department, University of California, Davis

14 Dec 2013

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 in the operation of thealgorithm.
1. 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.
2. D.H. Bailey and D.J. Broadhurst, “Parallel Integer Relation Detection: Techniques and Applications,” Mathematics of Computation, vol. 70,
no. 236 (Oct 2000), pg. 1719–1736.

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 wasdesigned for parallel systems, but runs faster even on 1 CPU.

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

The first major PSLQ discovery: The BBP formula for π
In 1996, a PSLQ program discovered this new formula for π:π =
∞
X
1
4
2
1
1−−−
16k 8k + 1 8k + 4 8k + 5 8k + 6n=0

This formula permits one to compute binary (or hexadecimal) digits of π beginning at
an arbitrary starting position, using a very simple scheme that requires only standard
64-bit or 128-bit arithmetic.
In 2004, Borwein, Galway and Borwein proved that no base-n formulas of this type
exist for π, except when n = 2m .
BBP-type formulas (discovered with PSLQ) are now known for numerous othermathematical constants.
1. D.H. Bailey, P.B. Borwein and S. Plouffe, “On the rapid computation of various polylogarithmic constants,” Mathematics of Computation,
vol. 66, no. 218 (Apr 1997), pg. 903–913.
2. J.M. Borwein, W.F. Galway and D. Borwein, “Finding and excluding b-ary Machin-type BBP formulae,” Canadian Journal of Mathematics,
vol. 56 (2004), pg. 1339–1342.

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

ζ(3)=2167254
9144−−−+(6k + 1)2(6k + 2)2(6k + 3)2(6k + 4)2(6k + 5)2
k=0

∞
2 X 1
405
81
27
243−−−
27
729k (12k + 1)2(12k + 2)2(12k + 4)2(27k + 5)2
k=0

995
1
72−
−
−+−
(12k + 6)2(12k + 7)2(12k + 8)2(12k + 10)2(12k + 11)2

∞
1 X 1
6144
43008
24576
30720
1536
−++−
1792
212k (24k + 1)3(24k + 2)3(24k + 3)3(24k + 4)3(24k + 5)3
1X 1
864kk=0307276830722688
192
1536+−−−−+(24k + 6)3(24k + 7)3(24k + 9)3(24k + 10)3(24k + 11)3(24k + 12)396
6723842448
12−−−++−(24k + 13)3
(24k + 14)3(24k + 15)3(24k + 17)3(24k + 18)3(24k + 19)3

120
4842
3++−+(24k + 20)3(24k + 21)3(24k + 22)3(24k + 23)3
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 American Mathematical Society, vol. 60 (2013), no. 7, pg. 844-854.

Algebraic numbers in Poisson potential functions and lattice sums
Lattice sums arising from the Poisson equation have been studied widely inmathematical physics and also in image processing. We numerically discovered, andthen proved, that for rational (x, y ), the two-dimensional Poisson potential function
satisfies
11 X cos(mπx) cos(nπy )=log αφ2 (x, y ) =222πm +nπm,n odd

where α is an algebraic number, i.e., the root of an integer polynomial
0 = a0 + a1 α + a2 α 2 + · · · + an α nThe minimal polynomials for these α were found by PSLQ calculations, with the(n + 1)-long vector (1, α, α2 , · · · , αn ) as input, where α = exp(8πφ2 (x, y )). PSLQ
returned the vector of integer coefficients (a0 , a1 , a2 , ·, an ) as output.
1. 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), pg. 115201, http://www.davidhbailey.com/dhbpapers/PoissonLattice.pdf.
2. D.H. Bailey and J.M. Borwein, “Compressed lattice sums arising from the Poisson equation: Dedicated to Professor Hari Sirvastava,”
Boundary Value Problems, vol. 75 (2013), DOI: 10.1186/1687-2770-2013-75,
http://www.boundaryvalueproblems.com/content/2013/1/75.

Samples of minimal polynomials found by PSLQ
k5678910

Minimal polynomial for exp(8πφ2 (1/k, 1/k))
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 + α8The minimal polynomial for exp(8πφ2 (1/32, 1/32)) has degree 128, with individual
coefficients ranging from 1 to over 1056 . This PSLQ computation required 10,000-digit
precision. See next page.
Other polynomials required up to 50,000-digit precision.

Degree-128 minimal polynomial for exp(8πφ2 (1/32, 1/32))

High-precision tanh-sinh numerical integration
Given f (x) defined on (−1, 1), define g (t) = tanh(π/2 sinh t). Then setting x = g (t)yieldsZ ∞Z 1
f (x) dx−1=−∞f (g (t))g 0 (t) dt ≈ hNX

wj f (xj ),j=−Nwhere xj = g (hj ) and wj = g 0 (hj ). Since g 0 (t) goes to zero very rapidly for large t, theproduct f (g (t))g 0 (t) typically is a nice bell-shaped function, so that the simple
summation above converges very rapidly. Reducing h by half typically doubles thenumber of correct digits.
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.
1. D.H. Bailey, X.S. Li and K. Jeyabalan, “A Comparison of Three High-Precision Quadrature Schemes,” Experimental Mathematics, vol. 14(2005), no. 3, pg. 317–329.
2. H. Takahasi and M. Mori, “Double Exponential Formulas for Numerical Integration,” Publications of RIMS, Kyoto University, vol. 9 (1974),
pg. 721-741.

Extended Mordell-Tornheim-Witten functions
In this study, we examineωs1 , . . . , sM
d1 , . . . , dM

| t1 , . . . , tN
| e1 , . . . eN

X:=MNY(− log mj )dj Y (− log nk )ekmj sj

m1 ,...,mM ,n1 ,...,nN > 0 j=1PMPm = N
nj=1 jk=1 k=

12πZ 2π YM
0j=1(d )Lisj j



eiθ

N
Y

k=1(e )

Litk knk tk

e −iθ dθ,k=1

(d)∂ dLis (z).
where the s-th outer derivative of a polylogarithm is denoted Lis (z) := ∂sThese are a generalization of the class of Mordell-Tornheim-Witten sums, which appear
in many contexts, including combinatorics, number theory and mathematical physics.
1. D.H. Bailey and J.M. Borwein, “On Eulerian log-gamma integrals and Tornheim-Witten zeta functions,” Ramanujan Journal, 27 Feb 2013,
available at http://www.davidhbailey.com/dhbpapers/log-gamma.pdf.
2. D.H. Bailey, J.M. Borwein and A. Mattingly, “Computation and theory of extended Mordell-Tornheim-Witten sums II,” manuscript, available
from author.

Computation of polylogarithm outer derivatives to high precision
Fix k = 0, 1, 2 . . . and m = 1, 2 . . . . For | log z| < 2π and L = log(− log z)(m)

Lik+1 (z) =P0≤n6=k ζnklog z(m) (k + 1 − n) log z + m! ck,m (L) k! .
n!

Here, for k ≥ 1(−1)jγj − bk,j+1 (L)j!
X L p Γ(t) (1)bk,m (L) :=(−1)t fk,q ,p!
t!
p+t+q=jck,j (L)=p,t,q≥0

where γj are Stieltjes gamma constants, and fk,q are definedrecursively byPf0,0 = 1, f0,q = 0 (q > 0), fk,0 = 1 (k > 0) and fk,q = qh=0 (−1)h fk−1,q−h /k h .
Additional computational machinery is required to compute derivatives of zeta function
at negative integers and derivatives of gamma function at 1.

Why not just use Maple or Mathematica?
Performing the computations just described is very demanding. We use custom-written
code and the ARPREC high-precision arithmetic package (due to DHB and others).
Neither Mathematica nor Maple is up to the task:
1. Maple cannot compute numerical values of outer derivatives of polylogarithms.
2. Mathematica computes numerical values of outer derivatives of polylogarithms,
but in some cases they are not correct to the specified number of digits.
3. Mathematica refuses to compute ζ (m) (0) for m >= 4. It returns the expression
unevaluated. Wolfram Alpha produces analytic expressions for these, but will only
give numerical to a few digits unless one signs up for Wolfram Alpha Pro.
4. Mathematica’s computations of ζ (m) (−n) to high precision often produce strange
internal errors.
5. Our custom-written code, using ARPREC high-precision software is many times
faster than either Maple or Mathematica for performing high-precision numericalintegration (quadrature).

Order three omega constants
For nonnegative a, b, c, q, r , s, we focus on this particular class of omega constants:
ωa,b,c (q, r , s) = ωq, r | sa, b | cZ ∞  s−1 (c)x(a)(b)=Liq (e −x ) Lir (e −x ) dx,Γ(s)0which is valid when q ≥ 0, r ≥ 0, s > 0, with q + r + s > 2, and a ≥ 0, b ≥ 0, c ≥ 0.
Here the notation (·)(c) denotes the c-th partial derivative of the expression inparentheses with respect to s. For computational purposes, it is necessary to break theintegral into two pieces, and then using the substitution u = e −x for the second:Z 1  s−1 (c)x(a)(b)ωa,b,c (q, r , s) =Liq (e −x ) Lir (e −x ) dx
Γ(s)0
(c)
Z 1/e (− log u)s−1
du(a)
(b)
=+Liq (u) Lir (u) .
Γ(s)u0Small sample of relations discovered by PSLQ among omega constants
0 = ω0,0,3 (0, 1, 2) − ω0,3,0 (0, 1, 2) − ω0,3,0 (1, 0, 2) + ζ (3) (3)0 = ω0,1,1 (0, 1, 3) + ω0,1,1 (1, 0, 3) + ω0,1,1 (2, 0, 2) − ω0,1,1 (2, 1, 1)
0 = −ω0,2,1 (0, 1, 3) − ω0,2,1 (0, 2, 2) − ω0,2,1 (1, 0, 3) + ω0,2,1 (1, 2, 1)?

0 = ω0,0,3 (0, 1, 3) − ω0,1,2 (0, 2, 2) − ω0,2,1 (0, 2, 2) − ω0,3,0 (0, 1, 3) − ω0,3,0 (0, 2, 2)−ω0,3,0 (1, 0, 3) + 2ω1,2,0 (0, 1, 3) + ω1,2,0 (0, 2, 2) + 2ω1,2,0 (1, 0, 3) + ω1,2,0 (2, 0, 2)
0 = ω1,3,0 (0, 1, 3) + ω1,3,0 (1, 0, 3) + ω1,3,0 (2, 0, 2) − ω1,3,0 (2, 1, 1)?

0 = −2ω0,0,4 (0, 1, 3) − ω0,0,4 (0, 2, 2) + 4ω0,4,0 (0, 1, 3) + 2ω0,4,0 (0, 2, 2) + 4ω0,4,0 (1, 0, 3)+ω0,4,0 (2, 0, 2) − 3ζ (4) (4)?

0 = −ω1,2,0 (0, 1, 4) − ω1,2,0 (1, 0, 4) − ω1,2,0 (2, 0, 3) + ω1,2,0 (2, 1, 2)?

0 = −3ω0,0,4 (0, 1, 4) − ω0,0,4 (0, 2, 3) + 6ω0,4,0 (0, 1, 4) + 5ω0,4,0 (0, 2, 3) + 3ω0,4,0 (0, 3, 2)+6ω0,4,0 (1, 0, 4) + ω0,4,0 (2, 0, 3) − 4ζ (4) (5)Provable relations among the order-three omegas
The omegas satisfy relations such asωa,b,c (q, r , s) = ωb,a,c (r , q, s)ωa,b,c (r , s, t − 1) = ωa,b,c (r − 1, s, t) + ωa,b,c (r , s − 1, t)
ωa,b,0 (s, t, 0) − ωa,0,b (t, 0, s) − ωa,0,b (s, 0, t) = ζ (a+b) (t + s)
r Xr +s −i −1
ωa,b,c (r , s, t) =ωa,b,c (i, 0, N − i)s −1i=1
s X
r +s −i −1+
ωa,b,c (0, i, N − i)
r −1
i=1

We conjecture that all relations found by PSLQ are generated by these relations.

Summary

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