Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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[RiemannR[1000]]
N[Sum[MoebiusMu[n]*LogIntegral[1000^(1/n)]/n, {n, 1, Infinity}]]

in Mathematica, I get different numerical values:

168.359
168.341

What am I doing wrong?

share|improve this question
1  
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. –  KennyColnago Jun 7 '13 at 2:34
    
@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.) –  J. M. Jun 7 '13 at 3:04

1 Answer 1

up vote 7 down vote accepted

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"]]
   168.35915686601484

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

% - RiemannR[1000]
   -0.000289415

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[].)

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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