# [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)

Here is another prime sieve. It is about half the length of the fastest contributed code (ONeillPrimes) and nearly as fast until it blows up on garbage collection:
```
```
```% cat ONeillPrimes.hs | grep -v "^--" | wc
185    1105    6306
% cat BayerPrimes.hs  | grep -v "^--" | wc
85     566    2418```
```
[Integer] -O   1*10^6   2*10^6   3*10^6   4*10^6   5*10^6
---------------------------------------------------------
ONeillPrimes |  3.555 |  7.798 | 12.622 | 18.927 | 23.529
BayerPrimes  |  3.999 |  8.895 | 18.003 | 22.977 | 38.053
```
```
```
I wrote this as an experiment in coding data structures in Haskell's lazy evaluation model, rather than as explicit data. The majority of the work done by this code is done by "merge"; the multiples of each prime percolate up through a tournament consisting of a balanced tree of suspended "merge" function calls. In order to create an infinite lazy balanced tree of lists, the datatype
```
```
```data List a = A a (List a) | B [a]
```
```
```
is used as scaffolding. One thinks of the root of the infinite tree as starting at the leftmost child, and crawling up the left spine as necessary.
```
What I'm calling a "venturi"

```
```venturi :: Ord a => [[a]] -> [a]
```
```
```
merges an infinite list of infinite lists into one list, under the assumption that each list, and the heads of the lists, are in increasing order. This could be a generally useful function. If one can think of a better way to write "venturi", swapping in your code would in particular yield a faster prime sieve.
```
```
I found that a tertiary merge tree was faster than a binary merge tree, because this leads to fewer suspensions. One can speed this code up a bit by interleaving strict and lazy calls, but I prefer to leave the code short and readable.
```
It is bizarre that

```
```trim p = let f m x = mod x m /= 0 in filter (f p)
```
```
```
lurks in the prime sieve code, but it is only used with primes of size up to the square root of the largest output prime. I tried more thoughtful alternatives, and they all slowed down the sieve. Sometimes dumb is beautiful.
```
```
Thanks to apfelmus for various helpful remarks that lead me to think of this approach.
```
Here is the code:

```
```module BayerPrimes (venturi,primes) where

-- Code for venturi :: Ord a => [[a]] -> [a]

merge :: Ord a => [a] -> [a] -> [a] -> [a]
merge xs@(x:xt) ys@(y:yt) zs@(z:zt)
| x <= y = if x <= z
then x : (merge xt ys zs)
else z : (merge xs ys zt)
| otherwise = if y <= z
then y : (merge xs yt zs)
else z : (merge xs ys zt)
merge _ _ _ = undefined

data List a = A a (List a) | B [a]

mergeA :: Ord a => List a -> List a -> List a -> List a
mergeA (A x xt) ys zs = A x (mergeA xt ys zs)
mergeA (B xs)   ys zs = mergeB xs ys zs

mergeB :: Ord a => [a] -> List a -> List a -> List a
mergeB xs@(x:xt) ys@(A y yt) zs = case compare x y of
LT -> A x (mergeB xt ys zs)
EQ -> A x (mergeB xt yt zs)
GT -> A y (mergeB xs yt zs)
mergeB xs (B ys) zs = mergeC xs ys zs
mergeB _ _ _ = undefined

mergeC :: Ord a => [a] -> [a] -> List a -> List a
mergeC xs@(x:xt) ys@(y:yt) zs@(A z zt)
| x < y = if x < z
then A x (mergeC xt ys zs)
else A z (mergeC xs ys zt)
| otherwise = if y < z
then A y (mergeC xs yt zs)
else A z (mergeC xs ys zt)
mergeC xs ys (B zs) = B \$ merge xs ys zs
mergeC _ _ _ = undefined

root :: Ord a => List a -> [List a] -> [a]
root (A x xt) yss       = x : (root xt yss)
root (B xs) (ys:zs:yst) = root (mergeB xs ys zs) yst
root _ _ = undefined

wrap :: [a] -> List a
wrap [] = B []
wrap (x:xt) = A x \$ B xt

triple :: Ord a => [List a] -> [List a]
triple (x:y:z:xs) = mergeA x y z : (triple xs)
triple _ = undefined

group :: Ord a => [List a] -> [List a]
group (x:y:xt) = x : y : (group \$ triple xt)
group _ = undefined

venturi :: Ord a => [[a]] -> [a]
venturi (x:xt) = root (wrap x) \$ group \$ map wrap xt
venturi _ = undefined

-- Code for primes :: Integral a => [a]

diff  :: Ord a => [a] -> [a] -> [a]
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)
diff _ _ = undefined

trim :: Integral a => a -> [a] -> [a]
trim p = let f m x = mod x m /= 0 in filter (f p)

seed :: Integral a => [a]
seed = [2,3,5,7,11,13,17]

wheel :: Integral a => [a]
wheel = drop 1 [ m*j + k | j <- [0..], k <- ws ]
where m  = foldr1 (*) seed
ws = foldr trim [1..m] seed

multiples :: Integral a => [a] -> [[a]]
multiples ws = map fst \$ tail \$ iterate g ([], ws)
where g (_,ps@(p:pt)) = ([ m*p | m <- ps ], trim p pt)
g _ = undefined

primes :: Integral a => [a]
primes = seed ++ (diff wheel \$ venturi \$ multiples wheel)
```
```
_______________________________________________