Post date: Mar 17, 2017 7:17:43 PM
The second Project Euler problem asks us to find the sum of the even Fibonacci numbers less than 4000000. As such, the simple generate, filter, and sum solution will look like that of Problem 1, but with a couple twists:
(defn euler-2 [] (->> (fib-seq) (take-while #(< % 4000000)) (filter even?) (reduce +)))
Here, we take from the sequence of Fibonacci numbers in order until the next Fibonacci number is greater than or equal to 4000000, filter the resulting collection to get the even elements, and sum them by reducing - a straightforward translation of the problem. The main trick would seem to be writing the function fib-seq. Our first attempt is as follows:
(defn fib-seq-bad "Returns a sequence of Fibonacci numbers." ([] (fib-seq 0N 1N)) ([a b] (cons b (fib-seq-bad b (+ a b)))))
In this code, a represents the previous Fibonacci number and b represents the current Fibonacci number in the sequence. Assume we're at Fk in the sequence. The next time we call fib-seq-bad, Fk, the current value, will be the next previous value, and Fk+1=Fk+Fk-1(the current previous value)=a+b will be the next value. This gives us the recursion relation shown in the function. The zero arity function seeds the 2 arity function with initial values of 0 and 1 (previous and current values for the first time the function is called). The list is built by consing each current value to the list generated by the next values. The only problem comes when we actually try to run the function - it fails with a stack overflow error. The reason why is inherent in the solution - in order to cons the current Fibonacci number onto the tail of the sequence, you first need to build that sequence. And to build this infinite sequence, you must recurse infinitely. The good news about Clojure is that instead of forming this infinite list, we can instead use lazy-seq to form a lazy sequence.
So what is lazy-seq?
According to its Clojure documentation, lazy-seq:
Takes a body of expressions that returns an ISeq or nil, and yields a Seqable object that will invoke the body only the first time seq is called, and will cache the result and return it on all subsequent seq calls.
In our case, we'd use it like this:
(defn fib-seq "Returns a lazy sequence of Fibonacci numbers" ([] (fib-seq 0N 1N)) ([a b] (lazy-seq (cons b (fib-seq b (+ a b))))))
In this code, the first time seq is called to get the contents of the list, the body is evaluated, leading to a list holding the current Fibonacci number, followed by the rest of the list. But when we attempt to build the rest of the list, our call to fib-seq immediately hits a lazy-seq call, which does not get invoked until the next seq is called. Each time this happens, another item is consed onto the end of the list, allowing us to realize the list in order and only as needed. As such, we've created a lazy sequence of Fibonacci numbers which can be processed by our solution. This is not the only lazy sequence we shall be building. Problem 3 also utilizes lazy-seq to generate a sequence of prime numbers.
In this problem, a key piece is generating a lazy sequence of primes. To do this, we use a sieve - in a map we hold associations between the next number that is going to be sieved out by the current set of primes and the primes that number will be sieved out by. To find the next prime, we look in the table at the next number, d. If we find d in the table, we know it is not prime and that the next number(s) sieved by the algorithm will be the sum of that number and the prime(s) that sieves the number. If d's not in the table, then it's prime. We'll step through the algorithm to demonstrate.
Start with an empty map and initial number d=2. Since 2 is not a key in the map, 2 must not be sieved yet and is thus a prime - add it to the list of primes and add to the map [d2, (d)]. The map now holds {4 (2)}. The next number is 3. We check and it is not in the map - again this number is prime so add it to the list of primes. Add [9 (3)] to the map, giving a map of {4 (2), 9 (3)}. The next number is 4. It is in the map so it is not prime. For each of the sieve factors of d, add the sieve factor for the number and update the map. This gives us a map of {6 (2), 9 (3)}. The next number, 5, is not in the map and prime, leading to the resulting map {6 (2), 9 (3), 25 (5)}. Six is in the map and is not prime. We increment it in the sieve, resulting in the map {8 (2), 9 (3), 25 (5)}. We continue this way until we come to the number 10, when our map contains {10 (2), 12 (3), 25 (5), 49 (7) }. Ten is not prime, but when we add it's sieving factor to it, we find that the number 12 is already in the map. To fix this, we add both sieve factors under 12, giving the map {12 (2 3), 25 (5), 49 (7)}. Eleven is prime so we augment the map again giving {12 (2 3), 25 (5), 49 (7), 121 (11)}. Twelve, the next number, is not prime and we increment the number by both factors in the sieve, re-splitting the factors giving: {14 (2), 15 (3), 25 (5), 49 (7), 121 (11)}. We continue, adding primes to our list and updating the sieve properly depending on whether the currently examined number is prime or not.. The algorithm is shown here:
(defn gen-primes "Generates an infinite, lazy sequence of prime numbers" [] (let [reinsert (fn [table x prime] (update-in table [(+ prime x)] conj prime))] (defn primes-step [table d] (if-let [factors (get table d)] (recur (reduce #(reinsert %1 d %2) (dissoc table d) factors) (inc d)) (lazy-seq (cons d (primes-step (assoc table (* d d) (list d)) (inc d)))))) (primes-step {} 2)))
In this code, reinsert does the job of updating the table with the sieved numbers and their sieving primes. The call to recur (when the number is in the table) handles the case where the number is not prime using reinsert to update the table (repeatedly for each prime sieving factor, using reduce), while the cons branch of the code handles the case when the current number is prime, adding of the current number being checked to the list of primes and updating the sieve with the square of the prime. The lazy-seq call in the cons branch ensures that the list will be delivered lazily. And, of course, we start at 2 with an empty map.
Why is the first sievable number for a given prime initialized to that prime's square? One can see that earlier multiples of this prime will be caught by composites and primes less than this number. So 2*5 will be caught by 2's sieve, 3*5 by 3's etc. The first number not caught by a number's sieve is that number squared (and after that, its increments take care of the sieving).
So what does the problem ultimately ask? It wants the maximal prime factor of 600851475143. How do we find the maximal factor? Assume a number N can be written as a product of prime numbers p0ip1j...pnk, with factors ordered from low to high. By dividing by p0, we obtain a smaller number. Continuing removing smallest factors in this manner, we are eventually left with N=pk, the maximal prime factor. To find the smallest prime factor of the number at any iteration, test against the sequence of primes until the first one is found that divides N. This will be the smallest prime factor.
The code to accomplish all this is:
(defn smallest-prime-factor "Iterate through the sequence of primes until you find a divisor." [N] (first (drop-while #(not (factor? % N)) (gen-primes)))) (defn max-prime-factor "Divide by smaller factors until you're left with the largest." [N] (let [p (smallest-prime-factor N)] (if (= N p) p (recur (/ N p)))))
This allows the solution to Problem 3 to be written simply as:
(defn euler-3 [] (max-prime-factor 600851475143))
This problem (like several from the Project Euler corpus) deals with palindromic numbers (numbers whose digits read the same backward and forward). A prime requirement is to test whether a given number is a palindrome. For this, we check if the reverse of the number's string equals the number's string:
(defn palindrome? [N] (= (str N) (apply str (reverse (str N)))))
The problem asks for the maximal palindromic number formed from the product of two three digit numbers. The set of three digit numbers can be generated by:
(def three-digits (range 100 1000))
This returns the numbers from 100 to 999. The remainder of the problem can be written as follows:
(defn euler-4 [] (->> (combo/cartesian-product three-digits three-digits) (filter (fn [[i j]] #(>= i j))) (map (fn [[i j]] (* i j))) (filter palindrome?) (apply max)))
In it, we generate all pairs of three digit numbers using combo/cartesian-product. We take those pairs and, since i*j=j*i, we get rid of half of them by trying only the ones where i>=j. Next we map these pairs to their products and filter them on whether the products are palindromic. Finally, we take the max of all of the palindromic products to find the answer.
Next post, we'll tackle problems five through nine.