Project Euler Answers

To start off, and to make our code more clear, we define a small predicate called divisible? which checks if a number is divisible by some other number. While it doesn't really shorten the code all that much, it adds clarity to the code

(defn divisible? [num div] (zero? (mod num div)))

The most elegant solution would be the sieve of Eratosthenes. The logic behind it is fairly simple:

1. all natural numbers are either primes or a composite of a number of primes
2. a number is prime if no prime divisor less than itself exists
3. therefore, if one creates a list from \(2...n\) and takes the head of the list and filter out all divisors, add the head to a second list, and continue until the list is emptied the second list should contain the primes from \(2...n\)

As you can see, it follows that the list generated will be entirely comprised of prime numbers.

To solve this we can implement a lazy sieve of Eratosthenes. Lazy in that the resulting sequence will be infinite.

(defn seive
  ([] (seive (drop 2 (range))))
  ([l]
   (lazy-seq
    (let [[x & xs] l]
      (if (empty? xs)
        [x]
        (cons x (seive (filter #(not (divisible? % x)) xs))))))))

(def primes (seive))

Unfortunately this does not work for what we want. It is very beautiful, and I wish it did work, but because of the specific features of TCO in Clojure ² 2 Even with the recur special form. it doesn't actually produce TCO code, and, because it is recursive, will end up causing a stack overflow.

However, because of the lazy list semantics you can take a nearly infinite series of primes from it, assuming you only increase your requested list size by one each time you take from it. This is interesting, though unfortunately not something I can really use.

So in order to solve this we can rewrite the is-prime function to, rather than generating a list from 0 to \(n\) and then filtering it, simply loop until a value at the range of the prime number is found. This ensures that the space complexity of it is \(O(\log(n))\) rather than \(O(n)\), though does not change the time complexity.

First we define an atom to memoize our primes called found-primes, and initialize it with the first prime number, 2.

(def ^:private found-primes
  "A map of all primes that have been found thus far by the system"
  (atom #{2}))
(def ^:private max-prime (atom 2))

Next we define an is-prime? predicate which will test if a value is prime by iterating over the list of primes that are less than the square root of the new value. This is basically the same as the sieve of Eratosthenes applied to a single element. Below is the function itself, don't worry, we will break it down.

(defn is-prime?
  ([n] (is-prime? n false))
  ([n contiguous?]
   (let [max-div (int (inc (Math/sqrt n)))]
     (letfn [(div? [n]
               (->> @found-primes
                    (take-while #(> max-div %))
                    (some #(divisible? n %))
                    not))]
       (cond
         (contains? @found-primes n)
         true

         (> @max-prime n)
         false

         (< @max-prime max-div)
         (do (doall (map #(is-prime? % true) (range @max-prime max-div)))
             (if (div? n)
               (do (swap! found-primes conj n)
                   true)
               false))

         (div? n)
         (do (if contiguous?
               (do (swap! found-primes conj n)
                   (reset! max-prime n))
               (swap! found-primes conj n))
             true)

         :else false)))))

There are, however, a number of fun tricks being used to speed up the computation. First we generate the maximum divisor using this little bit of code. This isn't really the maximum size of the prime divisor, but rather the maximum size of the minimum element of the pair of divisors produced.

(int (inc (Math/sqrt n)))

Next we have the conditional. First it checks if the number is in the list of found primes, and if it is then it simply returns true. This means that it has \(O(1)\) runtime should the prime be found already. ³ 3 Though in Clojure, since sets are trees, it has more \(O(log(n))\) runtime.

(contains? @found-primes n)
      true

But, if it is both not in the list of found primes and is less than the largest prime, it must be the case that it is not a prime, and therefore it returns false.

(> @max-prime n)
false

Following this we check the case that the max-prime is smaller than max-div. In this case we check the primality of every prime from the last largest prime to the max-div. This places each prime number in the set of found primes and ensures that it can run and check if the number is prime. Following that we check that n is prime.

We do add it to the list of primes, but we don't set the max-prime, even though it far exceeds it. This is because the value max-prime represents the maximum prime for which all primes lower than it are known. In this case not all primes lower than it are known.

Unfortunately it is not possible to, for any \(n\) know that it is the next prime without testing the range between n and the previous prime. Therefore, for now, this part of it is the only way to check if previous elements are prime.

(< @max-prime max-div)
(do (doall (map #(is-prime? % true) (range @max-prime max-div)))
    (if (div? n)
      (do (swap! found-primes conj n)
          true)
      false))

Now we check if a number is or is not in fact prime, and in this case we do set the prime in the event it is contiguous because the number is the next number in the set of found primes.

(div? n)
(do (if contiguous?
      (do (swap! found-primes conj n)
          (reset! max-prime n))
      (swap! found-primes conj n))
    true)

Once we have this function, we can then begin to look into generating the primes to some value. Like before we want to ensure that we can use as many numbers form our prime list as possible, therefore we implement a simple function that will first attempt to return all the values from our found-primes list. should they not exist, it will then use the side effects of the is-prime? function to generate it.

<<divisible?>>

(defn prime-range
  ([end] (prime-range 2 end))
  ([start end]

   (if (@max-prime)
     (drop-while))
   
   (->> (range (max start 2) end)
        (filter is-prime?))))

(defn prime-factors [x]
  (->> (prime-range (inc (Math/sqrt x)))
       (filter #(divisible? x %))))


(last (prime-factors 600851475143))