A couple weeks ago I published a blog post on computing prime numbers with Clojure. I got some good feedback on it, including pointers to other implementation techniques I'd overlooked.

In this post, I want to explore some of those, notably two based on the Java standard library.

Recap

First off, let's recap where we left off in the first post: we had two approaches, one based on checking each integer for primality, and the other, based on the Sieve of Eratosthenes, based on eliminating non-primes so that we're left with only prime numbers (without having to check them).

In terms of implementation, here's what they might look like. First, the primality check on each number:

(defn check-primality
  []
  (let [p (promise)]
    (deliver p (cons 2
                     (->> (iterate #(+ 2 %) 3)
                          (remove (fn [n]
                                    (->> @p
                                         (take-while #(<= (* % %) n))
                                         (some #(zero? (rem n %)))))))))
    @p))

We want to produce a list of all prime numbers, so it has to be lazy. In this case, both iterate and remove are lazy. We know from the commutativity of multiplication that we don't need to check for divisors larger than the square root of the number we're testing if there were no divisor smaller than it (take-while); we also know that we don't need to check for non-prime divisors, which is why we need to introduce a promise there to allow for the remove operation to refer to the list of primes we are in the process of constructing. This does not result in an infinite loop because, at run time, the remove operation will never consume more primes than we've already produced (thanks to the take-while condition). Finally, we use some to stop checking for divisors as soon as we find one.

Second, we had an approach based on the Sieve of Eratosthenes:

(defn sieve-based
  [sieve]
  (letfn [(helper [known-primes already-produced sieve-bound]
            (concat (drop already-produced known-primes)
                    (lazy-seq
                      (loop [new-bound (* 10 sieve-bound)]
                        (let [new-primes (sieve new-bound)]
                          (if (> (count new-primes) (count known-primes))
                            (helper new-primes
                                    (count known-primes)
                                    new-bound)
                            (recur (* 10 new-bound))))))))]
    (helper (sieve 10) 0 10)))

(defn array-sieve
  [bound]
  (let [candidates (long-array bound)]
    ;; initialize array to contain [0 0 2 3 4 5 ...]
    (loop [idx 2]
      (when (< idx bound)
        (aset-long candidates idx idx)
        (recur (inc idx))))
    ;; eliminate known non-primes
    (loop [idx 2]
      (when (< idx bound)
        (let [next-candidate (aget candidates idx)]
          (when-not (zero? next-candidate)
            (loop [update-idx (* idx idx)]
              (when (< update-idx bound)
                (aset-long candidates update-idx 0)
                (recur (+ idx update-idx)))))
          (recur (inc idx)))))
    ;; return a Clojure vector with only primes in it
    (vec (remove zero? candidates))))

This is a bit more code, but the principle is fairly simple: we start with a chalkboard where we've written down all of the naturals starting at 2 and up to some limit (because chalkboards in Ancient Greece were not infinite), then declare 2 a prime, and remove all the multiples of 2 from the board. We don't need any test for that, we just increment our index by 2 each time. Then we take the next number still on the board (in our implementation, the next non-zero number), and keep going like that. We start the "erasure" process at (* n n) instead of (* 2 n) because, by construction, we've already eliminated all (* m n) numbers for any (< m n). It's very annoying that we have to pick a bound for the array-based sieve method to work, and there's going to be some repeated work as we increase the bound, but it's still very fast despite that.

Note: In my previous post, I screwed up the implementation of generate-sieve by hardcoding sieve-upto (instead of using the argument sieve-fn). That makes the switch from iterate to an explicit lazy seq arguably unnecessary, but I still like this form better.

Benchmarking

Just playing with various ways to define the prime number computation is fun, but I find it a lot more interesting to compare them for speed. As a first approximation, we could just time them, as I did in my previous post, but this time I'd like to go for something a bit more principled, and will use the criterium library instead, running in a standalone application (rather than directly at the REPL, which disables some JVM optimizations).

(defmacro bench
  [primes-fn]
  `(doseq [n# [100000 300000 1000000]]
     (println (format "%-28s %7d: %5.2f"
                      (pr-str '~primes-fn)
                      n#
                      (->> (crit/benchmark (nth ~primes-fn n#) {})
                           :mean first)))))

BigInteger

One possibly useful method from the Java standard library is BigInteger's isProbablePrime method. It's documented as:

Returns true if this BigInteger is probably prime, false if it's definitely composite. If certainty is ≤ 0, true is returned.

Given that this is a probabilistic approach, with a given certainty, I would expect this to be faster, for large non-primes, than actually checking for all divisors. And because it is guaranteed to return true for primes, we can use it as a filter with no semantic change to our existing check-based implementation:

(defn check-based-with-bigint
  [certainty]
  (let [p (promise)]
    (deliver p (cons 2
                     (->> (iterate #(+ 2 %) 3)
                          (filter (fn [n]
                                    (.isProbablePrime
                                      (java.math.BigInteger. n)
                                      certainty)))
                          (remove (fn [n]
                                    (->> @p
                                         (take-while #(<= (* % %) n))
                                         (some #(zero? (rem n %)))))))))
    @p))

Does it make it any faster?

t.core=> (time (nth (check-based) 1000000))
"Elapsed time: 107905.383266 msecs"
15485867
t.core=> (time (nth (check-based-with-bigint 2) 1000000))
"Elapsed time: 110825.211513 msecs"
15485867
t.core=> (time (nth (check-based-with-bigint 5) 1000000))
"Elapsed time: 119101.514406 msecs"
15485867
t.core=>

Oh well. They can't all be winners.

BitSet

Another possibly-useful class in the Java standard library is BitSet: in our Sieve implementation, we only really need to keep track of which elements are crossed; the actual value of the element is not useful because, if it's not zero, it's always equal to its index. So why bother storing an integer at all? We could just store a bit for that. That's exactly what BitSet lets us do.

Here is the sieve implemented on top of BitSet instead of an array:

(defn bitset-sieve
  [^long bound]
  (let [candidates (java.util.BitSet. bound)]
    ;; initialize array to contain [0 0 2 3 4 5 ...]
    (.set candidates 2 bound)
    ;; eliminate known non-primes
    (loop [idx 2]
      (when (< idx bound)
        (when-let [next-candidate (.get candidates idx)]
          (loop [update-idx (* idx idx)]
            (when (< update-idx bound)
              (.clear candidates update-idx)
              (recur (+ idx update-idx)))))
        (recur (inc idx))))
    ;; return a Clojure seq with only primes in it
    (take-while pos? (iterate #(.nextSetBit candidates (inc %)) 2))))

If we compare this to the array-based implementation, it's really not more complicated, just using a slightly less common API. How does it do?

t.core=> (time (last (array-sieve 10000000)))
"Elapsed time: 4214.152658 msecs"
9999991
t.core=> (time (last (bitset-sieve 10000000)))
"Elapsed time: 391.648236 msecs"
9999991
t.core=>

So this one does seem promising. Let's now run a real benchmark with all that.

Benchmark results

Here's the benchmark we're running, based on the definitions above:

(defn -main
  [& args]
  (bench (check-based))
  (bench (check-based-with-bigint 2))
  (bench (check-based-with-bigint 5))
  (bench (sieve-based array-sieve))
  (bench (sieve-based bitset-sieve)))

and here's the result on my machine:

$ java -server -jar target/uberjar/t-app-standalone.jar
(check-based)                 100000:  1.63
(check-based)                 300000:  7.73
(check-based)                1000000: 46.26
(check-based-with-bigint 2)   100000:  1.99
(check-based-with-bigint 2)   300000:  9.20
(check-based-with-bigint 2)  1000000: 47.43
(check-based-with-bigint 5)   100000:  2.32
(check-based-with-bigint 5)   300000:  9.77
(check-based-with-bigint 5)  1000000: 49.19
(sieve-based array-sieve)     100000:  2.87
(sieve-based array-sieve)     300000:  2.80
(sieve-based array-sieve)    1000000: 29.51
(sieve-based bitset-sieve)    100000:  0.18
(sieve-based bitset-sieve)    300000:  0.20
(sieve-based bitset-sieve)   1000000:  5.18
$

I must say I did not expect these numbers to be quite so much lower than the REPL ones, but the overall rankings seem to hold.

Conclusion

While the BigInteger route did not yield any fruit, the BitSet class gave us a very significant speed boost, at the cost of a very minor, local change to the code of array-sieve. While this post deals with an eminently useless problem, there's a very practical, broader lesson here. Clojure being built on the JVM was a very deliberate choice, and it was not just done to avoid having to reimplement a garbage collector. Integration with the JVM and its ecosystem is a cornerstone of Clojure's design, and it's often worth looking around in existing Java libraries for solutions to problems encountered in Clojure.