Project Euler Solution #9: Special Pythagorean triplet

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,

For example, 3² + 4² = 9 + 16 = 25 = 5².

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

Solution

A problem like this is perfect for Haskell, since we can transform the mathematical description par for par into a programm:

solution n = head [a * b * c |
    c <- [1..n],
    b <- [1..c],
    a <- [1..b],
    a^2 + b^2 == c^2,
    a + b + c == n]

However, it takes quite some time to find the solution with this brute force attempt. One improvement is to use the fact that \(a+b+c=n\) which let us construct the variable c as \(c=n - a - b\) (we could construct c as \(c=\sqrt{a^2+b^2}\), but this is computational much more expensive and of course we could construct a and b but this would require another check if the result becomes negative). This minimal change reduces the execution time tremendously, since we reduced the algorithm from \(O(n^3)\) to \(O(n^2)\):

solution n = head [a * b * c |
  b <- [1..n],
  a <- [1..b],
  let c = n - a - b,
  a^2 + b^2 == c^2]

In order to mind the given relation \(a<b<c\), we would need to add a check of \(b<c\), but as \(a^2+b^2=c^2 \Rightarrow b<c\) for \(a>0\), we can skip this test. But based on the fact that \(a<b<c\) and \(a+b+c=n\), we can conclude even more: The maximum possible \(a\) and maximum possible \(b\) is:

\[a<\frac{n}{3} - 1\]

\[b<\frac{n}{2} - 1\]

Since \(a\) as this smallest number plus the sum of two other numbers shall not exceed \(n\). The same takes effect for \(b\) as this smallest number plus \(c\) shall not exceed \(n\). With this information, we can further reduce the search space:

solution n = head [a * b * c |
  a <- [1..quot n 3],
  b <- [a..quot n 2],
  let c = n - a - b,
  a^2 + b^2 == c^2]

We learned quite a lot already, but the search space is still too large, so back to the drawing board. We know \(a+b=n-c\) and \(a^2+b^2=c^2\), which lets us formulate

\[\begin{array}{rrl} & a^2 + b^2 =& c^2\\ \Leftrightarrow & a^2 + 2ab + b^2 =& c^2 + 2ab\\ \Leftrightarrow & (a+b)^2 =& c^2 + 2ab\\ \Leftrightarrow & 2ab =& (\underbrace{n-c}_{a+b})^2 - c^2 \end{array}\]

When subtracting the resulting equations from the starting equation, we get

\[\begin{array}{rrl} & a^2 - 2ab + b^2 =& c^2 - (n-c)^2 + c^2\\ \Leftrightarrow & (a-b)^2 =& c^2 - n^2 + 2nc - c^2 + c^2\\ \Leftrightarrow & (a-b)^2 =& c^2 - n^2 + 2nc \end{array}\]

Since \((a-b)^2 = c^2 - n^2 + 2nc\), we know the right part must be square. We achieved, that we completely removed a and b! As we constructed the triplet this way, we find a solution if the integer square root exists. Implementing all this knowledge leads to a solution which only depends on \(c\) and calculates in linear time \(O(n)\):

solution n = head [ a * b * c |
  c <- [1 + quot n 3 .. quot n 2],
  let sqa_b = c * c - n * n + 2 * n * c,
  let a_b = floor(sqrt(fromIntegral sqa_b)),
  let b = quot (n - c + a_b) 2,
  let a = n - b - c,
  a_b * a_b == sqa_b]

Or the same implemented in JavaScript:

function solution(n) {

  for (var c = Math.floor(n / 3 + 1); c < n / 2; c++) {

    var sqa_b = c * c - n * n + 2 * n * c
    var a_b = Math.floor(Math.sqrt(sqa_b));

    if (a_b * a_b == sqa_b) {
      var b = (n - c + a_b) / 2;
      var a = n - b - c;
      return a * b * c;
    }
  }
  return -1
}
solution(1000);