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

== Implicit Heap ==

== Implicit Heap ==

The following is an original implicit heap implementation for the sieve of

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 <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.

+

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.

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.

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: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

+

</haskell>

+

+

Here, <hask>primes'</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) [r' | k <- [0..(p-1)], r <- rs,

+

let r' = n*k+r, r' `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 how about adapting the wheel size while generating prime numbers? 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.

<code>foldi</code> is an infinitely right-deepening tree folding function found [[Fold#Tree-like_folds|here]]. <code>minus</code> of course is on the main page [[Prime_numbers#Initial_definition|here]].

+

+

[[Category:Code]]

Revision as of 22:18, 28 March 2014

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.

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:

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.

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]

We can create a larger wheel by rolling a smaller wheel of circumference

3 Using IntSet for a traditional sieve

module Sieve whereimportqualified 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)