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