# Prime numbers miscellaneous

### From HaskellWiki

(→One-liners: +1) |
(→One-liners: replace one) |
||

Line 140: | Line 140: | ||

Using | Using | ||

<haskell> | <haskell> | ||

− | (\\) = Data.List.\\ | + | (\\) = Data.List.\\ -- or Data.List.Ordered.minus, where appropriate |

_Y g = g (_Y g) | _Y g = g (_Y g) | ||

fix g = x where x = g x -- from Data.Function | fix g = x where x = g x -- from Data.Function | ||

Line 165: | Line 165: | ||

[n | let p n = []==[i | i<-[2..n-1], p i && rem n i==0], n<-[2..], p n] | [n | let p n = []==[i | i<-[2..n-1], p i && rem n i==0], n<-[2..], p n] | ||

-- (someone actually wrote this, though in C) | -- (someone actually wrote this, though in C) | ||

− | |||

− | |||

− | fix $ foldl (\\) [2..] . map (\p -> [2*p, 3*p..]) -- non-executable | + | fix $ foldl (\\) [2..] . map (\p -> [2*p, 3*p..]) -- non-executable specs |

+ | fix $ map head . scanl Data.List.Ordered.minus [2..] -- executable specs | ||

+ | . map (\p -> [p, p+p..]) | ||

[n | n<-[2..], all ((> 0).rem n) [2..floor.sqrt.fromIntegral$n]] | [n | n<-[2..], all ((> 0).rem n) [2..floor.sqrt.fromIntegral$n]] | ||

Line 300: | Line 300: | ||

<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> and <code>xs `union` ys == nub . sort $ xs ++ ys</code> for any finite, increasing lists <code>xs</code> and <code>ys</code>. | <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> and <code>xs `union` ys == nub . sort $ xs ++ ys</code> for any finite, increasing lists <code>xs</code> and <code>ys</code>. | ||

− | + | Some 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). | |

==A Tale of Sieves== | ==A Tale of Sieves== |

## Revision as of 18:06, 1 December 2017

For a context to this, please see Prime numbers.

## Contents |

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

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

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 [1]

nextSize (Wheel n rs) p = Wheel (p*n) [r2 | k <- [0..(p-1)], r <- rs, let r2 = n*k+r, r2 `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.

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 `(p-1)/p`

. See Euler's Sieve, or the functional pearl titled Lazy wheel sieves and spirals of primes for more.

## 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)*.

## 4 One-liners

Produce an unbounded (for the most part) list of primes.

Using

(\\) = Data.List.\\ -- or Data.List.Ordered.minus, where appropriate _Y g = g (_Y g) fix g = x where x = g x -- from Data.Function f <$> (a,b) = (a, f b) -- from Data.Functor unfold f a | (xs,b) <- f a = xs ++ unfold f b

(so that `unfold f = concat . Data.List.unfoldr (Just . f)`

):

[n | n<-[2..], product [1..n-1] `rem` n == n-1] -- Wilson's theorem Data.List.nubBy (((>1).).gcd) [2..] [n | n<-[2..], all ((> 0).rem n) [2..n-1]] -- trial division [n | n<-[2..], []==[i | i<-[2..n-1], rem n i==0]] -- implementing `all` [n | n<-[2..], []==[i | i<-[2..n-1], last [n, n-i..0]==0]] -- and `rem` -- is this still a trial division? [n | n<-[2..], []==[i | i<-[2..n-1], j<-[0, i..n], j==n]] -- or a "forgetful" ... sieve of Eratosthenes? [n | n<-[2..], []==[i | i<-[2..n-1], j<-[i*i, i*i+i..n], j==n]] [n | let p n = []==[i | i<-[2..n-1], rem n i==0], n<-[2..], p n] -- testing only by primes should be faster, right? [n | let p n = []==[i | i<-[2..n-1], p i && rem n i==0], n<-[2..], p n] -- (someone actually wrote this, though in C) fix $ foldl (\\) [2..] . map (\p -> [2*p, 3*p..]) -- non-executable specs fix $ map head . scanl Data.List.Ordered.minus [2..] -- executable specs . map (\p -> [p, p+p..]) [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]] -- optimal trial division fix (\xs-> 2 : [n | n<-[3..], all ((> 0).rem n) $ takeWhile ((<= n).(^2)) xs]) 2 : fix (\xs-> 3 : [n | n <- [5,7..], -- the wheel: 2, 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]), -- 3, 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]), -- 5 foldr (\p r-> p*p>n || (rem n p>0 && r)) True xs]) foldr (\x r-> x : filter ((> 0).(`rem`x)) r) [] [2..] -- 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) -- by spans between primes' squares . 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..]) -- segment-wise _Y (\ps-> 2 : [n | (r:q:_, px) <- (zip . tails . (2:) . map (^2) <*> inits) ps, n <- [r+1..q-1], all ((> 0).rem n) px]) _Y (\ps-> 2 : [n | (r:q:_, px) <- (zip . tails . (2:) . map (^2) <*> inits) ps, n <- [r+1..q-1] \\ [m | p <- px, m <- [0,p..q-1]]]) -- a sieve finds primes by eliminating multiples, isn't it? -- (found this one in the wild, too) [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..] fix (\sv -> dropWhile ((> 0).snd) >>> \((p,_):xs) -> -- counting by 1s p : sv (zipWith (fmap . max) (cycle $ [0 | _<-[2..p]] ++ [1]) xs)) . (`zip` repeat 0) $ [2..] 2 : fix ((3:) . unfold (\(x,u:us)-> ([x|u==0], (x+2,us))) -- counting by 2s . (,) 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 n = 2 : foldr (\r z-> if (head r^2) <= n then head r : z else r) [] (fix $ \rs-> [3,5..n] : [minus t [p*p, p*p+2*p..] | (p:t)<-rs]) primesTo = foldr (\n ps-> foldl (\a p-> -- with repeated square root minus a [p*p, p*p+p..n]) [2..n] ps) [] . takeWhile (> 1) . iterate (floor . sqrt . fromIntegral) concatMap (\(a:b:_)-> drop (length a) b) -- its dual, with repeated . tails . scanl (\ps n-> foldl (\a p-> -- squaring minus a [p*p, p*p+p..n]) [2..n] ps) [] $ iterate (^2) 2 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] } -- 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) 2 : minus [3,5..] (foldi (\(x:xs)-> (x:).union xs) [] -- sieving by all odds $ 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], xs `minus` [p, p+p..])) [2..] -- the basic sieve unfold (\xs@(p:_) -> (`minus` [p, p+p..]) <$> splitAt 1 xs) [2..] fix $ (2:) . unfold (\(p:ps,xs) -> (ps ,) . -- postponed sieve (`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-shaped 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:) -- with a wheel . minus (scanl (+) 13 $ tail wh11) . unionAll . map (\p-> map (p*) . dropWhile (< p) $ scanl (+) (p - rem (p-11) 210) wh11) ) [2,3,5,7] ++ _Y ( (11:) -- same, but 1.4x . minus (scanl (+) 13 $ tail wh11) -- slower than . unionAll -- the above . map (\(w,p)-> scanl (\c d-> c + p*d) (p*p) w) . isectBy (compare . snd) (tails wh11 `zip` scanl (+) 11 wh11) ) 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: -- 210-wheel, 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 -- from 11 -- cycle $ zipWith (-) =<< tail $ [i | i <- [11..221], gcd i 210 == 1]

`foldi`

is an infinitely right-deepening tree folding function found here. `minus`

and `union`

of course are on the main page here (and `merge`

a simpler, duplicates-ignoring variant of `union`

), such that `xs `minus` ys == xs Data.List.\\ ys`

and `xs `union` ys == nub . sort $ xs ++ ys`

for any finite, increasing lists `xs`

and `ys`

.

Some definitions use functions from the `data-ordlist`

package, and the 2-3-5-7 wheel `wh11`

; they might be leaking space unless `minus`

is fused with its input (into `gaps`

/`gapsW`

from the main page).

## 5 A Tale of Sieves

Old stuff, just put together in one place to better see the whole picture at once.

{-# LANGUAGE ViewPatterns #-} import Data.List ( (\\), tails, inits ) import Data.List.Ordered ( minus, unionAll ) import Data.Array.Unboxed primes = sieve [2..] where sieve (p:t) = [p] ++ sieve [n | n <- t, rem n p > 0] primes = 2 : sieve primes [3..] where sieve (p:ps) (span (< p*p) -> (h,t)) = h ++ sieve ps [n | n <- t, rem n p > 0] primes = 2 : [n | (r:q:_, px) <- (zip . tails . (2:) . map (^2) <*> inits) primes, n <- [r+1..q-1], all ((> 0).rem n) px] {- n <- [r+1..q-1] \\ concat [ [0,p..q] | p <- px]] n <- [r+1..q-1] `minus` unionAll [ dropWhile (<= r) [0,p..q-1] | p <- px]] (n,True) <- assocs ( accumArray (\_ _ -> False) True (r+1, q-1) [(m,()) | p <- px, let s = (r+p)`div`p*p, m <- [s,s+p..q-1]] :: UArray Int Bool )] -} primes = map head $ scanl minus [2..] [[p, p+p..] | p <- primes] -- = fix $ map head . scanl minus [2..] . map (\p -> [p, p+p..]) primes = 2 : 3 : [5,7..] `minus` unionAll [[p*p, p*p+2*p..] | p <- tail primes] --