Verify the Riemann Hypothesis for first 1000 zeros.

Question

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.

Answer 1

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)