1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
-- |
-- 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
|
