Riemann Zeta Function – Number of Zeros

riemann-zetazeta-functions

I want to write a program that calculates the number of zeros (It is not necessary to identify them, just the number of them) between 0 and x for the Riemann Zeta function, being x the imaginary part of z: 1/2 + ix

Is there an algorithm out there for doing what I need?

I modelled the Riemann-siegel function because I saw that I need it for doing what I want, but I don't know how to go on.

Sorry for my poor english :/

Best Answer

You may appreciate the following '$\zeta$ zeros counting function' if $t$ is the imaginary part : $$f(t) = \frac{1}{\pi} \Im\left(\ln(\Gamma\left(\frac 14 + \frac{i\,t}2\right)\right) - \frac{t}{2\pi}\ln(\pi) + \frac{1}{\pi} \Im\left(\ln\left(\zeta\left(\frac 12 + i\,t\right)\right)\right) + 1$$ Explication : if we multiply this by $\pi$ we note that the first imaginary term combined to the next term returns the actual argument of $\zeta$ (this is the well known Riemann–Siegel theta function) that grows with high regularity while the imaginary term at the right gives the argument of $\zeta$ modulo $2\pi$ which will change sign for every multiple of $\pi$ (that is at every zero of $\zeta$ !).

More about the way this formula was obtained here.

Illustration using the pari/gp code at the end (zoom possible) :

zeros counting function

The first step is between $14.1$ and $14.2$, the second between $21$ and $21.1$ and so on so that we are really counting the zeros of $\zeta$...

The picture is not always so nice and, from time to time (starting after $t=415$ I think), there will be a non-propagating error of $\,\pm 2\,$ when $\,\displaystyle\frac 1{\pi}\rm{Arg}\;\zeta(1/2+it)\;$ crosses the upper bound $+1$ or lower bound $-1$ (so that this appears rather when two consecutive zeros are distant i.e. for large loops as displayed in the other answer).

This formula was in this MO thread with further references to Guinand's article : 'A summation formula in the theory of prime numbers'.

pari/gp code used for the picture :

f(t)=imag(lngamma(1/4+I*t/2))/Pi-t/(2*Pi)*log(Pi)+imag(log(zeta(1/2+I*t)))/Pi+1