-- |
-- 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 :: a
approxPrimeCountOverestimateLimit = a
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 :: a -> a
approxPrimeCount = Double -> a
forall a b. (RealFrac a, Integral b) => a -> b
truncate (Double -> a) -> (a -> Double) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double -> Double) -> (a -> Double) -> a -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
appi (Double -> Double) -> (a -> Double) -> a -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral

-- | Following property holds:
--
-- > nthPrimeApprox n <= nthPrime n || n >= nthPrimeApproxUnderestimateLimit
nthPrimeApproxUnderestimateLimit :: Integer
nthPrimeApproxUnderestimateLimit :: Integer
nthPrimeApproxUnderestimateLimit = Integer
1000000000000

-- | @'nthPrimeApprox' n@ gives an
--   approximation to the n-th prime. The approximation
--   is fairly good for @n@ large enough.
nthPrimeApprox :: Integer -> Integer
nthPrimeApprox :: Integer -> Integer
nthPrimeApprox = Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
1 (Integer -> Integer) -> (Integer -> Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
truncate (Double -> Integer) -> (Integer -> Double) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
nthApp (Double -> Double) -> (Integer -> Double) -> Integer -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Double) -> (Integer -> Integer) -> Integer -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
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 :: Double -> Double
appi Double
x = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
yDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
300000 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
7Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ll
  where
    y :: Double
y = Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h))
    w :: Double
w = Double -> Double
forall a. Floating a => a -> a
log Double
x
    l :: Double
l = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
w
    ll :: Double
ll = Double -> Double
forall a. Floating a => a -> a
log Double
w
    h :: Double
h | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
10000000    = Double
2.5625
      | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
50000000    = Double
2.5
      | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
120000000   = Double
617Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
256
      | Bool
otherwise       = Double
2.0625 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
13.25Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
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 :: Double -> Double
nthApp Double
x = Double
a
  where
    l :: Double
l  = Double -> Double
forall a. Floating a => a -> a
log Double
x
    ll :: Double
ll = Double -> Double
forall a. Floating a => a -> a
log Double
l
    li :: Double
li = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
l
    l2 :: Double
l2 = Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ll
    a :: Double
a  = Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
liDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
liDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
0.3Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
liDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
0.02970812Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
l2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
li))
                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
8.725Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
2.749)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
llDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
3.892)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
li))) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
lDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ll Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
35