1
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Feb 1 '21 at 18:40
  • 1
    $\begingroup$ RiemannR[] is now built-in. $\endgroup$ Feb 1 '21 at 20:24
  • $\begingroup$ @KennyColnago and please read my edit. The graph won't change. $\endgroup$ Feb 1 '21 at 23:37
4
$\begingroup$

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

poor approximation

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

better approximation

$\endgroup$
2
$\begingroup$

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

enter image description here

$\endgroup$
5
  • $\begingroup$ 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}] $\endgroup$ Feb 1 '21 at 20:32
  • $\begingroup$ 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... $\endgroup$ Feb 1 '21 at 20:44
  • $\begingroup$ @vitamin, use CondtionalExpression[(* stuff *), ! (0 < x < 7/4)] for that. $\endgroup$ Feb 2 '21 at 2:04
  • $\begingroup$ @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. $\endgroup$ Feb 2 '21 at 12:39
  • $\begingroup$ @vitamin, that should definitely be a new question (but you can link to this one for context). $\endgroup$ Feb 2 '21 at 12:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.