I am trying to implement an algorithm to compute the following function:
$$\pi_0(x)=\sum_{n=1}^{\infty}\frac{\mu(n)}{n}f(x^{1/n}),$$
where $\mu$ is the Möebius function and
$$f(x)= li(x)-\sum_{\rho}^{\infty}li(x^{\rho})-\ln(2)+\int_{x}^{\infty}\frac{dt}{t(t^2-1)\ln(t)} $$
where $li=\int_{0}^{x}\frac{dt}{\ln(t)}$ is the logrithmic integral function and the sum is over all the non trivial zeros of the Riemann Zeta function. See https://repository.tudelft.nl/islandora/object/uuid:bdce6768-a0f4-463c-9ec7-3ac0a233b51e/datastream/OBJ/download, pages 49-50.
The thing is that i am trying to implement this algorithm since two days ago, and it does not work. I am considering, for example, the first $10$ zeros of the Riemann Zeta function (which I have in the list "zeros"), but it does not converge to the original prime counting function $\pi(x)$. I read that the sum involving the zeros of the Zeta function is conditionally convergent, but I am doing the sum in order, first $1/2+14.13i$, then $1/2-14.13i$, then $1/2+21.02i$ and so on.
For example $\pi_0(100)$ gives me around $-344$ with $100$ roots, which has no sense. What is happening here? Is it because of the conditionally convergence? Or maybe am I losing precision?
This is the code in Python3:
def f(N):
#f function
nroots=100
sumaceros=0
for j in range(1,nroots+1):
sumaceros=sumaceros+li(N**(1/2+1j*(zeros[j][1])))+li(N**(1/2-1j*(zeros[j][1])))
suma=li(N)-sumaceros-np.log(2)+expint(N)
return suma.real
def pi0(N):
#pi0 function
suma=0
for n in range(1,100):
suma=suma+mobius(n)/nf(N*(1/n)
return suma
Thank you very much.
