# Mathematica code to compute a constant

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)d\psi(x)-\int_0^{1/m_{q-r}}\varphi(x)\frac{dx}{x^2}\right)$$ where $$\psi(x)=\frac{\Gamma'(x)}{\Gamma(x)}$$ is the logarithmic derivative of the Gamma function. 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: I need a code in Mathematica to find the approx value of $$C_1$$.

I tried the following code on Mathematica

r = 3;
q = 13;
eta = {91, 27, 27, 27, 29, 30, 31, 32,
33, 34, 35, 36, 37, 38};

phi[x_] := Sum[Ceiling[y] +
Ceiling[eta[[1]] x - y] - Ceiling[y -
eta[[j]] x] -
Ceiling[(eta[[1]] - eta[[j]]) x - y] -
2*Ceiling[eta[[j]] x], {j, 1, r}] +
Sum[Ceiling[(eta[[1]] - 2
eta[[j]]) x] - Ceiling[y - eta[[j]] x]
-
Ceiling[(eta[[1]] - eta[[j]]) x
- y], {j, r + 1, q}];

psi[x_] := Gamma'[x]/Gamma[x];

C0 = -Re[f0[tau0]];
C1 = r*m1 + Sum[m[j], {j, 2, q - r}]
-
(Integrate[phi[x], {x, 0, 1}] -
Integrate[phi[x]/x^2, {x, 0, 1/m[q -
r]}]);

{C0, C1}


But I am getting output as:

 {-Re[f0[tau0]],3 m1+m[2]+m[3]+m[4]+m[5]+m[6]+m[7]+m[8]+m[9]+m[10]}


Any help would be highly appreciated. Thank you!

• f0[] is undefined and so is m[]. The line breaks look funny. I'm not sure if that's because of how you typed the code in the notebook or how you copied the code into this post. Commented Apr 1 at 11:27
– Max
Commented Apr 1 at 11:35
• You need to add to this script code to calculate the $m_j$ from $\eta_i$. The formula is in the text following displayed equation (2). Actually I am quite curious if you manage to reproduce the constants (I read your MO post) but at the moment I do not have the concentration to complete your code. Commented Apr 1 at 20:19
• Addendum: I just looked more carefully into your code: It contains also other omissions and lacunas - I will try to correct it but it can take a bit of time. Commented Apr 1 at 20:33
• @JürgenBöhm Please correct it. I will be highly obliged. Thank you.
– Max
Commented Apr 2 at 6:10

The following code should compute you the C1 sought after from the paper.

Unfortunately the running times for the NIntegrate commands are very long and I could not compute values in accordance with the result in the paper because with standard options the integrations fail with too big error margins.

But with patience and trying out options like

WorkingPrecision -> 20 (or even 40)


and

Method -> "LocalAdaptive"


you might get it through.

Now the code (I use Floor now instead of Ceiling before, because in the original paper the author writes $$\lfloor\cdot\rfloor$$. The results I now get seem to support this choice)

r = 3;

q = 13;

etalis = {91, 27, 27, 27, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38};

etasl = Table[eta[i] -> etalis[[i + 1]], {i, 0, 13}]

msl = Table[
m[j] -> (Max[eta[r], eta[0] - 2 eta[r + 1],
eta[0] - eta[1] - eta[r + j]]) /. etasl, {j, 1, q - r}]

phi1[x_, y_] :=
(Sum[
(Floor[y] + Floor[eta[0] x - y] - Floor[y - eta[j] x] -
Floor[(eta[0] - eta[j]) x - y] -
2 Floor[eta[j] x]) /. etasl, {j, 1, r}]
+
Sum[
(Floor[(eta[0] - 2 eta[j]) x] - Floor[y - eta[j] x] -
Floor[(eta[0] - eta[j]) x - y]) /. etasl,
{j, r + 1, q}
])

(*
You can try this (terminated PiecewiseExpand after endless running time). Should in
principle vastly improve the computation.

phi11 = PiecewiseExpand[phi1[x, y], {0 <= x <= 1, 0 <= y < 1}]
phi[x_] := Minimize[{phi11[x, y], 0 <= y < 1}, y][[1]]

*)

phi[x_?NumericQ] := Minimize[{phi1[x, y], 0 <= y < 1}, y][[1]]

(*

C11 agrees with the value from the sage computation you elicited
in math.stackexchange.

*)

C11 = (r m[1] + Sum[m[j], {j, 2, q - r}]) /. msl

(*
Try to fiddle around with the options of NIntegrate and
be prepared for long running times and up to 32 GB or more memory
consumption.
*)

C12 = NIntegrate[phi[x] LogGamma''[x], {x, (1/m[q - r]) /. msl, 1},
WorkingPrecision -> 10]

C13 = NIntegrate[
phi[x] LogGamma''[x] - phi[x]/x^2, {x, 0, (1/m[q - r]) /. msl},
WorkingPrecision -> 10]

C1 = C11 - (C12 + C13)

• ($+1$) Thanks. What is the output of $C_1$ that you get?
– Max
Commented Apr 3 at 5:20
• Is this code on Mathematica?
– Max
Commented Apr 3 at 5:26