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 18354 days, 1 hours, 10 minutes ago.