Why is my Clojure prime number lazy sequence so slow?

Question

I'm doing problem 7 of Project Euler (calculate the 10001st prime). I have coded a solution in the form of a lazy sequence, but it is super slow, whereas another solution I found on the web (link below) and which does essentially the same thing takes less than a second.

I'm new to clojure and lazy sequences, so my usage of take-while, lazy-cat rest or map may be the culprits. Could you PLEASE look at my code and tell me if you see anything?

The solution that runs under a second is here: https://zach.se/project-euler-solutions/7/

It doesn't use lazy sequences. I'd like to know why it's so fast while mine is so slow (the process they follow is similar).

My solution which is super slow:

(def primes 
  (letfn [(getnextprime [largestprimesofar]
    (let [primessofar (concat (take-while #(not= largestprimesofar %) primes) [largestprimesofar])]
      (loop [n (+ (last primessofar) 2)]
          (if
            (loop [primessofarnottriedyet (rest primessofar)]
              (if (= 0 (count primessofarnottriedyet))
                true
                (if (= 0 (rem n (first primessofarnottriedyet)))
                  false
                  (recur (rest primessofarnottriedyet)))))
            n
            (recur (+ n 2))))))]
    (lazy-cat '(2 3) (map getnextprime (rest primes)))))

To try it, just load it and run something like (take 10000 primes), but use Ctrl+C to kill the process, because it is too slow. However, if you try (take 100 primes), you should get an instant answer.

Answer 1

I will continue speeding it up, based on @manutter's solution.

(declare primes)

(defn is-relatively-prime? [n candidates]
  (if-let [current (first candidates)]
    (if (zero? (rem n current))
      false
      (recur n (rest candidates)))
    true))

(defn get-next-prime [largest-prime-so-far]
  (let [primes-so-far (concat (take-while #(not= largest-prime-so-far %) primes) [largest-prime-so-far])]
    (loop [n (+ (last primes-so-far) 2)]
      (if
        (is-relatively-prime? n (rest primes-so-far))
        n
        (recur (+ n 2))))))

(def primes
  (lazy-cat '(2 3) (map get-next-prime (rest primes))))

(time (first (drop 10000 primes)))

"Elapsed time: 14092.414513 msecs"

Ok. First of all let's add this current^2 > n optimization:

(defn get-next-prime [largest-prime-so-far]
  (let [primes-so-far (concat (take-while #(not= largest-prime-so-far %) primes) [largest-prime-so-far])]
    (loop [n (+ (last primes-so-far) 2)]
      (if
          (is-relatively-prime? n
                                (take-while #(<= (* % %) n)
                                            (rest primes-so-far)))
        n
        (recur (+ n 2))))))

user> (time (first (drop 10000 primes)))
"Elapsed time: 10564.470626 msecs"
104743

Nice. Now let's look closer at the get-next-prime:

if you check the algorithm carefully, you will notice that this: (concat (take-while #(not= largest-prime-so-far %) primes) [largest-prime-so-far]) really equals to all the primes we've found so far, and (last primes-so-far) is really the largest-prime-so-far. So let's rewrite it a little:

(defn get-next-prime [largest-prime-so-far]
  (loop [n (+ largest-prime-so-far 2)]
    (if (is-relatively-prime? n
          (take-while #(<= (* % %) n) (rest primes)))
      n
      (recur (+ n 2)))))

user> (time (first (drop 10000 primes)))
"Elapsed time: 142.676634 msecs"
104743

let's add one more order of magnitude:

user> (time (first (drop 100000 primes)))
"Elapsed time: 2615.910723 msecs"
1299721

Wow! it's just mind blowing!

but that's not all. let's take a look at is-relatively-prime function: it just checks that none of the candidates evenly divides the number. So it is really what not-any? library function does. So let's just replace it in get-next-prime.

(declare primes)

(defn get-next-prime [largest-prime-so-far]
  (loop [n (+ largest-prime-so-far 2)]
    (if (not-any? #(zero? (rem n %))
                  (take-while #(<= (* % %) n)
                              (rest primes)))
      n
      (recur (+ n 2)))))

(def primes
  (lazy-cat '(2 3) (map get-next-prime (rest primes))))

it is a bit more productive

user> (time (first (drop 100000 primes)))
"Elapsed time: 2493.291323 msecs"
1299721

and obviously much cleaner and shorter.

Answer 2

Let me re-write your code just a bit to break it down into pieces that will be easier to discuss. I'm using your same algorithm, I'm just splitting out some of the inner forms into separate functions.

(declare primes)   ;; declare this up front so we can refer to it below

(defn is-relatively-prime? [n candidates]
  (if (= 0 (count candidates))
    true
    (if (zero? (rem n (first candidates)))
      false
      (is-relatively-prime? n (rest candidates)))))

(defn get-next-prime [largest-prime-so-far]
  (let [primes-so-far (concat (take-while #(not= largest-prime-so-far %) primes) [largest-prime-so-far])]
    (loop [n (+ (last primes-so-far) 2)]
      (if
        (is-relatively-prime? n (rest primes-so-far))
        n
        (recur (+ n 2))))))

(def primes
  (lazy-cat '(2 3) (map get-next-prime (rest primes))))

(time (let [p (doall (take 200 primes))]))

That last line is just to make it easier to get some really rough benchmarks in the REPL. By making the timing statement part of the source file, I can keep re-loading the source, and get a fresh benchmark each time. If I just load the file once, and keep trying to do (take 500 primes) the benchmark will be skewed because primes will hold on to the primes it has already calculated. I also need the doall because I'm pulling my prime numbers inside a let statement, and if I don't use doall, it will just store the lazy sequence in p, instead of actually calculating the primes.

Now, let's get some base values. On my PC, I get this:

Loading src/scratch_clojure/core.clj... done
"Elapsed time: 274.492597 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 293.673962 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 322.035034 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 285.29596 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 224.311828 msecs"

So about 275 milliseconds, give or take 50. My first suspicion is how we're getting primes-so-far in the let statement inside get-next-prime. We're walking through the complete list of primes (as far as we have it) until we get to one that's equal to the largest prime so far. The way we've structured our code, however, all the primes are already in order, so we're effectively walking thru all the primes except the last, and then concatenating the last value. We end up with exactly the same values as have been realized so far in the primes sequence, so we can skip that whole step and just use primes. That should save us something.

My next suspicion is the call to (last primes-so-far) in the loop. When we use the last function on a sequence, it also walks the list from the head down to the tail (or at least, that's my understanding -- I wouldn't put it past the Clojure compiler writers to have snuck in some special-case code to speed things up.) But again, we don't need it. We're calling get-next-prime with largest-prime-so-far, and since our primes are in order, that's already the last of the primes as far as we've realized them, so we can just use largest-prime-so-far instead of (last primes). That will give us this:

(defn get-next-prime [largest-prime-so-far]
  ; deleted the let statement since we don't need it
  (loop [n (+ largest-prime-so-far 2)]
    (if
      (is-relatively-prime? n (rest primes))
      n
      (recur (+ n 2)))))

That seems like it should speed things up, since we've eliminated two complete walks through the primes sequence. Let's try it.

Loading src/scratch_clojure/core.clj... done
"Elapsed time: 242.130691 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 223.200787 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 287.63579 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 244.927825 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 274.146199 msecs"

Hmm, maybe slightly better (?), but not nearly the improvement I expected. Let's look at the code for is-relatively-prime? (as I've re-written it). And the first thing that jumps out at me is the count function. The primes sequence is a sequence, not a vector, which means the count function also has to walk the complete list to add up how many elements are in it. What's worse, if we start with a list of, say, 10 candidates, it walks all ten the first time through the loop, then walks the nine remaining candidates on the next loop, then the 8 remaining, and so on. As the number of primes gets larger, we're going to spend more and more time in the count function, so maybe that's our bottleneck.

We want to get rid of that count, and that suggests a more idiomatic way we could do the loop, using if-let. Like this:

(defn is-relatively-prime? [n candidates]
  (if-let [current (first candidates)]
    (if (zero? (rem n current))
      false
      (recur n (rest candidates)))
    true))

The (first candidates) function will return nil if the candidates list is empty, and if that happens, the if-let function will notice, and automatically jump to the else clause, which in this case is our return result of "true." Otherwise, we'll execute the "then" clause, and can test for whether n is evenly divisible by the current candidate. If it is, we return false, otherwise we recur back with the rest of the candidates. I also took advantage of the zero? function just because I could. Let's see what this gets us.

Loading src/scratch_clojure/core.clj... done
"Elapsed time: 9.981985 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 8.011646 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 8.154197 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 9.905292 msecs"
Loading src/scratch_clojure/core.clj... done
"Elapsed time: 8.215208 msecs"

Pretty dramatic, eh? I'm an intermediate-level Clojure coder with a pretty sketchy understanding of the internals, so take my analysis with a grain of salt, but based on those numbers, I'd guess you were getting bitten by the count.

There's one other optimization the "fast" code is using that yours isn't, and that's bailing out on the is-relatively-prime? test whenever current squared is greater than n---you might speed up your code some more if you can throw that in. But I think count is the main thing you're looking for.