Plotting Riemann Prime Counting

Question

I want to program the formula in this page. Didn't take me long to get to my code below. But 1) it only plots the prime counting function and 2) evaluation takes incredibly long. How can I fix it?

R[x_] := Sum[MoebiusMu[n]/n*LogIntegral[x^(1/n)],{n,100}]
F[x_] := -2*Real[Sum[ExpIntegralEi[ZetaZero[k]*Log[x]],{k,10}]]
M[x_] := Sum[MoebiusMu[n]/n*F[x],{n,100}]
Plot[{PrimePi[t],R[t]-M[t]-1/Log[t]+1/Pi*ArcTan[Pi/Log[x]]},{t,0,10},
PlotRange ->All]

EDIT, major problem: If I am changing the amount of zeros, in other words k, the graph stays the same which is not good. It should slowly become a better and better approximation. Please help!!

Answer 1

Riemann's solution for the prime counting function is

Pi[x]=Sum[(MoebiusMu[n]/n) * J[x^(1/n)], {n,1,nmax}],

where J[x] contains several terms. Define the integral term in J[x].

RiemannJint[x_?NumericQ] := NIntegrate[1/(u (u^2 - 1) Log[u]), {u, x, Infinity}]

Using the first k zeta function zeros, the equation becomes

Sum[(MoebiusMu[n]/n) * 
    (
     LogIntegral[x^(1/n)]
     - Log[2]
     + RiemannJint[x^(1/n)]
     -2 * Total[Re[ExpIntegralEi[N[ZetaZero[Range[k]]]*Log[x]/n]]]
    ), {n, 1, nmax}]

The built-in RiemannR[x] covers the contribution of the LogIntegral term.

The (not-so-good) approximation ArcTan[Pi/Log[x]]/Pi covers the Log[2] and integral (RiemannJint) terms. This approximation is based on nmax=154. Thus,

approxPi[x_, k_Integer] :=
   RiemannR[x] +
   ArcTan[Pi/Log[x]]/Pi -
   2*Sum[(MoebiusMu[n]/n) *
     (Total[Re[ExpIntegralEi[N[ZetaZero[Range[k]]]*Log[x]/n]]]),
     {n, 1, 154}]

For example:

ListLinePlot[
   {Table[{x, PrimePi[x]}, {x, 2., 50., 0.5}],
    ParallelTable[{x, approxPi[x, 10]}, {x, 2., 50., 0.5}],
    ParallelTable[{x, approxPi[x, 100]}, {x, 2., 50., 0.5}]},
   ImageSize -> 500, Frame -> True, GridLines -> Automatic, 
   BaseStyle -> {FontSize -> 14},
   PlotRange -> {{2, 50}, {0, 16}},
   PlotLegends -> 
      Placed[{"True \[Pi][x]", "10 Zeros", "100 Zeros"}, {Left, Top}]]

The approximation ArcTan[Pi/Log[x]]/Pi has little effect for larger x. Removing the approximation improves the fit.

The following code avoids this problematic approximation to Log[2] and the integral.

RiemannJ[x_, k_Integer] := 0 /; x < 2
RiemannJ[x_, k_Integer] := 
   LogIntegral[x] - Log[2.] + RiemannJint[x] - 
   2*Total[Re[ExpIntegralEi[N[ZetaZero[Range[k]]]*Log[x]]]]
SetAttributes[RiemannJ, Listable]

RiemannPi[x_, k_Integer] := Sum[(MoebiusMu[n]/n) RiemannJ[x^(1/n), k], {n, 1, k}]

The corresponding plot compares Pi[x] to RiemannPi[x]. The plot shows near perfect agreement in the sense that RiemannPi[x] converges to Pi[x] when x is not prime, and converges to Pi[x]-1/2 when x is prime. Using more zeros improves the agreement.

ListStepPlot[{
   Table[{x, PrimePi[x]}, {x, 1, 50}],
   ParallelTable[{x, RiemannPi[x, 10]}, {x, 1, 50}],
   ParallelTable[{x, RiemannPi[x, 100]}, {x, 1, 50}]},
   Frame -> True, 
   FrameLabel -> {"Argument x", "PrimePi[x]   and   RiemanPi[x]"},
   PlotLabel -> "RiemannPi[x]", BaseStyle -> {FontSize -> 14},
   PlotRange -> {{0, 50}, {0, 16}}, ImageSize -> 500, 
   GridLines -> Automatic,
   PlotLegends -> 
      Placed[{"True \[Pi][x]", "RiemannPi with 10 Zeros",
             "RiemannPi with 100 Zeros"}, {Left, Top}]]

Answer 2

Using immediate assignments and correcting the mistakes pointed out by @KennyColnago, plotting is quick:

R[x_] = Sum[MoebiusMu[n]/n*LogIntegral[x^(1/n)], {n, 100}];
F[x_] = -2*Re[Sum[ExpIntegralEi[ZetaZero[k]*Log[x]], {k, 10}]];
M[x_] = Sum[MoebiusMu[n]/n*F[x], {n, 100}];
Plot[{PrimePi[t], R[t] - M[t] - 1/Log[t] + 1/Pi*ArcTan[Pi/Log[t]]}, {t, 0, 10}]