View source for Prime numbers miscellaneous

For a context to this, please see [[Prime numbers#Implicit_Heap | Prime numbers]].

== Implicit Heap ==

The following is an original implicit heap implementation for the sieve of
Eratosthenes, kept here for historical record. Also, it implements more sophisticated, lazier scheduling. The [[Prime_numbers#Tree merging with Wheel]] section simplifies it, removing the <code>People a</code> structure altogether, and improves upon it by using a folding tree structure better adjusted for primes processing, and a [[Prime_numbers#Euler.27s_Sieve | wheel]] optimization.

See also the message threads [http://thread.gmane.org/gmane.comp.lang.haskell.cafe/25270/focus=25312 Re: "no-coding" functional data structures via lazyness] for more about how merging ordered lists amounts to creating an implicit heap and [http://thread.gmane.org/gmane.comp.lang.haskell.cafe/26426/focus=26493 Re: Code and Perf. Data for Prime Finders] for an explanation of the <code>People a</code> structure that makes it work.

<haskell>
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
</haskell>

<code>nonprimes</code> effectively implements a heap, exploiting lazy evaluation.

== 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 <math>6k+1</math> or <math>6k+5</math>. Thus, we only need to test these numbers:

<haskell>
primes :: [Integer]
primes = 2:3:prs
  where
    1:p:candidates = [6*k+r | k <- [0..], r <- [1,5]]
    prs            = p : filter isPrime candidates
    isPrime n      = all (not . divides n)
                       $ takeWhile (\p -> p*p <= n) prs
    divides n p    = n `mod` p == 0
</haskell>

Here, <hask>prs</hask> is the list of primes greater than 3 and <hask>isPrime</hask> does not test for divisibility by 2 or 3 because the <hask>candidates</hask> by construction don't have these numbers as factors. We also need to exclude 1 from the candidates and mark the next one as prime to start the recursion.

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.
<haskell>
data Wheel = Wheel Integer [Integer]
</haskell>
We prick out numbers by rolling the wheel.
<haskell>
roll (Wheel n rs) = [n*k+r | k <- [0..], r <- rs]
</haskell>
The smallest wheel is the unit wheel with one spike, it will prick out every number.
<haskell>
w0 = Wheel 1 [1]
</haskell>
We can create a larger wheel by rolling a smaller wheel of circumference <hask>n</hask> along a rim of circumference <hask>p*n</hask> while excluding spike positions at multiples of <hask>p</hask>.

<haskell>
nextSize (Wheel n rs) p =
  Wheel (p*n) [r2 | k <- [0..(p-1)], r <- rs,
                    let r2 = n*k+r, r2 `mod` p /= 0]
</haskell>

Combining both, we can make wheels that prick out numbers that avoid a given list <hask>ds</hask> of divisors.
<haskell>
mkWheel ds = foldl nextSize w0 ds
</haskell>

Now, we can generate prime numbers with a wheel that for instance avoids all multiples of 2, 3, 5 and 7.
<haskell>
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
</haskell>
It's a pretty big wheel with a circumference of 210 and allows us to calculate the first 10000 primes in convenient time.

A fixed size wheel is fine, but adapting the wheel size while generating prime numbers quickly becomes impractical, because the circumference grows very fast, as primorial, but the returns quickly diminish, the improvement being just <code>(p-1)/p</code>. See [[Prime_numbers#Euler.27s_Sieve | Euler's Sieve]], or the [[Research papers/Functional pearls|functional pearl]] titled [http://citeseer.ist.psu.edu/runciman97lazy.html Lazy wheel sieves and spirals of primes] for more.

== Using IntSet for a traditional sieve ==
<haskell>

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)
</haskell>

''(doesn't look like it runs very efficiently)''.



== One-liners ==
Produce an unbounded (for the most part) list of primes. 

Using <code>unfold f x | (a,y) <- f x = a ++ unfold f y</code>, <code>fix g = x where x = g x</code>, and <code>_Y g = g (_Y g)</code>  

(so that <code>unfold f = concat . Data.List.unfoldr (Just . f)</code> ).

<haskell>[n | n<-[2..], product [1..n-1] `rem` n == n-1] -- Wilson's theorem

nubBy (((>1).).gcd) [2..]

[n | n<-[2..], all ((> 0).rem n) [2..n-1]]
[n | n<-[2..], []==[i | i<-[2..n-1], rem n i==0]]
[n | n<-[2..], []==[j | i<-[2..n-1], j<-[i*i, i*i+i..n], j==n]]
[n | n<-[2..], all ((> 0).rem n) [2..floor.sqrt.fromIntegral$n]]
2 : [n | n<-[3,5..], all ((> 0).rem n) [3,5..floor.sqrt.fromIntegral$n]]

fix (\xs-> 2 : [n | n<-[3..], all ((> 0).rem n) $ takeWhile ((<= n).(^2)) xs])
2 : fix (\xs-> 3 : [n | n <- [5,7..], 
                            foldr (\p r-> p*p>n || (rem n p>0 && r)) True xs])
2:3 : fix (\xs-> 5 : [n | n <- scanl (+) 7 (cycle [4,2]), 
                            foldr (\p r-> p*p>n || (rem n p>0 && r)) True xs])
2:3:5 : fix (\xs-> 7 : [n | n <- scanl (+) 11 (cycle [2,4,2,4,6,2,6,4]), 
                            foldr (\p r-> p*p>n || (rem n p>0 && r)) True xs])

foldr (\x r-> x : filter ((> 0).(`rem`x)) r) [] [2..]    -- wrong, bottom-up
nub . map head . scanl (\xs x-> filter ((> 0).(`rem`x)) xs) [2..] $ [2..]
map head . iterate (\(x:xs)-> filter ((> 0).(`rem`x)) xs) $ [2..]
unfold (\(x:xs)-> ([x], filter ((> 0).(`rem`x)) xs)) [2..]  -- top-down

fix $ concatMap (fst.snd)
    . iterate (\(p:t,(h,xs)) -> (t,span (< head t^2) [y | y<-xs, rem y p>0]))
    . (, ([2,3],[4..]))                         --  (minus xs [p*p, p*p+p..])

[n | n<-[2..], not $ elem n [j*k | j<-[2..n-1], k<-[2..n-1]]]
[n | n<-[2..], not $ elem n [j*k | j<-[2..n`div`2], k<-[2..n`div`j]]]

zipWith (flip (!!)) [0..] . scanl1 minus              -- APL-style
                          . scanl1 (zipWith (+)) $ repeat [2..] 
tail . unfold (\(a:b:t) -> (:t) . (`minus` b) <$> span (< head b) a) 
     . scanl1 (zipWith (+) . tail) $ tails [1..]

2 : fix ((3:) . unfold (\(x,u:us)-> ([x|u==0], (x+2,us))) -- counting by 1s
              . (,) 5 . ([0,0] ++) 
              . foldi (\(x:xs) ys -> let n=div (head ys - x) 2 - 1 in
                            x:take n xs ++ zipWith max ys (drop n xs)) []
              . map (\p-> (p*p :) . tail . cycle $ 1 : replicate (p-1) 0) )

primesTo n = foldl (\a i-> minus a [i*i, i*i+2*i..n]) (2:[3,5..n]) 
                                       [3,5..(floor . sqrt . fromIntegral) n]
primesTo  = foldr (\n ps-> foldl (\a p-> minus a [p*p, p*p+p..n]) [2..n] ps)
                  [] . takeWhile (> 1) . iterate (floor . sqrt . fromIntegral)

2 : unfold (\(xs,p:ps)-> let (h,t)=span (< p*p) xs in 
                     (h, (filter ((> 0).(`rem`p)) t, ps))) ([3,5..],[3,5..])
2 : fix (\xs-> 3 : unfold (\(xs,p:ps)-> let (h,t)=span (< p*p) xs in
                     (h, ((`minus` [p*p, p*p+2*p..]) t, ps))) ([5,7..], xs))
2 : _Y (\ps-> concatMap snd $ iterate (\((fs:ft, x, p:t),_) ->       
                  ((ft,p*p+2,t), [x | x <- [x, x+2 .. p*p-2],
                   all ((/= 0).rem x) fs])) ((inits ps, 5, ps), [3]) ) 

2 : _Y (\ps-> concatMap snd $ iterate (\((fs:ft, x, p:t),_) ->       
                  ((ft,p*p+2,t),     minus [x, x+2 .. p*p-2] 
                   $ foldi union [] [[o, o+2*i .. p*p-2] | i <- fs, 
                   let o=x+mod(i-x)(2*i)])) ((inits ps, 5, ps), [3]) ) 

let { sieve (x:xs) = x : [n | n <- sieve xs, rem n x > 0] } -- wrong, bottom-up
           in sieve [2..] 
let { sieve (x:xs) = x : sieve [n | n <- xs, rem n x > 0] } -- top-down -
           in sieve [2..]                                     -- Turner's sieve

fix $ \xs-> let { sieve xs (p:ps) = let (h,t)=span (< p*p) xs in  -- postponed
                                   h ++ sieve (filter ((> 0).(`rem`p)) t) ps }
           in 2 : 3 : sieve [5,7..] (tail xs)
fix $ \xs-> let { sieve xs (p:ps) = let (h,t)=span (< p*p) xs in 
                                   h ++ sieve (t `minus` [p*p, p*p+2*p..]) ps }
           in 2 : 3 : sieve [5,7..] (tail xs)
                                           -- producing only unique composites:
fix $ \xs-> 2 : minus [3..] (foldr (\(p:ps) r-> fix $ ((p*p) :) .
                           merge r . map (p*) . merge ps) [] $ tails xs)
                                                    -- sieving by all odds:
2 : minus [3,5..] (foldi (\(x:xs)-> (x:).union xs) [] 
                              $ map (\x->[x*x, x*x+2*x..]) [3,5..])

unfold (\a@(1:p:_) -> ([p], a `minus` map (*p) a)) [1..]  -- Euler's sieve
unfold (\(p:xs) -> ([p], xs `minus` map (*p) (p:xs))) [2..]  

unfold (\(p:xs) -> ([p], minus xs [p, p+p..])) [2..]   -- the basic sieve
unfold (\xs@(p:_) -> (`minus` [p, p+p..]) <$> splitAt 1 xs) [2..] 
                                                       -- postponed
fix $ (2:) . unfold (\(p:ps,xs) -> ((ps ,) . (`minus` [p*p, p*p+p..]))  
                                   <$> span (< p*p) xs) . (, [3..]) 

                                    -- Richard Bird's combined sieve
fix $ (2:) . minus [3..] . foldr (\(x:xs)-> (x:) . union xs) [] 
                          . map (\p-> [p*p, p*p+p..])

2 : _Y ( (3:)                        -- unbounded tree-like folding
                . minus [5,7..]
                 . unionAll    -- ~= foldi (\(x:xs) ys-> x:union xs ys) []
                  . map (\p-> [p*p, p*p+2*p..]) )

[2,3,5,7] ++ _Y ( (11:)                                -- using a wheel;
                    . minus (scanl (+) 13 $ tail wh11) -- 1.4x faster than 
                     . unionAll                        --     the next one
                      . map (\p-> map (p*) . dropWhile (< p) $ 
                                    scanl (+) (p - rem (p-11) 210) wh11) )
[2,3,5,7] ++ _Y ( (11:)                   
                         . minus (scanl (+) 13 $ tail wh11)
                          . unionAll
                           . map (\(w,p)-> scanl (\c d-> c + p*d) (p*p) w)
                            . isectBy (compare . snd)
                                      (tails wh11 `zip` scanl (+) 11 wh11) )
_Y g = g (_Y g)
-- fix g = x where x = g x  -- from Data.Function
unfold f x | (a,y) <- f x = a ++ unfold f y
-- f <$> (a,b) = (a, f b)   -- from Data.Functor
wh11 = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:  
       4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wh11
</haskell>

<code>foldi</code> is an infinitely right-deepening tree folding function found [[Fold#Tree-like_folds|here]]. <code>minus</code> and <code>union</code> of course are on the main page [[Prime_numbers#Initial_definition|here]] (and <code>merge</code> a simpler, duplicates-ignoring variant of <code>union</code>), such that <code>xs `minus` ys == xs Data.List.\\ ys</code> for any finite increasing <code>xs</code> and <code>ys</code>.

The few last definitions use functions from the [http://hackage.haskell.org/package/data-ordlist <code>data-ordlist</code>] package, and the 2-3-5-7 wheel <code>wh11</code>; they might be leaking space unless <code>minus</code> is fused with its input (into [[Prime_numbers#Tree merging |<code>gaps</code>]]/[[Prime_numbers#Tree merging with Wheel|<code>gapsW</code>]] from the main page).

[[Category:Code]]