69 lines
2.5 KiB
Haskell
69 lines
2.5 KiB
Haskell
{-# LANGUAGE PostfixOperators #-}
|
||
{-# LANGUAGE DataKinds #-}
|
||
{-# LANGUAGE FlexibleInstances #-}
|
||
|
||
import Data.Numbers.Primes
|
||
import Data.Array.Unboxed hiding ((!))
|
||
import qualified Data.Array.Unboxed as Array
|
||
import Data.CReal
|
||
import Data.CReal.Internal
|
||
import GHC.TypeLits
|
||
|
||
instance KnownNat n => Enum (CReal n) where
|
||
toEnum i = fromIntegral i
|
||
fromEnum _ = error "Cannot fromEnum CReal"
|
||
enumFrom = iterate (+ 1)
|
||
enumFromTo n e = takeWhile (<= e) $ iterate (+ 1)n
|
||
enumFromThen n m = iterate (+(m-n)) n
|
||
enumFromThenTo n m e = if m >= n then takeWhile (<= e) $ iterate (+(m-n)) n
|
||
else takeWhile (>= e) $ iterate (+(m-n)) n
|
||
|
||
|
||
-- partial_sum x y a b = (p,q) where
|
||
-- p/q = sum_{a<i<=b} x(i) / poduct_{a<j<=j} y(j)
|
||
-- The complexity of partial_sum x y 0 n is O(n log n)
|
||
partial_sum x y = pq where
|
||
pq a b = if a>=b then (0,1)
|
||
else if a==b-1 then (fromIntegral $ x b, fromIntegral $ y b )
|
||
else (p_ab,q_ab)
|
||
where
|
||
c=(a+b) `div` 2
|
||
(p_ac,q_ac) = pq a c
|
||
(p_cb,q_cb) = pq c b
|
||
p_ab = p_cb + q_cb*p_ac
|
||
q_ab = q_ac*q_cb
|
||
|
||
-- c is the real constant that is used in the formula for primes
|
||
-- c = sum_{1<i} p_i / (2i+1)!
|
||
-- where p_i is i-th prime.
|
||
-- This will work for any sequence of integers p, where |p_n| < 2n(2n+1) * 0.375
|
||
c = crMemoize f where
|
||
f n = 2^n * p `div` q where
|
||
n' = fromIntegral n
|
||
u = head [ceiling (x) | x<-[(n' * log 2/ (log n'-1)/2 ) ..] , 2*x*log (2*x) - 2*x > n'*log 2]
|
||
-- Invariant: (2u+1)! > 2^n
|
||
ar :: UArray Int Int
|
||
ar = listArray (1,u) $ primes
|
||
(p,q) = partial_sum (ar Array.!) (\n-> 2*n*(2*n+1) ) 0 u
|
||
|
||
|
||
-- Fractorial part of x
|
||
-- By definition it is in the interval [-0.5; 0.5]
|
||
-- But it gurantes to work corectly if fractional part of x is in (-0.375; 0.375)
|
||
fract x = x - fromIntegral (round (x :: CReal 3))
|
||
|
||
-- Factorial.
|
||
-- The complexity of (n!) is O(n log n) (which is better than O(n^2) for product [1..n] )
|
||
(!) :: (RealFrac a, Num b) => a -> b
|
||
(!) = fromIntegral . snd . partial_sum (const 0) id 0 . round
|
||
|
||
-- Analytic function for n-th prime.
|
||
-- NB. Strictly speaking this function is not analytic, because it uses factorial, fractional part and round functions
|
||
-- To make it truly analytic you need to replace
|
||
-- fract x = acos (cos (2*pi*x)) / (2*pi)
|
||
-- round x = x - fract x
|
||
-- and use the Gamma function instead of factorial.
|
||
-- Then you will get analytic function prime :: CReal 0 -> CReal 0
|
||
|
||
prime n = round( 2*n*(2*n+1) * fract ( c * ((2*n-1)!)))
|