# One of the numbers $\zeta(5), \zeta(7), \zeta(9), \zeta(11)$ is irrational

I am reading an interesting paper One of the numbers ζ(5), ζ(7), ζ(9), ζ(11) is irrational by Zudilin. We fix odd numbers $$q$$ and $$r$$, $$q\geq r+4$$ and a tuple $$\eta_0,\eta_1,...,\eta_q$$ of positive integer parameters satisfying the conditions $$\eta_1\leq \eta_2\leq...\leq \eta_q<\eta_0/2$$ and $$\eta_1+\eta_2+...+\eta_q\leq \eta_0\left(\frac{q-r}{2}\right)\tag{1}$$ Define $$F_n:=\frac{1}{(r-1)!}\sum_{t=0}^\infty R_n^{(r-1)}(t)\tag{2}$$ and note that $$R_n(t)=O(t^{-2})$$. We put $$m_j=\max\{\eta_r,\eta_0-2\eta_{r+1},\eta_0-\eta_1-\eta_{r+j}\}$$ for $$j=1,2,...,q-r$$ and define the integer $$\Phi_n:=\prod_{\sqrt{\eta_0 n} where only primes enter the product and $$\varphi(x)=\min_{0\leq y<1}\left(\sum_{j=1}^{r}([y]+[\eta_0x-y]-[y-\eta_j x]-[(\eta_0-\eta_j)x-y]-2[\eta_j x])+\sum_{j=r+1}^{q}([(\eta_0-2\eta_j)x]-[y-\eta_j x]-[(\eta_0-\eta_j)x-y])\right)$$ where [.] denotes the ceiling function. Let $$D_N$$ denote the lcm of $$1,2,...,N$$.

Lemma $$1$$: ($$2$$) defines a linear form of $$\zeta(r+2),\zeta(r+4),...,\zeta(q-2)$$ with rational coefficients; moreover, $$D_{m_1n}^r D_{m_2n... D_{m_{q-r}n}}.\Phi_n^{-1}.F_n\in\mathbb{Z}+\mathbb{Z}\zeta(r+2)+\mathbb{Z}\zeta(r+4)+...+\mathbb{Z}\zeta(q-2) \tag{3}$$ By Prime Number Theorem, $$\lim_{n\to\infty}\frac{\log D_{m_j n}}{n}=m_j,\ \ \ j=1,...,q-r.$$ We introduce the auxiliary function $$f_0(\tau)=r\eta_0\log(\eta_0-\tau)+\sum_{j=1}^{q} (\eta_j\log(\tau-\eta_j)-(\eta_0-\eta_j)\log(\tau-\eta_0+\eta_j)) -2\sum_{j=1}^r \eta_j\log \eta_j+\sum_{j=r+1}^q (\eta_0-2\eta_j)\log(\eta_0-2\eta_j)$$ defined in the $$\tau$$-plane with the cuts $$(-\infty,\eta_0-\eta_1]$$ and $$[\eta_0,+\infty)$$

Lemma $$2$$: Let $$r=3$$ and $$\tau_0$$ be a zero of the polynomial $$(\tau-\eta_0)^r(\tau-\eta_1)...(\tau-\eta_q)-\tau^r(\tau-\eta_0+\eta_1)...(\tau-\eta_0+\eta_q)$$ with Im $$\tau_0>0$$ and the maximum possible value of Re $$\tau_0$$. Assume Re $$\tau_0<\eta_0$$ and Im $$f_0(\tau_0)\notin \pi \mathbb{Z}.$$ Then $$\overline{\lim}_{n \to \infty} \frac{\log |F_n|}{n}=Re f_0(\tau_0)$$ If the sequence of linear forms on the left side of ($$3$$) assumes non-zero arbitrarily small values as $$n$$ increases, then in the case $$r=3$$ there are irrational numbers among $$\zeta(5),\zeta(7),...,\zeta(q-4),\zeta(q-2) \tag{4}$$ Therefore the following holds:

Lemma $$3$$: Suppose that $$r=3$$ and in the above notation $$C_0=-Re f_0(\tau_0)$$, $$C_1=rm_1+m_2+...+m_{q-r}-\left(\int_0^1\varphi(x)\frac{\Gamma'(x)}{\Gamma(x)}dx-\int_0^{1/m_{q-r}}\varphi(x)\frac{dx}{x^2}\right)$$ If $$C_0>C_1$$, then at least of the numbers ($$4$$) is irrational. The line below Lemma $$3$$ reads: we put, $$r=3,q=13$$, $$\eta_0=91,\eta_1=\eta_2=\eta_3=27,\eta_4=29,\eta_5=30,\eta_6=31,...,\eta_{12}=37,\eta_{13}=38.$$ Then $$C_0=227.58019641...,C_1=226.24944266...$$

Question: Is the above result true only for $$r=3$$ or for other odd values of $$r>3$$ also? For example can we take $$r=5,q=13$$ and if we show that $$C_0>C_1$$ and with suitable choice of $$\eta's$$, can we say at least one of the four numbers $$\zeta(7),\zeta(9),\zeta(11),\zeta(13)$$ is irrational? Rephrasing my question, Is Proposition $$5$$ in Arithmetic of linear forms involving odd zeta values true for $$r=5$$ and $$q=13$$?

Any help would be highly appreciated. Thank you!

• Have you tried Bernstein polynomial ? Commented Mar 30 at 8:29

There isn't a lot of thought that goes into the verification. It's pure computation at this point.

I'll include code that computes $$C_0$$. Computing $$C_1$$ is similarly a direct computation (but it requires integrating the digamma function).

# This is sage code, for the Sagemath computer algebra system.
# It's very similar to python, but with extra batteries included.

etas = [91, 27, 27, 27, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]
r = 3
q = 13

term1(x) = 1
term2(x) = 1

for eta in etas:
term1 *= (x - eta)
term1 *= (x - etas[0])^2

for eta in etas[1:]:
term2 *= (x - etas[0] + eta)
term2 *= x^3


Now term1 and term2 are the two polynomials

$$(x - 38) \cdot (x - 37) \cdot (x - 36) \cdot (x - 35) \cdot (x - 34) \cdot (x - 33) \cdot (x - 32) \cdot (x - 31) \cdot (x - 30) \cdot (x - 29) \cdot (x - 91)^{3} \cdot (x - 27)^{3}$$ and $$(x - 62) \cdot (x - 61) \cdot (x - 60) \cdot (x - 59) \cdot (x - 58) \cdot (x - 57) \cdot (x - 56) \cdot (x - 55) \cdot (x - 54) \cdot (x - 53) \cdot (x - 64)^{3} \cdot x^{3}.$$ (I just had sage output the tex for these polynomials for me).

We continue.

f = term1 - term2
f = CC[x](f)  # interpret f as a polynomial over complex numbers

roots = f.roots()
# print(roots)


The roots are

 (-14.1951485567264),
(45.4999978623522),
(105.195148556996),
(3.52099458165416 - 3.32820690552685*I),
(3.52099458165416 + 3.32820690552685*I),
(45.4999998458723 - 7.26440744326658*I),
(45.4999998458723 + 7.26440744326658*I),
(45.4999999754805 - 13.0899280357995*I),
(45.4999999754805 + 13.0899280357995*I),
(45.5000000003980 - 24.9895408343146*I),
(45.5000000003980 + 24.9895408343146*I),
(45.5000012455796 - 3.29700826688036*I),
(45.5000012455796 + 3.29700826688036*I),
(87.4790054197046 - 3.32820691055980*I),
(87.4790054197046 + 3.32820691055980*I).


Continuing,

cand_roots = [r for r, deg in roots if imag(r) > 0]
cand_roots.sort(key = lambda z: real(z))
tau0 = cand_roots[-1]


The root $$\tau_0$$ is given by

$$\tau_0 = 87.4790054197046 + 3.32820691055980i.$$

Finally, we directly compute $$f_0$$.

def f0(tau):
ret = r * etas[0] * log(etas[0] - tau)
for j in range(1, q + 1):
term = etas[j] * log(tau - etas[j])
term -= (etas[0] - etas[j]) * log(tau - etas[0] + etas[j])
ret += term
for j in range(1, r + 1):
term = 2 * etas[j] * log(etas[j])
ret -= term
for j in range(r + 1, q + 1):
term = (etas[0] - 2 * etas[j]) * log(etas[0] - 2 * etas[j])
ret  += term
return ret

print(-real(f0(tau0).n()))
--
227.580196314623


I observe that there are some differences in the less significant digits. If I use 200 bits of precision (instead of the default 53 bit double precision implicit above) ((second parenthetical: this can be done by using f = ComplexField(200)[x](f) instead of CC[x](f) above)), then I get

$$C_0 \approx 227.58019641270392421956170302923769732274712647908924133109.$$

This agrees with the claimed $$C_0$$. It is also possible to bound the numerical error due to precision, but that's a separate topic.

I looked at computing $$C_1$$ again today and give my incomplete code. When I wrote the first version of this answer, I tried to be clever about computing $$\varphi(x)$$ (which is merely a somewhat complicated step function, and which we can in principle identify what intervals it takes which values).

Today, I'm not clever and I'll give a brute force answer. I actually get something very close to what I had last week (perhaps I'm consistently wrong somewhere, in which case I hope this partial code leads someone to identify my mistake).

To compute $$\varphi(x)$$, I look on a mesh of $$10000$$ different $$y$$ values and identify the minimum.

ms = [0]
for j in range(1, q - r + 1):
mj = max(etas[r],
etas[0] - 2 * etas[r + 1],
etas[0] - etas[1] - etas[r + j])
ms.append(mj)

first_term = r * ms[1]
for j in range(2, q - r + 1):
first_term += ms[j]

def brutephi(x):
def innerphi(y):
inner_ret = 0
for j in range(1, r + 1):
term = floor(y) + floor(etas[0] * x - y) - floor(y - etas[j] * x) \
- floor( (etas[0] - etas[j])*x - y) - 2 * floor(etas[j] * x)
inner_ret += term
for j in range(r + 1, q + 1):
term = floor((etas[0] - 2 * etas[j]) * x) - floor(y - etas[j]*x) \
- floor((etas[0] - etas[j]) * x - y)
inner_ret += term
return inner_ret
m = innerphi(0.5)
LIMIT = 100**2
for y in range(LIMIT):
cand = innerphi(y / LIMIT)
m = min(m, cand)
if m == 0:
return 0
return m


This is just an approximation. (Actually, I populated a table with precomputed values for numerical integration. I used something equivalent to brutephi to populate the table).

For the first integral, we note that $$\psi$$ is an increasing, differentiable function, and thus

$$\int_0^1 \varphi(x) d \psi = \int_0^1 \varphi(x) \psi'(x) dx.$$

The correct way to do this would be to break $$[0, 1]$$ into the intervals where $$\varphi(x)$$ is constant and explicitly integrate $$\psi'(x)$$ on each of those intervals (which is just $$\psi(x)$$). But again, today we use brute force and I use numerical Riemann-Stieltjes integration. (Aside: I also implemented the trapezoidal rule in Riemann integration, and it gave me the same answer).

RRMY = RealField(200)
def integral1brute(TERMS):
ret = RRMY(0)
for i in range(1, TERMS - 1):
j = RRMY(i + 1)
xi = RRMY(1. * i / TERMS)
xj = RRMY(1. * j / TERMS)
mid = RRMY((xi + xj)/2.)
height = RRMY(brutephi(mid))
width = RRMY(psi(xj) - psi(xi))
ret += RRMY(n(height * width))
return ret


The second integral is similar. I use the trapezoidal rule for numerical integration, explicitly.

def integrand2(x):
x = RRMY(x)
return brutephi(x) / RRMY(x^2)

def integral2(TERMS):
ret = RRMY(0)
width = RRMY(((1./ms[q-r]) - 0) / TERMS / 2.)
# 0th is 0, so we ignore it
for i in range(1, TERMS + 1):
x = RRMY(0 + (1 - (1/ms[q-r])) * i / TERMS)
height = integrand2(x)
ret += RRMY(n(2 * height))
ret += integrand2(1./ms[q-r])
ret *= RRMY(width)
return ret


We put this all together.

def compC1(TERMS):
ft = first_term
i1 = integral1brute(TERMS)
i2 = integral2(TERMS)
print(ft, i1, i2, "C1 is approx", ft - i1 + i2)

compC1(10000)
# 403
# 226.88252
# 6.9549373
# C1 is approx 183.07241


This approximation says that $$C_1 \approx 183$$, which is rather far from the claimed value. I note that brute-forcing $$\varphi$$ in the way I've done it should overestimate the first integral, which should result in subtracting too large of a number, and so any errors from $$\varphi$$ computation all make the resulting $$C_1$$ estimate too small.

I would be surprised if the brute approximation to $$\varphi$$ was enough to cause such a big difference. Nonetheless, I hope that this incomplete solution is helpful.

• ($+1$). Thank you so much for the answer. Actually I was more curious to calculate $C_1$. How do we find $C_1$? Please edit your answer. I will with utmost respect accept it.
– Max
Commented Mar 22 at 1:46
• I tried your code on CoCalc online. I am getting the following error: ERROR: "failed to get sage socket". I am also getting "NameError: name 'tau0' is not defined". I am also getting "SyntaxError: cannot assign to function call here. Maybe you meant '==' instead of '='?"
– Max
Commented Mar 22 at 4:05
• failed to get sage socket is a cocalc problem. Not having tau0 is a sign you didn't include the line tau0 = cand_roots[-1]. I don't know what line you're referring to for the other. I'm using sage 10.2 for this, and I just reran the code here to make sure it works. If you have a different computing environment, it's probably not so hard to make appropriate adjustments to make it work. Commented Mar 22 at 17:21
• Computing $C_1$ is not so bad (I think), but I just wrote some code to do it and I currently don't get the same answer as claimed. I'll look at it again later. Commented Mar 22 at 17:31
• I need to compute $C_1$. I am not that friendly with sage. Any help would be highly appreciated.
– Max
Commented Mar 28 at 12:48

Too long for a comment. I have entered the code given by @davidlowryduda on sagemath as:

# This is sage code, for the Sagemath computer algebra system.
# It's very similar to python, but with extra batteries included.
etas = [91, 27, 27, 27, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]
r = 3
q = 13

term1(x) = 1
term2(x) = 1

for eta in etas:
term1 *= (x - eta)
term1 *= (x - etas[0])^2

for eta in etas[1:]:
term2 *= (x - etas[0] + eta)
term2 *= x^3
f = term1 - term2
f = CC[x](f)  # interpret f as a polynomial over complex numbers

roots = f.roots()
# print(roots)
cand_roots = [r for r, deg in roots if imag(r) > 0]
cand_roots.sort(key = lambda z: real(z))
tau0 = cand_roots[-1]
def f0(tau):
ret = r * etas[0] * log(etas[0] - tau)
for j in range(1, q + 1):
term = etas[j] * log(tau - etas[j])
term -= (etas[0] - etas[j]) * log(tau - etas[0] + etas[j])
ret += term
for j in range(1, r+1 ):
term = 2 * etas[j] * log(etas[j])
ret -= term
for j in range(r + 1, q+1):
term = (etas[0] - 2 * etas[j]) * log(etas[0] - 2 * etas[j])
ret  += term
return ret

print(-real(f0(tau0).n()))
ms = [0]
for j in range(1, q - r + 1):
mj = max(etas[r],
etas[0] - 2 * etas[r + 1],
etas[0] - etas[1] - etas[r + j])
ms.append(mj)

first_term = r * ms[1]
for j in range(2, q - r + 1):
first_term += ms[j]

def brutephi(x):
def innerphi(y):
inner_ret = 0
for j in range(1, r + 1):
term = floor(y) + floor(etas[0] * x - y) - floor(y - etas[j] * x) \
- floor( (etas[0] - etas[j])*x - y) - 2 * floor(etas[j] * x)
inner_ret += term
for j in range(r + 1, q + 1):
term = floor((etas[0] - 2 * etas[j]) * x) - floor(y - etas[j]*x) \
- floor((etas[0] - etas[j]) * x - y)
inner_ret += term
return inner_ret
m = innerphi(0.5)
LIMIT = 100**2
for y in range(LIMIT):
cand = innerphi(y / LIMIT)
m = min(m, cand)
if m == 0:
return 0
return m
RRMY = RealField(200)
def integral1brute(TERMS):
ret = RRMY(0)
for i in range(1, TERMS - 1):
j = RRMY(i + 1)
xi = RRMY(1. * i / TERMS)
xj = RRMY(1. * j / TERMS)
mid = RRMY((xi + xj)/2.)
height = RRMY(brutephi(mid))
width = RRMY(psi(xj) - psi(xi))
ret += RRMY(n(height * width))
return ret
def integrand2(x):
x = RRMY(x)
return brutephi(x) / RRMY(x^2)

def integral2(TERMS):
ret = RRMY(0)
width = RRMY(((1./ms[q-r]) - 0) / TERMS / 2.)
# 0th is 0, so we ignore it
for i in range(1, TERMS + 1):
x = RRMY(0 + (1 - (1/ms[q-r])) * i / TERMS)
height = integrand2(x)
ret += RRMY(n(2 * height))
ret += integrand2(1./ms[q-r])
ret *= RRMY(width)
return ret
def compC1(TERMS):
ft = first_term
i1 = integral1brute(TERMS)
i2 = integral2(TERMS)
print(ft, i1, i2, "C1 is approx", ft - i1 + i2)

compC1(10000)
# 403
# 226.88252
# 6.9549373
# C1 is approx 183.07241


But I am getting the output as: $$227.5801...$$ which is the value of $$C_0$$.

• I let my code run for a very long time (brute force power). I'm going to guess that you just haven't let it terminate yet. Commented Mar 30 at 11:17
• @davidlowryduda Orry, I could not get you. I let the above code run but it is only showing $227.58...$ as an answer.
– Max
Commented Mar 30 at 11:19

Not an answer just some speculation about gamma function and Bernstein's polynomial which is of independent interest here but could help .

Using fallaciously the inverse function of $$f(x)$$:

$$f\left(x\right)=e^{x^{4}-\ln\left(\frac{1}{y}\right)}!,x!=\Gamma(x+1),0

Then using Bernstein form it seems we have $$\forall x>0$$:

$$\lim_{n\to \infty}\sum_{k=0}^{n}\frac{f\left(\frac{k}{n}\right)n!}{k!\left(n-k\right)!}\left(\ln\left(x^{\frac{1}{6}}+1\right)\right)^{k}\left(1-\ln\left(1+x^{\frac{1}{6}}\right)\right)^{\left(n-k\right)}=y!$$

The only advantage I see we have a local inversion of the Gamma function and so the digamma function .

Perhaps I'm wrong .

• (+$1$) Thanks. Do you have a sagemath code for my question? It would be highly appreciated
– Max
Commented Mar 30 at 11:14
• @Max thanks ! I'm really not good for code but perhaps irrelevant I find github.com/paulmasson/sagemath-docs go to psi function . Commented Mar 30 at 12:17