# Plotting Riemann Prime Counting

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!!

• Instead of "Real" in the definition of F[x], use "Re", and I think replace Log[x] in the Plot command with Log[t]. The Plot limits should be, say, 2 to 10. The function blows up at 1. Commented Feb 1, 2021 at 18:40
• RiemannR[] is now built-in. Commented Feb 1, 2021 at 20:24
• @KennyColnago and please read my edit. The graph won't change. Commented Feb 1, 2021 at 23:37

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) *
(
- 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}]]


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}]


• Making the replacement I suggested in the comments: Plot[{PrimePi[t], RiemannR[t] - M[t] - 1/Log[t] + ArcTan[Pi/Log[t]]/Pi}, {t, 0, 10}] Commented Feb 1, 2021 at 20:32
• Is it possible to remove the yellow graph in the intervall [0,1.75]? And thank you so much guys! You literally saved my life, I have a presentation tomorrow, phew... Commented Feb 1, 2021 at 20:44
• @vitamin, use CondtionalExpression[(* stuff *), ! (0 < x < 7/4)] for that. Commented Feb 2, 2021 at 2:04
• @J.M. okay thanks. Did you read my other comment above? If you don't want to answer it, no problem. Can ask a new question. Commented Feb 2, 2021 at 12:39
• @vitamin, that should definitely be a new question (but you can link to this one for context). Commented Feb 2, 2021 at 12:40