Computation of the omega function
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Laboratory (retired)
University of California, Davis, Department of Computer Science1 / 17

Mordell-Tornheim-Witten sums and the omega function
Mordell-Tornheim-Witten sums are defined as:
W (r1 , r2 , . . . , rn , t) =X

1.
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
12nI 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.

2 / 17

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) = −4nXk=0!

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

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

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). This can be written,avoiding singularities, as:1
Γ(2s − 1)Γ(1 − s)ζ(3s − 1)
ω(s) =
Γ(s)++M−1X(−1)kΓ(s + k)ζ(2s + k)ζ(s − k)k!
k=012π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.
5 / 17

The Riemann zeta function in the complex plane

6 / 17

The omega function on [x + iy , 0 ≤ x ≤ 2, 1/8 ≤ y ≤ 32]

7 / 17

The omega function on [x + iy , −2 ≤ x ≤ 2, 1/8 ≤ y ≤ 64]

8 / 17

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 positiveintegers? 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?

9 / 17

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

10 / 17

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.
11 / 17

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+kk=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:∞X1π s/2
Γ(s/2, πn2 ) π s−1/2 Γ((1 − s)/2, πn2 )ζ(s) =++Γ(s/2) s(s − 1) n=1
nsn1−s"!#,where Γ(s, x ) is the incomplete gamma function:s −xΓ(s, x ) = Γ(s) 1 − x e∞X!

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

12 / 17

Computing the integral in Romik’s formula 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 1f (x ) dx =−10f (g(t))g (t) dt ≈ h−∞NX

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 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.
13 / 17

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

14 / 17What has been done30

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 allreal 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 rectangular10contour (0, 1/8) to (2, 1/8) to (2, 64) to (0, 64), andback to (0, 1/8), all using 70-digit precision.
50.5

1.0

1.5

2.0

15 / 17

Zeroes of omega in the region [x + iy , −2 ≤ x ≤ 2, 1/8 ≤ y ≤ 64]
6050
4030
20
100-2 -101216 / 17

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?

Contact me at: david@davidhbailey.com
This talk is available at:
http://www.davidhbailey.com/dhbtalks/dhb-wcnt-2017.pdf
17 / 17