RiemannPrime.frink

View or download 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
    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
}


View or download 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 18233 days, 19 hours, 13 minutes ago.