# Verify the Riemann Hypothesis for first 1000 zeros.

I'm trying to verify the Riemann hypothesis for the first 1000 zeros using the Euler-MacLaurin expansion of $\zeta(s)$.

But here's where my problems begin, firstly my python code doesn't seem to be working, I'm essentially using the method outlined in Mike Rubinstein - "Computational Methods and Experiments in Analytic Number Theory" [https://arxiv.org/pdf/math/0412181v1.pdf] and here is the code I have:

def CalculateRZ(s, tol):
A = int(0)
while (2*A - 1 < tol + 0.5 * math.log10(abs(s + 2*A - 2))):
A = A + 1
B = int(2 * A)
N = int(math.ceil((10 / 2 * pi) * abs(s + B - 2)))

sum1 = 0.0
for n in range(1, N + 1):
sum1 = sum1 + n ** (-s)

ksum = 0.0
for k in range(1, B + 1):
ksum = ksum + Binomial(s + k - 2, k - 1) * (bernoulli(k) / k) * N ** (-s - k - 1)

sum2 = (N ** (1 - s)) / (s - 1) + ksum

EulerApx = sum1 + sum2

return EulerApx


where Binomial(n, k), and bernoulli(k) are exactly what they look like.

This seems perfectly correct to me, however when I set

x = 1/2 + 14.134725141734693790457251983562470270784257115699243175685567460149 * 1j


(i.e., the first zero) and run

CalculateRZ(x, 100)


I get as an output $\zeta(x) = -0.008282095160420333 + 0.009227302793957781j$ which may be kind of close to zero but certainly not to 100 digits of accuracy and likely not enough accuracy to verify the Riemann Hypothesis.

Further to this, once I have this code working I'm unsure how to go about verifying the Riemann Hypothesis afterwards.

Would I be better off creating a new function for $Z(t)$ then using Newton Raphson to count the zeros? If so how would I go about computing $Z(t)$? Should I just make it as a function return $\zeta(1/2+1t)e^{i\theta(t)}$ using my version of zeta and a function that just straight up calculates $\theta(t)$? Where should I set my start points for N-R method to ensure I hit every zero once and only once?

Once the zeros are counted I understand I can then integrate $\zeta$ over some rectangle, and divide that number by $2\pi$ which will give me the number of singularities which should match up to the number of zeros found using Newton Raphson, therefore showing that all zeros lie on the line $Z(t)$, and therefore verifying the RH. But how would this integral be calculated?

I know this is a huuge question so any help would be appreciated even if it's not an answer to the full question.

• 10 / 2 * pi is wrong (should be 10/(2*pi)), but that can't be what is causing your problem – Matthew Towers Jan 28 '18 at 13:47
• The exponent of N in ksum is wrong too, it should be -s-k+1. That might be it. – Matthew Towers Jan 28 '18 at 13:54
• Right you are, I now get -0.07756554873082433 - 0.006647684611503335j – kingee Jan 28 '18 at 14:12
• Unfortunately, Newton-Raphson gives you no kind of guarantee that you found all the roots. – Yves Daoust Jan 28 '18 at 16:19
• Damn... know any way I can gurantee to find all roots? Or do you think running it and getting as many as I can would be sufficient to obtain the result I need? – kingee Jan 28 '18 at 17:04

The biggest problem here is that $1/2$ in x = 1/2 + 14.1347 * 1j gets interpreted as integer division so you are trying to compute $\zeta(x)$ for x = 14.1347 * 1j. Write it as 0.5 instead.

Then you also need to fix the typo mentioned in the comments $N^{-s-k-1} \to N^{-s-k+1}$. Running your code with these fixes gives -1.17435267217e-14 - 3.36278987909e-16j which is pretty close to zero (in double precision numbers, which is what you are working with).

Also note that since you are not using arbitrary precision numbers in the calculations there is really no point of running with tol greater than ~$20$.

For completeness here is a fully working code

import math
from sympy import bernoulli

# Binomial coefficients for arbitrary n and integer k
def Binomial(n,k):
res = 1.0
for i in range(0,k):
res *= (n-i)/(k-i)
return res

def CalculateRZ(s, tol):
A = int(0)
while (2*A - 1 < tol + 0.5 * math.log10(abs(s + 2*A - 2))):
A = A + 1
B = int(2 * A)
N = int(math.ceil((10 / (2 * math.pi)) * abs(s + B - 2)))

sum1 = 0.0
for n in range(1, N + 1):
sum1 = sum1 + n ** (-s)

ksum = 0.0
for k in range(1, B + 1):
b = Binomial(s + k - 2, k - 1)
ksum = ksum + b * (bernoulli(k) / k) * N ** (-s - k + 1)

sum2 = (N ** (1 - s)) / (s - 1) + ksum

EulerApx = sum1 + sum2

return complex(EulerApx)

s = 0.5 +  14.1347251417j
tol = 20
print CalculateRZ(s, tol)

• Wow that looks great thanks, however when I run that exact code I get the solution as: -0.047483292047607095 - 0.15081557275065474j what does this mean? Something wrong with my python? – kingee Jan 28 '18 at 17:07
• Using python 3.6... – kingee Jan 28 '18 at 17:14
• You know what it may be the binomial function: def bernoulli(n): A = [0] * (n + 1) for m in range(n + 1): A[m] = Fr(1, m + 1) for j in range(m, 0, -1): A[j - 1] = j * (A[j - 1] - A[j]) return A[0] # (which is Bn) I just got this online and have a feeling its not doing what I want it to correctly... – kingee Jan 28 '18 at 17:21