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

RiemannR[]
is now built-in. $\endgroup$