/** 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 }