# How can I program the RiemannR function using the LogIntegral command?

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":

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


in Mathematica, I get different numerical values:

168.359
168.341


What am I doing wrong?

-
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.) –  Ｊ. Ｍ. 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},
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 NumberTheoryNumberTheoryFunctionsDumpPlainRiSeries[].)

-