I would like to program the RiemannR function using the LogIntegral command because I would like to later experiment with a double integral instead of the LogIntegral function.

I have tried to program the RiemannR function starting with the LogIntegral command as defined in the Stein-Mazur paper called "Primes":

RiemannR function as in Stein Mazur paper on primes

That agrees almost completely with the definition given in the Mathematica help page about the RiemannR command:

$$R(x)=\sum_n^\infty \mu(n)\mathrm{li}(x^{1/n})/n$$

The only differences I see are that the summation limit $n=1$ is given in the Stein-Mazur paper while in the Mathematica help page it is not. Also the Stein-Mazur paper uses $\mathrm{Li}$ for the logarithmic integral while the Mathematica help page uses $\mathrm{li}$. However in the Mathematica help page for the LogIntegral command the lower case $\mathrm{li}$ is used.

My problem is that when I enter the following lines:

N[Sum[MoebiusMu[n]*LogIntegral[1000^(1/n)]/n, {n, 1, Infinity}]]

in Mathematica, I get different numerical values:


What am I doing wrong?

  • 2
    $\begingroup$ If you are investigating the prime number theorem, please note the chapter on the zeta function in Wagon's Mathematica In Action, where he refers to a book by H.M. Edwards called Riemann's Zeta Function. Wagon cautions that ExpIntegralEi[r Log[x]] must be used instead of LogIntegral[x^r] when x is real and r is complex. This is because the complex log function returns r Log[x]-2 Pi n i, with n chosen to make the imaginary part lie between -Pi and Pi, which leads to incorrect sums. $\endgroup$ Jun 7 '13 at 2:34
  • $\begingroup$ @Kenny, yes, that's another good reason to use ExpIntegralEi[], at least in this context. (I had forgotten to point out that Stan Wagon indeed discussed Riemann's function in his book.) $\endgroup$
    – J. M.'s torpor
    Jun 7 '13 at 3:04

I've found that NSum[] takes a bit too long here to compute Riemann's prime-counting function, so I've resorted to generating the terms and summing them:

With[{n = 8*^3},
     Total[MoebiusMu[Range[n]] N[LogIntegral[1000^(1/Range[n])]/Range[n]],
           Method -> "CompensatedSummation"]]

which gives a result close to that of RiemannR[]:

% - RiemannR[1000]

It is, however, faster (at least in some of the tests I did) to exploit an identity connecting the logarithmic integral and the exponential integral $\mathrm{Ei}(z)$:

With[{n = 8*^3}, 
     Total[MoebiusMu[Range[n]] N[ExpIntegralEi[Log[1000]/Range[n]]/Range[n]],
           Method -> "CompensatedSummation"]]

But, the best method to compute Riemann's function would be to use Gram's series:

$$R(x)=1+\sum_{k=1}^\infty \frac{(\log\,x)^k}{k k!\zeta(k+1)}$$

In Mathematica:

1 + NSum[(Log[1000]^k)/(k k! Zeta[k + 1]), {k, 1, ∞}]

As it turns out, if you inspect how Mathematica internally implements RiemannR[], you'll find that Gram's series is indeed used for large arguments. (In particular, look at NumberTheory`NumberTheoryFunctionsDump`PlainRiSeries[].)


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.