-- |
-- 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 :: forall a. Integral a => 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 :: forall a. Integral a => a -> a
approxPrimeCount = forall a b. (RealFrac a, Integral b) => a -> b
truncate forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Ord a => a -> a -> a
max Double
0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
appi forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = forall a. Ord a => a -> a -> a
max Integer
1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (RealFrac a, Integral b) => a -> b
truncate forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
nthApp forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 forall a. Num a => a -> a -> a
- Double
yforall a. Fractional a => a -> a -> a
/Double
300000 forall a. Num a => a -> a -> a
+ Double
7forall a. Num a => a -> a -> a
*Double
ll
  where
    y :: Double
y = Double
xforall a. Num a => a -> a -> a
*Double
lforall a. Num a => a -> a -> a
*(Double
1forall a. Num a => a -> a -> a
+Double
lforall a. Num a => a -> a -> a
*(Double
1forall a. Num a => a -> a -> a
+Double
lforall a. Num a => a -> a -> a
*Double
h))
    w :: Double
w = forall a. Floating a => a -> a
log Double
x
    l :: Double
l = Double
1forall a. Fractional a => a -> a -> a
/Double
w
    ll :: Double
ll = forall a. Floating a => a -> a
log Double
w
    h :: Double
h | Double
x forall a. Ord a => a -> a -> Bool
< Double
10000000    = Double
2.5625
      | Double
x forall a. Ord a => a -> a -> Bool
< Double
50000000    = Double
2.5
      | Double
x forall a. Ord a => a -> a -> Bool
< Double
120000000   = Double
617forall a. Fractional a => a -> a -> a
/Double
256
      | Bool
otherwise       = Double
2.0625 forall a. Num a => a -> a -> a
+ Double
lforall a. Num a => a -> a -> a
*(Double
3forall a. Num a => a -> a -> a
+Double
llforall a. Num a => a -> a -> a
*Double
lforall a. Num a => a -> a -> a
*(Double
13.25forall a. Num a => a -> a -> a
+Double
llforall a. Num a => a -> a -> a
*Double
lforall 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  = forall a. Floating a => a -> a
log Double
x
    ll :: Double
ll = forall a. Floating a => a -> a
log Double
l
    li :: Double
li = Double
1forall a. Fractional a => a -> a -> a
/Double
l
    l2 :: Double
l2 = Double
llforall a. Num a => a -> a -> a
*Double
ll
    a :: Double
a  = Double
xforall a. Num a => a -> a -> a
*(Double
lforall a. Num a => a -> a -> a
+Double
llforall a. Num a => a -> a -> a
-Double
1forall a. Num a => a -> a -> a
+Double
liforall a. Num a => a -> a -> a
*(Double
llforall a. Num a => a -> a -> a
-Double
2forall a. Num a => a -> a -> a
-Double
liforall a. Num a => a -> a -> a
*(Double
llforall a. Num a => a -> a -> a
*(Double
0.3forall a. Num a => a -> a -> a
+Double
liforall a. Num a => a -> a -> a
*(Double
1forall a. Num a => a -> a -> a
+Double
0.02970812forall a. Num a => a -> a -> a
*Double
l2forall a. Num a => a -> a -> a
*Double
l2forall a. Num a => a -> a -> a
*Double
l2forall a. Num a => a -> a -> a
*Double
li))
                forall a. Num a => a -> a -> a
+ Double
8.725forall a. Num a => a -> a -> a
*(Double
llforall a. Num a => a -> a -> a
-Double
2.749)forall a. Num a => a -> a -> a
*(Double
llforall a. Num a => a -> a -> a
-Double
3.892)forall a. Num a => a -> a -> a
*Double
li))) forall a. Num a => a -> a -> a
+ Double
lforall a. Num a => a -> a -> a
*Double
ll forall a. Num a => a -> a -> a
+ Double
35