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
use ArbitraryPrecision.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 (which gives an approximation to the nthPrime
function) 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
n = 455052511
secantInvert[getFunction["primeCount", 1], n, 1/2 n ln[n], 2 n ln[n]]
TODO: Improve the hard-coded bounds passed to the secantInvert function
above, maybe by using one of the bounds given here:
https://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number
which states that by Rosser's theorem, nthPrime[n] > n ln[n] so the lower
bound above can be improved.
https://en.wikipedia.org/wiki/Prime_number_theorem#Prime-counting_function_in_terms_of_the_logarithmic_integral
https://en.wikipedia.org/wiki/Prime_number_theorem#Non-asymptotic_bounds_on_the_prime-counting_function
*/
primeCount[x] :=
{
sum = 1
k = 1
term = 0
lnx = arbitraryLn[x]
lnxterm = 1
lastsum = sum
do
{
lnxterm = lnxterm * lnx
term = lnxterm / RiemannPrimeDenomCache.get[k]
lastsum = sum
sum = sum + term
//println["$k\t$sum"]
k = k + 1
} while lastsum != sum and abs[term] > 1e-20
return sum
}
/** This class caches values used in the primeCount function. */
class RiemannPrimeDenomCache
{
class var cache = new dict
class get[k] :=
{
val = cache@k
if val == undef
val = cache@k = k k! RiemannZeta.zeta[k+1] // Calc and store into cache
return val
}
}
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, eliasen@mindspring.com