-- |
-- Module: Math.NumberTheory.Primes.Counting.Approximate
-- Copyright: (c) 2011 Daniel Fischer
-- Licence: MIT
-- Maintainer: Daniel Fischer
-- Stability: Provisional
-- Portability: portable
--
-- Approximations to the number of primes below a limit and the
-- n-th prime.
--
{-# OPTIONS_HADDOCK hide #-}
module Math.NumberTheory.Primes.Counting.Approximate
( approxPrimeCount
, nthPrimeApprox
) where
-- | @'approxPrimeCount' n@ gives (for @n > 0@) an
-- approximation of the number of primes not exceeding
-- @n@. The approximation is fairly good for @n@ large enough.
-- The number of primes should be slightly overestimated
-- (so it is suitable for allocation of storage) and is
-- never underestimated for @n <= 10^12@.
approxPrimeCount :: Integral a => a -> a
approxPrimeCount = truncate . appi . fromIntegral
-- | @'nthPrimeApprox' n@ gives (for @n > 0@) an
-- approximation to the n-th prime. The approximation
-- is fairly good for @n@ large enough. Dual to
-- @'approxPrimeCount'@, this estimate should err
-- on the low side (and does for @n < 10^12@).
nthPrimeApprox :: Integral a => a -> a
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