4

I'm learning clojure by going through project euler and am working on problem number 10 (find the sum of all the prime number below two million. I implemented a pretty literal algorithm for the sieve of eratosthenes but it works far too slowly to be useful for up to two million. I tried implementing it with loop-recur to not create any new frames but that didn't have a big impact on performance.

(defn find-primes-sum [last-prime nums]
    (loop [p last-prime n nums sum 0]
        (println p)
        (if (empty? n)
        sum
            (recur (first n) (doall (remove #(zero? (mod % (first n))) n)) (+ sum (first n))))))

(defn sieve-primes-until [limit]
    (find-primes-sum 2 (filter odd? (range 2 (inc limit)))))

(println (sieve-primes-until 2000000))
Scott
  • 41
  • 1

2 Answers2

2
(set! *unchecked-math* true)

(defmacro iloop [[b t n] & body]
  `(loop [~@b]
     (when ~t
       ~@body
       (recur ~n))))

(defn count-primes [^long n]
  (let [c (inc n)
        ^booleans prime? (make-array Boolean/TYPE c)]
    (iloop [(i 2) (<= i n) (inc i)]
      (aset prime? i true))
    (iloop [(i 2) (<= (* i i) n) (inc i)]
      (if (aget prime? i)
        (iloop [(j i) (<= (* i j) n) (inc j)]
          (aset prime? (* i j) false))))
    (areduce prime? i r 0
      (if (aget prime? i)
        (inc r)
        r))))

This version targets Clojure 1.3.0 alpha. It counts primes up to 1e8 in 2 seconds on my machine. It can be easily altered to collect them. It was originally written to show that you can implement the sieve so that it runs as fast as the comparable Java.

http://dosync.posterous.com/lispers-know-the-value-of-everything-and-the

dnolen
  • 18,496
  • 4
  • 62
  • 71
1

Personally, I'd structure the code differently. I've got an implementation of this problem that generates a lazy seq of all primes first, then sums the first 2,000,000 items. Takes 16 seconds on clojure 1.2, but it's hardly optimized except for using recur.

You're also doing a lot of unnecessary comparing with your input range; (remove #(zero? (mod % (first n))) n) tests all candidates against all odd numbers, which is nonsense**, especially when you force it with doall. You actually only have to test candidate prime x against all known primes <= sqrt(x) and discard the candidate when you find your first match.

** I just noticed your algorithm is not that stupid, but I'd still recommend you rewrite your algorithm as a lazy seq finds "the next prime" given a seq of previous primes and a candidate number.

Joost Diepenmaat
  • 17,633
  • 3
  • 44
  • 53
  • Abstraction is good on the eye and good for optimization. It's more intuitive to write smaller and better functions than optimize a complex one. So as @Joost-Diepenmaat said, at least separate the prime generation and the summing. Don't dwell on this particular problem too much though. Some people write sieve algo for their PhD! – Paul Lam Jun 17 '11 at 13:53
  • I'm not sure if you understand what the code is trying to do (should've commented!). When it recurs, the range of numbers is replaced with one that has had all of the multiples of the last prime removed because those are obviously cannot be prime numbers or are useful as a factor. The `doall` is there to force evaluation of the `remove #(zero? (mod % (first n))) n)` to prevent a stack overflow on laziness. – Scott Jun 17 '11 at 14:56
  • Another note, if I re-wrote it to find the next prime given all previous numbers I would I need to test every number in the range or is there a way to restrict what numbers I need to test (besides obviously using a range with only even numbers) – Scott Jun 17 '11 at 15:00
  • I do understand it, I think, but as far as I can see, doing the algorithm this way around implies doing a lot more divisibility tests because you're starting out with a huge list of candidates that you slowly weed down. If you define a prime number X as one that's indivisible by all prime numbers <= sqrt X, you should end up with a faster algorithm. At least, my implementation is a lot quicker than yours and it doesn't do anything else that's particularly clever. – Joost Diepenmaat Jun 17 '11 at 15:09
  • Re: your second question: you only need to test the previous *primes*. – Joost Diepenmaat Jun 17 '11 at 15:09
  • Whoops, I meant given all previous _prime_ numbers. Do I have to test every number incrementally from the last prime I've found (like having previous primes {2, 3, 5, 7} I need to test 8, 9, 10, and 11 until I get the next prime) or is there a way to restrict what numbers I need to test for primeness? – Scott Jun 17 '11 at 15:18
  • Well, you can always go from the last prime + 2, since even numbers aside from 2 aren't prime. There's a few optimizations you can do there, to exclude a few more obvious non-primes, but the main advantage is that to test if that next number is prime, you only need to test if it's divisible by any of the previous primes (and you can stop testing if it is, so work from 2 upwards to sqrt X - the goal is always to quickly discard non-primes). – Joost Diepenmaat Jun 17 '11 at 16:24