Project Euler Problem 12

Highly Divisible Triangular Number

Project Euler presents the following Problem #12:

The sequence of triangle numbers is generated by adding the natural numbers. So the 7$^{th}$ 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,\,\dots$$

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

$$\begin{align*} \bold{1}&:\,1\\ \bold{3}&:\,1,\,3\\ \bold{6}&:\,1,\,2,\,3,\,6\\ \bold{10}&:\,1,\,2,\,5,\,10\\ \bold{15}&:\,1,\,3,\,5,\,15\\ \bold{21}&:\,1,\,3,\,7,\,21\\ \bold{28}&:\,1,\,2,\,4,\,7,\,14,\,28 \end{align*}$$

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?

Helper Functions

As mentioned above, and seen as early as Project Euler Problem #1, triangle numbers are generated using the formula:

1/**
2 * Calculates the sum of the first [this] natural numbers, based on the formula:
3 *
4 *      {n}Sigma{k=1} k = n * (n + 1) / 2
5 *
6 * Conversion of very large Doubles to Longs in this formula can lead to large rounding
7 * losses, so Integer division by 2 is replaced with a single bitwise right shift,
8 * as n >> 1 = n / (2^1).
9 */
10fun Int.gaussSum(): Long  = 1L * this * (this + 1) shr 1

While previously solving Problem #3, we implemented a function that returns the prime decomposition of a number:

1/**
2 * Prime decomposition repeatedly divides out all prime factors using a Direct Search
3 * Factorization algorithm without any optimization.
4 *
5 * @return map of prime factors (keys) and their exponents (values).
6 * @throws IllegalArgumentException if [n] <= 1.
7 */
8fun primeFactors(n: Long): Map<Long, Int> {
9    require(n > 1) { "Must provide a natural number greater than 1" }
10    var num = n
11    val primes = mutableMapOf<Long, Int>()
12    var factor = 2L
13    while (factor * factor <= num) {
14        while (num % factor == 0L && num != factor) {
15            primes[factor] = primes.getOrDefault(factor, 0) + 1
16            num /= factor
17        }
18        factor++
19    }
20    if (num > 1) primes[num] = primes.getOrDefault(num, 0) + 1
21    return primes
22}

Once a triangle number $t$ is generated and converted to its decomposed form, the exponents $e$ of each prime factor $p$ can be used to count the total number of divisors $d(t)$, as shown here:

1fun countDivisors(n: Int): Int {
2    return primeFactors(1L * n).values
3        .map { it + 1 }
4        .reduce { acc, v -> acc * v }
5}

With the help of these three functions, we can compare a few different approaches to brute forcing a solution.

Brute Force

The first approach consists of generating each triangle number and counting its divisors until a count is found that exceeds the limit, $N$:

1fun firstTriangleBruteA(limit: Int): Int {
2    var t = 3  // first triangle number with >1 divisor
3    if (limit == 1) return t
4    var n = 2
5    var count = 2
6    while (count <= limit) {
7        n++
8        t = n.gaussSum().toInt()
9        count = countDivisors(t)
10    }
11    return t
12}

A second approach avoids using the function gaussSum() by taking advantage of Kotlin's sequences. Since we don't need the result of the entire generated triangle number progression, we can use a sequence to lazily process each element until the first to match the predicate is found:

1fun firstTriangleBruteB(limit: Int): Int {
2    return generateSequence(Pair(3, 3)) { Pair(it.first + it.second, it.second + 1) }
3        .first { countDivisors(it.first) > limit }
4        .first
5}

Line 2: The sequence of triangle numbers is generated by cumulatively adding each natural number, using a starting value seed, with the result and the next natural number stored as a Pair. For further clarification, here's an example of generateSequence() being used to return a Fibonacci sequence.

A third brute force solution avoids counting the number of divisors for very large triangle numbers by instead focusing on the following observation: given that $n$ and $n+1$ must be coprime, each triangle number can be reduced to two smaller factors. This is based on the equation:

The count of divisors is calculated for each smaller factor by cycling through a single new value dependent on the polarity of the current $n$:

1fun firstTriangleBruteC(limit: Int): Int {
2    if (limit == 1) return 3
3    var n = 2 // D(2) = D(1) * D(3)
4    var isEven = true
5    var dn1 = 2 // D(3) = 2
6    var count = 2
7    while (count <= limit) {
8        n++
9        isEven = !isEven
10        val dn2 = if (isEven) countDivisors(n + 1) else countDivisors((n + 1) / 2)
11        count = dn1 * dn2
12        dn1 = dn2
13    }
14    return n.gaussSum().toInt()
15}

Here's a breakdown of the average execution speeds on my system when $N=1000$:

It should be possible to appreciate that, while it improves on the first two solutions, firstTriangleBruteC() still has room for improvement as countDivisors() is occasionally called for repeated values. This also becomes more evident if the solution is used to ascertain the results of multiple values of $N$.

Solution With Two Caches

To reduce how many times countDivisors() is invoked, let's introduce a cache that stores the count of divisors.

The issue with implementing this is that we don't know the upper size limit to enforce. If we're interested in getting the first triangle number to have over 500 divisors, what is the first natural number to have over 500 divisors? Are these numbers even maybe equivalent? Since we don't have these answers, let's just choose an arbitrary value, 20,000 for now, and implement a fallback in the event that this value is exceeded.

⚠️ 20,000 is in fact too low and the fallback function will need to be called over 11,000 times.

1fun firstTrianglesCache(limit: Int) {
2    val eMax = 20_000
3    val divisorCount = IntArray(eMax) { 2 }.apply { this[1] = 1 }
4    for (divisor in 2 until eMax) {
5        for (num in divisor * 2 until eMax step divisor) {
6            divisorCount[num]++
7        }
8    }
9    // TODO
10}

Line 3: The array is initialized with all values equal to 2, since natural numbers will have at minimum two divisors (the number 1 and the number itself). The exception, 1, is adjusted appropriately and the value at index 0 is ignored as it will never be accessed.

Lines 4-8: Iterates over every number and increments the count for every multiple below the arbitrary limit, except the number multiplied by 1 (as this was already handled during initialization).

Since this solution is about caching, let's also store the first triangle number to exceed $n$ divisors for every $n<=N$, so that the returned array can be quickly index-accessed for multiple values:

1fun firstTrianglesCache(limit: Int): IntArray {
2    // ✂️ caching of divisor counts
3    val triangles = IntArray(limit + 1).apply {
4        this[0] = 1
5        this[1] = 3
6        this[2] = 6
7    }
8
9    var lastN = 3
10    var isEven = false
11    var lastDn1 = 2
12    var lastTriangle = 6
13    var lastCount = 4
14    for (i in 3..limit) {
15        if (i < lastCount) {
16            triangles[i] = lastTriangle
17            continue
18        }
19        // TODO
20    }
21    return triangles
22}

Line 3: Each index represents a divisor count $n$. For example, the first triangle number to have more than 2 divisors is 6 and is stored at index 2. Note that 6 is the 3rd triangle number, so lastN = 3.

Line 15: The only value of interest is the smallest triangle number to have a divisor count that exceeds each index value, so this conditional block avoids a lot of repeated counts. For example, the triangle number 120 is the correct result for 7 consecutive divisor counts.

From here on, the rest of the solution is similar to firstTriangleBruteC() except that the triangles cache and local variables are adjusted every time a new triangle number is encountered with a divisor count that exceeds the current index value. Here's the final solution:

1fun firstTrianglesCache(limit: Int): IntArray {
2    val eMax = 20_000
3    val divisorCount = IntArray(eMax) { 2 }.apply { this[1] = 1 }
4    for (divisor in 2 until eMax) {
5        for (num in divisor * 2 until eMax step divisor) {
6            divisorCount[num]++
7        }
8    }
9    val triangles = IntArray(limit + 1).apply {
10        this[0] = 1
11        this[1] = 3
12        this[2] = 6
13    }
14
15    var lastN = 3
16    var isEven = false
17    var lastDn1 = 2
18    var lastTriangle = 6
19    var lastCount = 4
20    for (i in 3..limit) {
21        if (i < lastCount) {
22            triangles[i] = lastTriangle
23            continue
24        }
25        var nextN = lastN + 1
26        isEven = !isEven
27        do {
28            val toCount = if (isEven) nextN + 1 else (nextN + 1) / 2
29            val dn2 = if (toCount < eMax) divisorCount[toCount] else countDivisors(toCount)
30            val count = lastDn1 * dn2
31            lastDn1 = dn2
32            if (i < count) {
33                val triangle = nextN.gaussSum().toInt()
34                triangles[i] = triangle
35                lastTriangle = triangle
36                lastN = nextN
37                lastCount = count
38                break
39            }
40            nextN++
41            isEven = !isEven
42        } while (true)
43    }
44    return triangles
45}

Line 29: This fallback checks if a previously determined divisor count could exist based on the chosen upper limit. If not, the helper function is invoked.

📊 This solution returns a quick-pick array of all solutions for $N<=1000$ in approximately 94.44ms on my system.

Solution With One Cache

Taking a closer look at firstTrianglesCache() shows that the majority of local variables exist for the sole purpose of generating the returned triangles cache. If we return to the original problem goal involving a single result, we can eliminate much of this code and bring it more in line with the third brute force solution, firstTrianglesBruteC(). The difference being that the divisor count will be both cached and calculated with every iteration.

Once again, an arbitrary size has to be chosen for the cache, but by now we've worked through enough solutions to exhaust all values of $N<=1000$ and can be certain that 41,100 will suffice. This limit is actually overkill for smaller values of $N$ as brute forcing all the ratios of triangle number to $N$ shows that the maximum index accessed will not exceed a factor of 53.

1fun firstTriangle(n: Int): Int {
2    val nMax = minOf(n * 53, 41100)
3    val divisorCount = IntArray(nMax)
4    var num = 0
5    var dT = 0
6    while (dT <= n) {
7        num++
8        for (i in num until nMax step num) {
9            divisorCount[i]++
10        }
11        dT = if (num % 2 == 0) {
12            divisorCount[num/2] * divisorCount[num-1]
13        } else {
14            divisorCount[num] * divisorCount[(num-1)/2]
15        }
16    }
17    return num * (num - 1) / 2
18}

This alternate solution also uses the cyclic formula from the previously detailed $equation (1)$, but looking forward to $n+1$ is not possible since the cache is only filled up to $n$ (and its multiples) with every iteration. Instead, iteration starts with a count of 0 divisors from the $0^{th}$ triangle number and treats that result as the first result, thereby allowing us to look backwards using $n-1$. For example, when num = 4, $d(4)=d(2)*d(3)=2*2=4$, which maps to the 3rd triangle number, 6.

📊 The average execution speed of this solution performs better than all the brute force solutions at 27.62ms when $N=1000$.

📚 While the triangle number sequence was mentioned explicitly in the problem outline, some Project Euler problems leave it up to the solver to deduce any potential patterns. If you are curious about whether some observed numbers are part of a larger documented sequence, consider performing a search using the very helpful On-Line Encyclopedia of Integer Sequences.

Here's the full source code for this problem (and its corresponding test suite).

This solution set is also available in Python and in C++ on my GitHub.

Problem content on the Project Euler site is licensed under a Creative Commons (CC) License: CC BY-NC-SA 4.0. Solutions to only the first hundred problems are encouraged to be distributed publicly.

Thanks for reading!🏝️🔺📐🚩