Project Euler Problem 9

Sep 29th, 2022

11 min read

Special Pythagorean Triplet

Project Euler presents the following Problem #9:

A Pythagorean triplet is a set of three natural numbers, $a < b < c$ , for which,
$a^2+b^2=c^2$
For example, $3^2+4^2=9+16=25=5^2$ .

There exists exactly one Pythagorean triplet for which $a+b+c=1000$ . Find the product $abc$ .

Definitions
Pythagorean Triplet: Refers to a set of three positive integers $(a,b,c)$ , such that a right-angled triangle exists with the largest value $c$ as the hypotenuse, thereby satisfying the equation above. The first five triplets are: (3, 4, 5), (6, 8, 10), (5, 12, 13), (9, 12, 15), (8, 15, 17). The term "triplet" is used interchangeably with "triple" depending on the source.
Primitive Pythagorean Triplet: Refers to a Pythagorean triplet, such that $\gcd(a,b,c)=1$ . The first five primitive triplets are: (3, 4, 5), (5, 12, 13), (8, 15, 17), (7, 24, 25), (20, 21, 29). All non-primitive triplets can be generated from a primitive one.
Coprime: Two integers are considered coprime if their only common positive divisor is 1. This property is also referred to as "relatively prime" or "mutually prime".

Pythagorean triplets have special properties (many that we won't cover):

If primitive, one of $a$ or $b$ must be even and the other odd, with $c$ always being odd. This means that all generated triplets will switch between having all their elements be even (from even factors) or all their elements following the same pattern as the original primitive (from odd factors). For example, (3, 4, 5) $\to$ (6, 8, 10) $\to$ (9, 12, 15) $\to$ (12, 16, 20) $\to\dots$
Property #1 then dictates that the perimeter $P$ of a Pythagorean triplet triangle must always be even, as the sum of all even numbers is also an even number, as is the sum of two odd numbers.
At least one element of every triplet is divisible by 3, 4, or 5, although an element may not be divisible by any of those factors.
The product of the two non-hypotenuse sides of a Pythagorean triplet triangle is always divisible by 12 and the product $abc$ is always divisible by 60.

Let's use a few of these properties to implement three functions that take a perimeter $P$ and return the largest product $abc$ if a triplet is found whose elements sum to $P$ ; otherwise, the functions will return -1. Note that, while the original problem makes it clear that only one triplet exists for $P=1000$ , triplets that share the same $P$ do exist.

Nested Loop Solution

This first solution iterates through combinations of $(a,b,c)$ with an outer loop representing values of $c$ and an inner loop representing $b$ , while using some of the above properties to limit the amount of considered triplets.

First, any value $P$ that is not even (property #2) is immediately rejected as it will not have a valid triplet.

Second, $c$ is constrained to being greater than $\dfrac{P}{3}$ , but not greater than $\dfrac{P}{2}-1$ . Its lower limit is based on the following:

\begin{align*} c>b>a&\implies c>b\;\And\;c>a\\[.5em] \therefore\; 2c&>b+a\\[.5em] 2c+c&>b+a+c\\[.5em] 3c&>P\\[.5em] c&>\dfrac{P}{3}\tag{1} \end{align*}

The upper limit of $c$ can be shown based on the following:

\begin{align*} \text{assuming }\;c &< \dfrac{P}{2}\\[.5em] 2c &< a+b+c\\[.5em] \therefore\; c &< a+b\\[.5em] \text{given}\;c&=\sqrt{a^2+b^2}\\[.5em] \sqrt{a^2+b^2} &< a+b\\[.5em] a^2+b^2 &< a^2+2ab+b^2\\[.5em] 0 &< 2ab \end{align*}

Third, $b$ must of course be smaller than $c$ , but should not drop below $\dfrac{P-c}{2}$ .

Fourth, $b$ must be greater than $a$ , so iterating backwards through these two loops allows us to break out of the inner loop early and start fresh with a new $c$ if $b$ ends up being too small in value.

Once a triplet has reached this point, it only needs to satisfy two scenarios: property #4 and the Pythagoras Theorem. hypot() from the kotlin.math package solves the formula $c=\sqrt{a^2+b^2}$ to validate the triplet as being Pythagorean. Once a valid triplet is found, there is no point in iterating over further values of $b$ as any future combinations would result in a product of lesser value.

1import kotlin.math.hypot
2
3fun pythagoreanTripletBC(p: Int): Int {
4    if (p % 2 != 0) return -1
5    
6    var product = 0
7    nextC@for (c in (p / 2 - 1) downTo (p / 3 + 1)) {
8        val diff = p - c
9        nextB@for (b in (c - 1) downTo (diff / 2)) {
10            val a = diff - b
11            if (b <= a) continue@nextC
12            if (a * b % 12 != 0) continue@nextB
13            if (hypot(a.toDouble(), b.toDouble()) == c.toDouble()) {
14                product = maxOf(product, a * b * c)
15                break
16            }
17        }
18    }
19
20    return product
21}

📊 pythagoreanTripletBC() checked 51,758 triplets by invoking hypot() when p = 3000. Let's see how a single loop solution can reduce this amount in the next section.

Single Loop Solution

The nested loop solution can be altered to use only one for loop by focusing on $a$ instead of $c$ .

As before, only even values of $P$ are considered. Then, given that (3,4,5) is the smallest existing triplet, $a$ must be greater than or equal to 3 and cannot exceed $\dfrac{P}{3}-1$ (based on a similar breakdown that reached $equation(1)$ in the previous section).

If $c=P-a-b$ is inserted into the Pythagoras Theorem, the formula reduces to:

\begin{align*} a^2+b^2&=(P-a-b)^2\\[.5em] a^2+b^2&=P^2+a^2+b^2-2aP-2bP+2ab\\[.5em] 0&=P^2-2aP-2bP+2ab\\[.5em] 2bP-2ab&=P^2-2aP\\[.5em] b&=\dfrac{P(P-2a)}{2(P-a)} \end{align*}

This formula can be used to calculate $b$ for every $a$ , then the same conditions as for the first solution can be applied to find a valid triplet. Exhaustive search shows that iterating backwards through potential values of $a$ produces the correct output as the first valid triplet, so the loop can be broken early.

1import kotlin.math.hypot
2
3fun pythagoreanTripletA(p: Int): Int {
4    if (p % 2 != 0) return -1
5
6    var product = 0
7    for (a in (p / 3 - 1) downTo 3) {
8        val b = p * (p - 2 * a) / (2 * (p - a))
9        if (a >= b || a * b % 12 != 0) continue
10        val c = p - a - b
11        if (hypot(a.toDouble(), b.toDouble()) == c.toDouble()) {
12            product = a * b * c
13            break
14        }
15    }
16
17    return product
18}

📊 pythagoreanTripletA() gives a significant performance boost by only checking 38 triplets using hypot() when p = 3000.

Alternative Solution

The single loop solution is more than performant enough for our problem goal, but it would be prudent to cover a less straightforward (but still efficient) approach, particularly since it will be a useful algorithm for future problems (as early as Problem 39).

This solution will make use of the so far unmentioned term "Primitive Pythagorean Triplet" that was defined at the beginning of this post.

Based on the fact that all Pythagorean triplets can be derived from a primitive triplet, we can reduce all triplets to their primitive form by dividing out the $\gcd(a,b,c)=d$ factor. Then each primitive triplet can be generated from a unique pair of positive coprime integers $m$ and $n$ based on Euclid's formula:

\begin{align*} \text{given }\;m>n>0,\;\;&\gcd(m,n)=1,\;\;m\nleftrightarrow n\pmod 2\\[.7em] a&=m^2-n^2\\[.5em] b&=2mn\\[.5em] c&=m^2+n^2 \end{align*}

We'll first implement a helper function based on these conditions and equations:

1import kotlin.math.abs
2
3/**
4 * Euclid's formula generates all Pythagorean triplets from 2 numbers, m and n.
5 *
6 * All triplets originate from a primitive one by multiplying them by d = gcd(a,b,c).
7 *
8 * @throws IllegalArgumentException if arguments do not follow m > n > 0, or if not exactly 
9 * one is even, or if they are not co-prime, i.e. gcd(m, n) != 1.
10 */
11fun pythagoreanTriplet(m: Int, n: Int, d: Int): Triple<Int, Int, Int> {
12    require(n in 1 until m) { "Positive integers assumed to be m > n > 0" }
13    require((m % 2 == 0) xor (n % 2 == 0)) { "Integers must be opposite parity" }
14    require(isCoPrime(m, n)) { "Positive integers must be co-prime" }
15    
16    val a = (m * m - n * n) * d
17    val b = 2 * m * n * d
18    val c = (m * m + n * n) * d
19    return Triple(minOf(a, b), maxOf(a, b), c)
20}
21
22fun isCoPrime(x: Int, y: Int): Boolean = gcd(x.toLong(), y.toLong()) == 1L
23
24/**
25 * Recursive calculation of the greatest common divisor of n1 and n2 using the 
26 * Euclidean algorithm.
27 *
28 *      gcd(x, y) = gcd(|y % x|, |x|); where |y| >= |x|
29 *
30 *      gcd(x, 0) = gcd(0, x) = |x|
31 */
32fun gcd(n1: Long, n2: Long): Long {
33    val x = abs(n1)
34    val y = abs(n2)
35    if (x == 0L || y == 0L) return x + y
36    val bigger = maxOf(x, y)
37    val smaller = minOf(x, y)
38    return gcd(bigger % smaller, smaller)
39}

Now that we know how to get a triplet from $(m,n,d)$ , we need to determine how to choose appropriate values for these arguments.

If we insert the three equations from Euclid's formula into the perimeter formula, the latter becomes:

\begin{align*} a+b+c&=P\\[.5em] m^2-n^2+2mn+m^2+n^2&=P\\[.5em] 2m^2+2mn&=P\\[.5em] m^2+mn&=\dfrac{P}{2}\tag{2} \end{align*}

So, at its minimum, $m$ must be greater than 1 (thereby allowing $n=1$ ) and cannot exceed $\Big\lceil\sqrt{\dfrac{P}{2}}\Big\rceil$ . If $equation(2)$ is rearranged to solve for $n$ , it becomes:

\begin{align*} m^2+mn&=\dfrac{P}{2}=L\\[.5em] m(m+n)&=L\tag{3}\\[.5em] n&=\dfrac{L}{m}-m \end{align*}

So, for $n$ to be a positive integer, we must disregard any value of $m$ that is not a factor of $\dfrac{P}{2}$ . This is what our solution looks like so far:

1import kotlin.math.ceil
2import kotlin.math.sqrt
3
4fun pythagoreanTripletAlt(p: Int): Int {
5    if (p % 2 != 0) return -1
6    
7    var product = 0
8    val limit = p / 2
9    val mMax = ceil(sqrt(1.0 * limit)).toInt()
10    for (m in 2 until mMax) {
11        if (limit % m == 0) {
12            // TODO()
13            // Find suitable n values
14            // Iterate through d values & generate triplets
15            // Keep track of max product for each triplet found
16        }
17    }
18    
19    return product
20}

For every valid $m$ , we have to find valid values of $n$ . If we consider $equation(3)$ , we can extract a new value $k=m+n$ as a starting point, which is adjusted based on the parity of $m$ to ensure an opposite parity for $n$ . $k$ will be incremented by 2 (thereby incrementing $n$ while ensuring its parity is retained) without exceeding either $2m$ or $\dfrac{L}{m}$ . The latter can be shown by $equation(3)$ , whereas the former is based on the following:

\begin{align*} \text{given }\;k=m&+n,\;n< m\\[.5em] k-m&< m\\[.5em] k&< 2m \end{align*}

Here is our updated solution:

1import kotlin.math.ceil
2import kotlin.math.sqrt
3
4fun pythagoreanTripletAlt(p: Int): Int {
5    if (p % 2 != 0) return -1
6    
7    var product = 0
8    val limit = p / 2
9    val mMax = ceil(sqrt(1.0 * limit)).toInt()
10    for (m in 2 until mMax) {
11        if (limit % m == 0) {
12            var k = if (m % 2 == 1) m + 2 else m + 1
13            var kMax = limit / m
14            while (kMax % 2 == 0) {
15                kMax /= 2
16            }
17            while (k < 2 * m && k <= kMax) {
18                // TODO()
19                k += 2
20            }
21        }
22    }
23    
24    return product
25}

If the non-primitive versions of the equations are used in the perimeter formula, we can see how a value of $d$ is found:

\begin{align*} a+b+c&=P\\[.5em] (m^2-n^2)d+2mnd+(m^2+n^2)d&=P\\[.5em] 2m^2d+2mnd&=P\\[.5em] d&=\dfrac{P}{2m(m+n)}=\dfrac{L}{mk} \end{align*}

So, to ensure a positive integer $d$ , $k$ needs to be a factor of $\dfrac{L}{m}$ . Also, for $m$ and $n$ to be coprime, $m$ and $k$ must also be coprime.

Finally, once a valid $k$ is found, pythagoreanTriplet() is called with $m$ , $n=k-m$ , and $d=\dfrac{L}{mk}$ , as shown in the final solution:

1import kotlin.math.ceil
2import kotlin.math.sqrt
3
4fun pythagoreanTripletAlt(p: Int): Int {
5    if (p % 2 != 0) return -1
6    
7    var product = 0
8    val limit = p / 2
9    val mMax = ceil(sqrt(1.0 * limit)).toInt()
10    for (m in 2 until mMax) {
11        if (limit % m == 0) {
12            var k = if (m % 2 == 1) m + 2 else m + 1
13            var kMax = limit / m
14            while (kMax % 2 == 0) {
15                kMax /= 2
16            }
17            while (k < 2 * m && k <= kMax) {
18                if (kMax % k == 0 && isCoPrime(k, m)) {
19                    val (a, b, c) = pythagoreanTriplet(m, k - m, limit / (m * k))
20                    product = maxOf(product, a * b * c)
21                }
22                k += 2
23            }
24        }
25    }
26    
27    return product
28}

📊 pythagoreanTripletAlt() generated only 3 Pythagorean triplets when p = 3000. A less direct approach, but well worth having in your repertoire!

Here's a breakdown of the average execution times on my system over 100 iterations at different values of $P$ :

p	Nested Loop	Single Loop (ns)	Alternative (ns)
1000	9.9e+05ns	6.1e+04	2.2e+04
3000	2.03ms	1.1e+05	5.1e+04

📚 Want to read more about Pythagorean triplets? Here are some references about other ways to generate triplets, as well as more in-depth proofs of the theorems used in this post.

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!🏝️📐📐📐