Problem 012: Highly divisible triangle number

From Project Euler:

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Let us list the factors of the first seven triangle numbers:

1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28
We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?

A way to solve this problem is fairly straight-forward: progress through all the triangle numbers, figure out their divisors, and count them. A simple brute force method is to just enumerate all the numbers up to $\sqrt{n}$, and for each one that divides $n$, add 2 to our count (one for this number, and one for $n$ divided by it). The corner case is if $n$ is square, just count the square root once:

(defn num-divisors [n]
  (loop [current 1
         total 0]
    (cond
      (> current (Math/sqrt n)) total
      (== current (Math/sqrt n)) (+ 1 total)
      :else (recur
              (+ 1 current)
              (+ total
                 (if (== 0 (mod n current)) 2 0))))))

We can then do a simple search through all the triangle numbers to get the first one that has at least 500 divisors:

(defn search-brute-force [limit]
  (loop [n 1
         triangle 1]
    (let [d (num-divisors n)]
     (if (> d limit)
       triangle
       (recur (+ 1 n)
              (+ triangle (+ 1 n)))))))

Unfortunately, this goes well past the one minute limit. Let's see if we can speed this up a bit. Take the formula for a triangle number $T_n$:
$$T_n = \sum_{i=1}^n = \frac{n(n + 1)}{2}$$
From this formula we can infer a few other interesting details. Whenever you have two sequential numbers one of them will be even and the other will be odd. Any even number will be divisible by 2, so the 2 in the denominator will get cancelled out by that and we'll be left with a divisor of $T_n$. The odd number left over will also be a divisor.

To write this out more explicitly, if $n$ is even then $\frac{n}{2}$ and $n + 1$ are divisors of $T_n$; if $n$ is odd then $n$ and $\frac{n + 1}{2}$ are divisors of $T_n$.

Next we'll want to dig into the properties of divisors. Any number can be broken down into a product of prime numbers, and some of those primes might be duplicated (for example, 4 = 2 * 2). So we can write any number $n$ like this:
$$n = p_1^{c_1} * p_2^{c_2} * ...$$
Then for each prime $p_i$, there's going to be $c_i$ factors of $n$ from that prime: $p_i$, $p_i^2$, $p_i^3$, all the way up to $c_i$. Each one of those can multiply with every other factor of $n$ to produce another factor (for example with $n = 28$ you can have $2 * 7$ to produce 14). So the number of divisors $d(n)$ is $c_i + 1$ multiplied by all the other possibilities:
$$d(n) = (c_1 + 1) * (c_2 + 1) * ...$$
Lastly, if you have two numbers $a$ and $b$ that are co-prime with one another, this means that they do not share any divisors other than 1. Therefore none of the $p_i$s are the same, so the number of divisors is just the product of the number of divisors for $a$ and $b$ separately. Since $n$ and $n + 1$ are necessarily coprime, we can just multiply the number of divisors of $n$ with the number of divisors of $n + 1$, dividing the even one by 2 to remove the denominator in the formula for $T_n$:
$$d(T_n) = \left\{ \begin{array}{ll} d(\frac{n}{2}) * d(n + 1) & n\ even \\ d(n) * d(\frac{n + 1}{2}) & n\ odd \end{array} \right.$$
These two numbers should be quite a lot smaller than $T_n$, so the search time will be much faster.

Let's write out a little function to calculate the number of divisors of a triangle number:

(defn num-triangle-divisors [n]
  (if (even? n)
    (* (num-divisors (quot n 2))
       (num-divisors (+ 1 n)))
    (* (num-divisors n)
       (num-divisors (quot (+ 1 n) 2)))))

And then tweak our search function to use this instead:

(defn search-triangle [limit]
  (loop [n 1
         triangle 1]
    (let [d (num-triangle-divisors n)]
     (if (> d limit)
       triangle
       (recur (+ 1 n)
              (+ triangle (+ 1 n)))))))

This one is much faster, even with the 3s or so it takes for lein to start up:

$ time lein run
76576500

real 0m3.848s
user 0m4.687s
sys 0m0.355s

Problem 011: Largest product in a grid

From Project Euler:

In the 20×20 grid below, four numbers along a diagonal line have been marked in red.

08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48

The product of these numbers is 26 × 63 × 78 × 14 = 1788696.

What is the greatest product of four adjacent numbers in the same direction (up, down, left, right, or diagonally) in the 20×20 grid?

This one doesn't really require any trickery, and can even be done by hand. That goes against our spirit of writing code to solve everything though, so here's a solution. We scan through each row except for the last 3, collecting the products in four directions: horizontally, vertically, and the diagonal directions left and right. Then we just take the max of all these products:

(def Grid
  [[ 8  2 22 97 38 15  0 40  0 75  4  5  7 78 52 12 50 77 91  8]
   [49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48  4 56 62  0]
   [81 49 31 73 55 79 14 29 93 71 40 67 53 88 30  3 49 13 36 65]
   [52 70 95 23  4 60 11 42 69 24 68 56  1 32 56 71 37  2 36 91]
   [22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80]
   [24 47 32 60 99  3 45  2 44 75 33 53 78 36 84 20 35 17 12 50]
   [32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70]
   [67 26 20 68  2 62 12 20 95 63 94 39 63  8 40 91 66 49 94 21]
   [24 55 58  5 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72]
   [21 36 23  9 75  0 76 44 20 45 35 14  0 61 33 97 34 31 33 95]
   [78 17 53 28 22 75 31 67 15 94  3 80  4 62 16 14  9 53 56 92]
   [16 39  5 42 96 35 31 47 55 58 88 24  0 17 54 24 36 29 85 57]
   [86 56  0 48 35 71 89  7  5 44 44 37 44 60 21 58 51 54 17 58]
   [19 80 81 68  5 94 47 69 28 73 92 13 86 52 17 77  4 89 55 40]
   [ 4 52  8 83 97 35 99 16  7 97 57 32 16 26 26 79 33 27 98 66]
   [88 36 68 87 57 62 20 72  3 46 33 67 46 55 12 32 63 93 53 69]
   [ 4 42 16 73 38 25 39 11 24 94 72 18  8 46 29 32 40 62 76 36]
   [20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74  4 36 16]
   [20 73 35 29 78 31 90  1 74 31 49 71 48 86 81 16 23 57  5 54]
   [ 1 70 54 71 83 51 54 69 16 92 33 48 61 43 52  1 89 19 67 48]])

(def N (count Grid))

(defn search []
  (apply max
    (for [x (range N)
          y (range (- N 4))
          :let
          [horiz (if (< x (- N 4))
                    (->> (range 4)
                         (map #(nth (nth Grid y) (+ x %)))
                         (reduce *))
                    0)
           vert (if (< y (- N 4))
                   (->> (range 4)
                        (map #(nth (nth Grid (+ y %)) x))
                        (reduce *))
                   0)
           ldiag (if (and (>= x 3)
                           (< y (- N 4)))
                    (->> (range 4)
                         (map #(nth (nth Grid (+ y %)) (- x %)))
                         (reduce *))
                   0)
           rdiag (if (and (< x (- N 4))
                           (< y (- N 4)))
                    (->> (range 4)
                         (map #(nth (nth Grid (+ y %)) (+ x %)))
                         (reduce *))
                   0)]]
      (max horiz vert ldiag rdiag))))

This runs and gives the correct result (takes about 5s for lein/JVM to get ready):

$ time lein run
70600674

real 0m5.800s
user 0m7.340s
sys 0m0.404s

Problem 010: Summation of primes

From Project Euler:

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.

My initial approach to this was to do something using the Sieve of Eratosthenes, and my solution was something like this:

(defn prime-sieve-sum [n]
  (loop [nums (vec (concat [2] (range 3 n 2)))
         total 0]
    (if (empty? nums)
      total
      (let [current (first nums)]
        (recur (vec (remove #(== 0 (mod % current)) (rest nums)))
               (+ current total))))))

This keeps a vector of all the possible primes, and filters out any multiples as it goes along. The downside is that this was very slow: it took about 10 minutes to run on my machine with $n$ at 2 million.

My guess is that since it is constantly creating new vectors, there is a lot of memory copying going on. A better approach might be to keep a large vector of bools and set them to false as each multiple is encountered. This doesn't jive well with Clojure's immutable vectors, but you can fall back to Java's ArrayList if you want to do this approach.

Instead I just decided to brute force it. Using the is-prime? function from Problem 003, I just scanned across:

(defn prime-sum [n]
  (loop [current 3
         total 2]
    (if (>= current n)
      total
      (recur (+ 2 current)
             (+ total (if (is-prime? current)
                        current
                        0))))))

This gets the job done much more quickly (about 5s overhead for lein/JVM):

$ time lein run
Running...
142913828922

real 0m21.508s
user 0m22.936s
sys 0m0.530s

Problem 009: Special Pythagorean Triplet

From Project Euler:

A Pythagorean triplet is a set of three natural numbers, $a < b < c$, for which, $$a^2 + b^2 = c^2$$ For example, 32 + 42 = 9 + 16 = 25 = 52. There exists exactly one Pythagorean triplet for which $a + b + c = 1000$. Find the product $abc$.

For this one, we could go ahead and dive into properties of Pythagorean triplets, of which there are many. You can read the Wikipedia page about them and get many hints on how to do this problem more efficiently than what I'm going to show you, and in fact in later articles we'll want to take advantage of these.

However as we say in software engineering: if the sub-optimal solution is good enough and more clear than more efficient solutions, then consider just using the sub-optimal solution. In this case, we can just use Clojure's for comprehension to give us an answer quickly:

(defn find-triplet [s]
  (first
    (for [a (range 1 s)
          b (range (+ a 1) (- s a))
          :let [c (- s a b)]
          :when (== (* c c)
                    (+ (* a a) (* b b)))]
      (list a b c))))

(println (reduce * (find-triplet 1000)))

This runs fast enough (remember that the JVM takes about 1.5s to start up on my machine):

$ time clj 9.clj
31875000

real 0m1.635s
user 0m3.341s
sys 0m0.084s

Problem 008: Largest product in a series

From Project Euler:

The four adjacent digits in the 1000-digit number that have the greatest product are 9 × 9 × 8 × 9 = 5832.

73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450

Find the thirteen adjacent digits in the 1000-digit number that have the greatest product. What is the value of this product?

This one was pretty straightforward. You could technically do it using string manipulation, but it's easy and fast to do it with numeric manipulation. First step was just to copy-paste the giant number into a file, and load it into a constant:

(def N
  (read-string (clojure.string/replace (slurp "big_number.txt") "\n" "")))

I also wrote a few helper functions. Some of these are available through libraries, but setting up a Leiningen project for everything is a bit tedious for functions that are small:

(defn pow [a b] (reduce * 1 (repeat b a)))

(defn prod [n]
  "Multiply the digits of `n` together."
  (if (< n 10)
    n
    (* (mod n 10) (prod (quot n 10)))))

Once you do that, you can just loop across $N$, chopping off the last $k$ digits and seeing if the product of those digits is bigger than anything you've seen before:

(defn top [_n num-digits]
  (let [mask (pow 10 num-digits)]
    (loop [n _n best 0]
      (if (== n 0)
        best
        (recur (quot n 10)
               (max (prod (mod n mask))
                    best))))))

(println (top N 13))

This runs pretty quick (remember it takes ~1.5s for the JVM to start up):

$ time clj 8.clj 
23514624000N

real 0m1.569s
user 0m3.196s
sys 0m0.111s

Problem 007: 10001st prime

From Project Euler:

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

What is the 10 001st prime number?

My initial approach to this was to use a Sieve of Eratosthenes to efficiently find all the prime numbers up to a given limit, however that ended up spending a lot of processing in manipulating lists. Turns out there are faster ways to do it.

I also discovered that Clojure does not have tail-call recursion, you need to use the loop and recur constructs to avoid getting stack overflows.

There are two optimizations here. The first is that we'll skip the first prime 2, since it's obviously not the 10001st prime and it is an anomaly that doesn't let us use the next optimization, which is that all primes greater than 3 have the form $6k \pm 1$. Instead of iterating through each number or even each odd number, we'll iterate through all numbers of that form.

Here's the code:

(defn nth-prime [n]
  (loop [i (- n 2)
         k 1
         s -1]
    (let [p (+ (* 6 k) s)
          next-s (if (= s 1) -1 1)
          next-k (if (= s 1) (+ k 1) k)]
      (cond
        (not (is-prime? p)) (recur i next-k next-s)
        (= i 0) p
        :else (recur (- i 1) next-k next-s)))))

Running this code is actually surprisingly slow:

rob@alien ~/code/euler $ time clj 7.clj 
104743

real 0m16.076s
user 0m17.930s
sys 0m0.349s

I'll have to do some profiling to figure out why this is the case. When translating the code as-is to Python it takes 0.176s, and Python isn't exactly known for speed.

Problem 006: Sum Square Difference

From Project Euler:

The sum of the squares of the first ten natural numbers is,

12 + 22 + ... + 102 = 385
The square of the sum of the first ten natural numbers is,

(1 + 2 + ... + 10)2 = 552 = 3025
Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025 − 385 = 2640.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

The brute-force algorithm for this is pretty simple, just run a loop and calculate the result:

(def N 101)

(defn square [x] (* x x))

(println (- (square (reduce + (range 1 N)))
            (->> (range 1 N)
                 (map square)
                 (reduce +))))

This runs pretty fast (note that the JVM takes 1.5s to start up on my machine):

$ time clj 6.clj
25164150

real 0m1.517s
user 0m3.001s
sys 0m0.133s

We can do better though. As we saw in Problem 001, we can use a formula to calculate the sum of numbers from 1 to $N$:
$$\sum^n_{i=1} i = \frac{n(n+1)}{2}$$
Or in Clojure:

(defn sumto [n] (quot (* n (+ n 1)) 2))

If we square this value, we get the square of the sums (the first value in our subtraction). This is O(1) rather than O(n).

Turns out there is also a formula for the sum of squares of a sequence:
$$\sum^n_{i=1} i^2 = \frac{n(n+1)(2n+1)}{6}$$
For more details, see here.

Writing this out in Clojure:

(defn sumto-sq [n] (quot (* n
                            (+ n 1)
                            (+ (* 2 n) 1))
                         6))

Combining them all together:

(def N 100)
(println (- (square (sumto N))
            (sumto-sq N)))

This also runs fast, and should run fast for any number that fits into a 64-bit integer¹:

$ time clj 6.clj
25164150

real 0m1.583s
user 0m3.043s
sys 0m0.140s

---

1 We may have to algebraically manipulate the formula a bit to avoid overflows but it should be doable.^↩