/** This attempts to implement the Riemann Prime Counting function.
http://mathworld.wolfram.com/RiemannPrimeCountingFunction.html
which is somewhat more accurate than the LogIntegral function, which
is asymptotically really accurate, see LogIntegral.frink
This uses the simpler Gram's formula (eq. 12 in the link above.)
TODO: Perform ln[x] and RiemannZeta[x] to arbitrary precision? Would it
improve anything?
*/
use RiemannZeta.frink
/** Return the Riemann Prime Counting Function's estimate of the number of
prime numbers less than or equal to x.
You can compare these estimates to the exact counts tabulated at:
http://mathworld.wolfram.com/PrimeCountingFunction.html
You can invert this count using something like the secant inversion
function. The following should return a number around 10^10, because
there are exactly 455052511 primes less than 10^10:
use secant.frink
use RiemannPrime.frink
secantInvert[getFunction["primeCount", 1], 455052511, 1, quadrillion]
*/
primeCount[x] :=
{
sum = 1
k = 1
term = 0
lnx = ln[x]
lnxterm = 1
lastsum = sum
do
{
lnxterm = lnxterm * lnx
term = lnxterm / (k k! RiemannZeta[k+1])
lastsum = sum
sum = sum + term
// println["$k\t$sum"]
k = k + 1
} while lastsum != sum and abs[term] > 1e-20
return sum
}