Programming Explicit Formula for Primes

Question

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)/n*f(N**(1/n)

        return suma

Thank you very much.

Answer 1

Most problems with such evaluations are detailed here.

Your first problem is probably the use of the logarithmic integral and you'll have to use instead : $$\operatorname{Li}(x^\rho)=\operatorname{Ei}(\rho\,\log\;x)$$ with $\operatorname{Ei}$ the exponential integral (the correct phase of $\rho$ was lost!).

See too this answer for visual results.