Lazy Sieve Eratosthenes in Scala and Clojure

Sieve Eratosthenes is a popular prime number generation algorithm.  Start with 2 and cross out all even numbers until a limit number. Then, doing the same for 3 to cross out all factors of 3, like 3, 6, 9.  The goal of the blog is to show some lazy sieve algorithms using clojure and scala and show the benchmark result of the algorithms. 
The complete code is available at https://github.com/anrizal06/rizier .


It is easy to implement eager Sieve algorithms, using scala, this can be achieved as follow:
1 def sieve(limit: Int) = {

   2    
   3    val composites = new Array[Boolean](limit)
   4    
   5    (2 to math.sqrt(limit).toInt).foreach(i => {
   6        if (!composites(i)) 
   7            (i * i  until limit by i).foreach(composites(_) = true)
   8    })
   9    
  10    for { i <- 2 until limit
  11          if (!composites(i))
  12    } yield i
  13 }

On of the problem with the solution above is that the limit parameter that must be used. The parameter is OK if we want to get prime numbers below a certain limit, but it is not OK when we want, for example to get the 1000th prime.

The solution for the problem above presented here is by using lazy evaluation provided in both scala and clojure. Using lazy algorithms, we generate a prime number, only when necessary. For example, if we want to get the 26th prime number, we only generate 26 prime number, not more.

The lazy algorithms for sieve are presented in the work of Melissa ONeill (2009) and the idea are adopted by for example clojure contrib API. The paper also showed that the naive functional solution (not exposed here) is very limited. Indeed, I made some tests with the solution and reached easily StackOverflow error for relatively small number.

Not So Naive Solution

Let's start first with the still naive, but optimized somehow, lazy sieve eratosthenes. First in clojure:

   1 (defn naive-primes
   2    "Implementation of lazy but naive sieve algorithms"
   3    []
   4   (letfn [
   5           (next-prime [x xs] 
   6              (if (some #(zero? (rem x %))
   7                     (take-while #(<= (* % %) x) xs))
   8              (recur (+ x 2) xs)
   9              (cons x (lazy-seq (next-prime (+ x 2) (conj xs x))))))]
  10    (cons 2 (lazy-seq (next-prime 3 [])))))

Line 10 starts the implementation by consing 2 to the lazy sequence of prime number starting from 3. Then, starting from 3, all odd numbers are examined if they are prime or not. This is done at line 6-7 of the code. If it is a composite number, then it recursively calls next prime to examine the next odd number (line 8) . Otherwise, if it is a prime, then the prime number is returned, and the list of prime numbers discovered so far is updated.

To get the 100th prime, one can use (nth (naive-primes) 100)

The equivalent Scala program is the following:

 1   def primes(): Stream[Int] = {
 2      lazy val ps = 2 #:: sieve(3)
 3      def sieve(p: Int): Stream[Int] = {
 4          p #:: sieve(
 5             Stream.from(p + 2, 2).
 6                 find(i=> ps.takeWhile(j => j * j <= i).
 7                 forall(i % _ > 0)).get)
 8      }
 9      ps
10   } 

Stream is the Scala equivalent of Clojure lazy-seq. Note the interesting use of lazy val. To get the 100th prime number, one can use primes().take(100).last

At stack overflow , I exposed another solution using mutable list as the alternative solution to the lazy stream above. It turned out that the mutable list solution was faster than the solution above.

With Wheel

The Clojure Contrib library implements Sieve algorithm using wheel. The idea of the algorithm is to discard as many as possible the composite number that is factor of first prime numbers (3, 5, 7 in addition to 2). Indeed, the first solution already discards 2. To discard the numbers, a circular list of step -- or wheel -- is usd. For (2, 3, 5, 7), the wheel is 

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

To make the list circular, clojure has a very convenient method call cycle that can then be used to create the circular list:

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


Scala does not have built-in cycle, but it is quite straight forward to implement cycle in scala:

 def cycle[A](xs: List[A]): Stream[A] = {
     lazy val result: Stream[A] = xs.toStream.append(result)
     result
   }

The complete implementation of the algorithm in Clojure is:

   1 (defn wheel-primes 
   2    "The optimization of lazy naive sieve algorithm 
   3    by discarding the number divisible by 2, 3, 5, and 7.
   4    This is done by adding the next number to be examined
   5    by the number given by the wheel."
   6    []
   7    (let [next-prime
   8          (fn next-prime [x xs [f & r]] 
   9              (if (some #(zero? (rem x %))
  10                     (take-while #(<= (* % %) x) xs))
  11              (recur (+ x f) xs r)
  12              (cons x (lazy-seq (next-prime (+ x f) (conj xs x) r)))))
  13          wheel (cycle [2 4 2 4 6 2 6 4 2 4 6 6 2 6  4  2
  14                         6 4 6 8 4 2 4 2 4 8 6 4 6 2  4  6
  15                         2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10])]
  16          (concat 
  17              [2 3 5 7]
  18              (lazy-seq (next-prime 11 [] wheel)))))

The second solution is then close to the first one, except that instead of incrementing with 2, the next number to be examined is obtained by incrementing with the number in the wheel.
The Scala version of the solution is: 

   1    def cycle[A](xs: List[A]): Stream[A] = {
   2      lazy val result: Stream[A] = xs.toStream.append(result)
   3      result
   4    }
   5    def wheel(): Stream[Int] = {
   6      cycle(List(
   7        2 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 
   8        4 , 6 , 6 , 2 , 6 , 4 , 2, 6 , 4 , 6 , 
   9        8 , 4 , 2 , 4 , 2 , 4 , 8 , 6 , 4 , 6 , 
  10        2 , 4 , 6, 2 , 6 , 6 , 4 , 2 , 4 , 6 , 
  11        2 , 6 , 4 , 2 , 4 , 2 , 10 , 2 , 10))
  12    }
  13  
  14    def from(x: Int, stepIt: Iterator[Int]): Stream[Int] = {
  15        x #:: from(x + stepIt.next, stepIt)
  16    }
  17    
  18    def primes(): Stream[Int] = {
  19        lazy val ps = 2 #:: 3 #:: 5 #:: 7 #:: sieve(11, wheel.iterator)
  20        def sieve(p: Int, it: Iterator[Int]): Stream[Int] = {
  21            p #:: sieve( from(p + it.next, it).
  22                            find(i => ps.takeWhile(j => j * j <= i).
  23                             forall(i % _ > 0)).get,
  24                          it)
  25        }
  26        ps
  27    }






Hash Map Based Solution
The problem with the first two solutions is that both solutions require the check whether a number is factor of the prime number discovered so far.  For example, when examining 121, the first two algorithms verify if 121 are factor 3, 5, and 7.  It turns out, that the check is not necessary, and this leads to the third solution of lazy sieve. 

Like the two solutions above, the algorithm receives a stream of odd number and returns a stream of prime number.  The algorithm uses a hash map that stores the last discovered composite number. Each composite number that is stored in the hash map is a factor of prime numbers discovered so far. The key of the hash map is the composite number and the value is the double of the corresponding prime number. 

When a number is under examination, it is checked against the hash map. 
If the number is in the hash map, obviously, it is a composite number. In this case, the number is removed from the hash map and replaced by new composite number that is factor of the corresponding prime. The new composite number must not be in the hash map, therefore it continues to increment the number until discovering a brand new composite number.

On the other hand, when the number is not in the hash map since the beginning, it is a prime. A new composite corresponding to the prime is then inserted to the hashmap; the number is the square of the prime. 


  1 (defn hashtable-primes 
  2    "An implementation of lazy sieve algorithm by storing in a hastable
  3    the next composite number corresponding to a prime number"
  4    []
  5   (letfn [
  6     (next-composite [x step sieve]
  7        (if (sieve x)
  8           ; the composite number has already been found previously
  9           (recur (+ x step) step sieve)
 10           ; brand new composite number
 11           (assoc sieve x step)))
 12     (next-prime [x sieve] 
 13       (let [step (sieve x)]         
 14        (if (sieve x)
 15           ; found a composite
 16           (recur (+ x 2) (next-composite (+ x step) step (dissoc sieve x)))
 17           ; found a prime
 18           (cons x 
 19             (lazy-seq 
 20                (next-prime (+ x 2) (assoc sieve (* x x) (* x 2))))))))]
 21   (cons 2 (lazy-seq (next-prime 3 {})))))


The Scala version of the solution is the following. I used mutable map to gain better performance without really sacrifying the functional nature of the code:






 1     import scala.collection.mutable.Map
 2     
 3     def primes() : Stream[Int] = {
 4         2 #:: sieve(3, Map{ 9 -> 6 })
 5     }
 6     private def sieve(p: Int, pQ: Map[Int, Int]): Stream[Int] = {
 7       p #:: sieve(nextPrime(p + 2, pQ), pQ )
 8     }
 9     private def nextCompositeNumber(
10           x:Int, step: Int, pQ: Map[Int, Int]): Unit = {
11       pQ.get(x) match {
12         case Some(_) => nextCompositeNumber(x + step, step, pQ)
13         case None => pQ(x) = step
14       } 
15     }
16     private def nextPrime(candidate: Int, pQ: Map[Int, Int]): Int = {
17       pQ.get(candidate) match {
18           case Some(step) =>
19               pQ -= candidate
20               nextCompositeNumber(candidate + step, step, pQ) 
21               nextPrime(candidate + 2, pQ)
22           case None =>
23               pQ(candidate * candidate) = candidate * 2
24               candidate 
25       }
26     }





Benchmark

To complete the implementation, I did a simple benchmark, mainly with the Scala versions, but Clojure version showed also the similar results. For the implementation, I used Scala 2.8.1 and Clojure 1.2 with Oracle JVM 1.6.0_14 version. The benchmark was done in my tiny laptop with 3GB of memory, Windows, and Intel Pentium Dual Core. 

The result is shown in the following graph: 






Where x-axis shows the limit in million and the y axis is the execution time in ms. 
The blue line is for the naive solution, and magenta one for wheel solution, and the yellow one for the hash map one. 

The same curve, but not the same execution time, is also observed in the implementation using Clojure. However, The Clojure implementation execution time is longer than Scala.

The wheel solution improves the performance of the naive one, but the implementation that uses hashmap performs much better. 

The maximum limit I could manage to have in my machine using the Scala implementation of the hash map algorithm was 150 000 000 that represented more than 8 million primes . The time of execution was around 270 s.






Reference
ONeill M (2009) The Genuine of Sieve Erathosthenes. Available at: 
http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf