content: -- |
  -- Module:      Math.NumberTheory.Primes.Counting.Approximate
  -- Copyright:   (c) 2011 Daniel Fischer
  -- Licence:     MIT
  -- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
  --
  -- Approximations to the number of primes below a limit and the
  -- n-th prime.
  --
  module Math.NumberTheory.Primes.Counting.Approximate
      ( approxPrimeCount
      , approxPrimeCountOverestimateLimit
      , nthPrimeApprox
      , nthPrimeApproxUnderestimateLimit
      ) where
  
  -- For prime p = 3742914359 we have
  --   approxPrimeCount p = 178317879
  --         primeCount p = 178317880
  
  -- | Following property holds:
  --
  -- > approxPrimeCount n >= primeCount n || n >= approxPrimeCountOverestimateLimit
  approxPrimeCountOverestimateLimit :: Integral a => a
  approxPrimeCountOverestimateLimit = 3742914359
  
  -- | @'approxPrimeCount' n@ gives an
  --   approximation of the number of primes not exceeding
  --   @n@. The approximation is fairly good for @n@ large enough.
  approxPrimeCount :: Integral a => a -> a
  approxPrimeCount = truncate . max 0 . appi . fromIntegral
  
  -- | Following property holds:
  --
  -- > nthPrimeApprox n <= nthPrime n || n >= nthPrimeApproxUnderestimateLimit
  nthPrimeApproxUnderestimateLimit :: Integer
  nthPrimeApproxUnderestimateLimit = 1000000000000
  
  -- | @'nthPrimeApprox' n@ gives an
  --   approximation to the n-th prime. The approximation
  --   is fairly good for @n@ large enough.
  nthPrimeApprox :: Integer -> Integer
  nthPrimeApprox = max 1 . truncate . nthApp . fromIntegral . max 3
  
  -- Basically the approximation of the prime count by Li(x),
  -- adjusted to give close but slightly too high estimates
  -- in the interesting range. The constants are empirically
  -- determined.
  appi :: Double -> Double
  appi x = y - y/300000 + 7*ll
    where
      y = x*l*(1+l*(1+l*h))
      w = log x
      l = 1/w
      ll = log w
      h | x < 10000000    = 2.5625
        | x < 50000000    = 2.5
        | x < 120000000   = 617/256
        | otherwise       = 2.0625 + l*(3+ll*l*(13.25+ll*l*57.75))
  
  -- Basically an approximation to the inverse of Li(x), with
  -- empirically determined constants to get close results
  -- in the interesting range.
  nthApp :: Double -> Double
  nthApp x = a
    where
      l  = log x
      ll = log l
      li = 1/l
      l2 = ll*ll
      a  = x*(l+ll-1+li*(ll-2-li*(ll*(0.3+li*(1+0.02970812*l2*l2*l2*li))
                  + 8.725*(ll-2.749)*(ll-3.892)*li))) + l*ll + 35
  