Download or view RiemannPrime.frink in plain text format
/** 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">secant.frink
use RiemannPrime.frink">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
}
Download or view RiemannPrime.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19963 days, 21 hours, 39 minutes ago.