Prime numbers miscellaneous
Revision as of 06:35, 7 June 2011
For a context to this, please see Prime numbers.
1 Implicit Heap
The following is an original implicit heap implementation for the sieve of
Eratosthenes, kept here for historical record. The Prime_numbers#Tree Merging with Wheel section simplifies it, removing the
People a structure altogether, and improves upon it by using a folding tree structure better adjusted for primes processing, and a wheel optimization.
See also the message threads Re: "no-coding" functional data structures via lazyness for more about how merging ordered lists amounts to creating an implicit heap and Re: Code and Perf. Data for Prime Finders for an explanation of the
People a structure that makes it work.
data People a = VIP a (People a) | Crowd [a] mergeP :: Ord a => People a -> People a -> People a mergeP (VIP x xt) ys = VIP x $ mergeP xt ys mergeP (Crowd xs) (Crowd ys) = Crowd $ merge xs ys mergeP xs@(Crowd (x:xt)) ys@(VIP y yt) = case compare x y of LT -> VIP x $ mergeP (Crowd xt) ys EQ -> VIP x $ mergeP (Crowd xt) yt GT -> VIP y $ mergeP xs yt merge :: Ord a => [a] -> [a] -> [a] merge xs@(x:xt) ys@(y:yt) = case compare x y of LT -> x : merge xt ys EQ -> x : merge xt yt GT -> y : merge xs yt diff xs@(x:xt) ys@(y:yt) = case compare x y of LT -> x : diff xt ys EQ -> diff xt yt GT -> diff xs yt foldTree :: (a -> a -> a) -> [a] -> a foldTree f ~(x:xs) = x `f` foldTree f (pairs xs) where pairs ~(x: ~(y:ys)) = f x y : pairs ys primes, nonprimes :: [Integer] primes = 2:3:diff [5,7..] nonprimes nonprimes = serve . foldTree mergeP . map multiples $ tail primes where multiples p = vip [p*p,p*p+2*p..] vip (x:xs) = VIP x $ Crowd xs serve (VIP x xs) = x:serve xs serve (Crowd xs) = xs
nonprimes effectively implements a heap, exploiting lazy evaluation.
2 Prime Wheels
The idea of only testing odd numbers can be extended further. For instance, it is a useful fact that every prime number other than 2 and 3 must be of the form 6k + 1 or 6k + 5. Thus, we only need to test these numbers:
primes :: [Integer] primes = 2:3:primes' where 1:p:candidates = [6*k+r | k <- [0..], r <- [1,5]] primes' = p : filter isPrime candidates isPrime n = all (not . divides n) $ takeWhile (\p -> p*p <= n) primes' divides n p = n `mod` p == 0
Such a scheme to generate candidate numbers first that avoid a given set of primes as divisors is called a prime wheel. Imagine that you had a wheel of circumference 6 to be rolled along the number line. With spikes positioned 1 and 5 units around the circumference, rolling the wheel will prick holes exactly in those positions on the line whose numbers are not divisible by 2 and 3.
A wheel can be represented by its circumference and the spiked positions.
data Wheel = Wheel Integer [Integer]
We prick out numbers by rolling the wheel.
roll (Wheel n rs) = [n*k+r | k <- [0..], r <- rs]
The smallest wheel is the unit wheel with one spike, it will prick out every number.
w0 = Wheel 1 
nextSize (Wheel n rs) p = Wheel (p*n) [r' | k <- [0..(p-1)], r <- rs, let r' = n*k+r, r' `mod` p /= 0]
mkWheel ds = foldl nextSize w0 ds
Now, we can generate prime numbers with a wheel that for instance avoids all multiples of 2, 3, 5 and 7.
primes :: [Integer] primes = small ++ large where 1:p:candidates = roll $ mkWheel small small = [2,3,5,7] large = p : filter isPrime candidates isPrime n = all (not . divides n) $ takeWhile (\p -> p*p <= n) large divides n p = n `mod` p == 0
It's a pretty big wheel with a circumference of 210 and allows us to calculate the first 10000 primes in convenient time.
3 Using IntSet for a traditional sieve
module Sieve where import qualified Data.IntSet as I -- findNext - finds the next member of an IntSet. findNext c is | I.member c is = c | c > I.findMax is = error "Ooops. No next number in set." | otherwise = findNext (c+1) is -- mark - delete all multiples of n from n*n to the end of the set mark n is = is I.\\ (I.fromAscList (takeWhile (<=end) (map (n*) [n..]))) where end = I.findMax is -- primes - gives all primes up to n primes n = worker 2 (I.fromAscList [2..n]) where worker x is | (x*x) > n = is | otherwise = worker (findNext (x+1) is) (mark x is)
(doesn't look like it runs very efficiently).