    -- |
    -- 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 = Double -> a
    forall b. Integral b => Double -> b
    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 b. Integral b => Double -> b
    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