The 1000-perimeter Pythagorean triplet: one loop, then none

There is exactly one right triangle with integer sides whose perimeter is 1000. Project Euler #9 asks for the product of those sides — and the honest fun is in refusing to search for it. This is a walk from a nested-loop scan, down to a single loop the algebra hands you for free, down to a divisor enumeration straight out of Euclid.

The problem, in my own words

A Pythagorean triplet is three whole numbers $a < b < c$ that satisfy $a^2 + b^2 = c^2$ — the side lengths of an integer right triangle. The classic example is $3, 4, 5$: $9 + 16 = 25$. Project Euler problem #9 fixes the perimeter and asks for the sides: there is exactly one such triplet with $a + b + c = 1000$, and the task is to compute $abc$.

Per Project Euler's request not to spoil the puzzle, I won't print the numeric product here. Every script below is complete and runs in microseconds — clone it and the answer is one command away.

Why Python

The whole computation is a handful of integer operations, so raw throughput is irrelevant — the bottleneck is thinking, not the CPU. What I do want is for the algebra I derive on paper to map one-to-one onto the code with no ceremony in between: exact integers, a % that means real divisibility rather than a float that's "close to" an integer, and math.gcd sitting in the standard library for the number-theoretic version. Python gives all three. There is no overflow bookkeeping to distract from the reasoning (even the three-nested-loop version never exceeds $10^6$), and the closed-form line reads almost exactly like the equation it came from. When a problem is small enough that clarity is the only thing you can optimise, optimise clarity.

The obvious approach, and why to leave it behind

The brute-force reading of the problem is a triple loop: try every $a$, every $b$, every $c$, keep the combination that is both a triplet and sums to 1000. That is $O(p^3)$ and immediately wasteful, because the perimeter constraint is not a filter to check at the end — it's a degree of freedom you've already spent. Once you pick $a$ and $b$, the third side is forced:

So the third loop collapses to an assignment, and you're left with $O(p^2)$: pick $a$, pick $b$, let $c$ fall out, test whether $a^2 + b^2 = c^2$. Two more observations trim the ranges. Since $a < b < c$ and they sum to $p$, the smallest side must be below $p/3$ (three copies of the smallest can't reach the total), so $a < p/3$. And $b < c$ means $b < (p-a)/2$. That's brute_quadratic below — correct, and a fine sanity anchor, but it still touches on the order of $p^2/6 \approx 1.7\times10^5$ candidate pairs.

Collapsing the second loop with algebra

We can do better by spending the second constraint the way we spent the first. We have two equations and, once $a$ is chosen, one unknown too many — so solve for $b$ directly. Substitute $c = p - a - b$ into $a^2 + b^2 = c^2$:

The $a^2$ and $b^2$ cancel on both sides, which is the whole point — the quadratic in $b$ is secretly linear:

Now there is no inner loop at all. Scan $a$ from 1 upward, compute the single candidate $b$, accept it only when the division is exact (an integer side) and the ordering $a < b < c$ holds. That's $O(p)$, and because the acceptance test is "does den divide num," it's pure integer arithmetic with no floating-point "is this close enough to a whole number" fudge.

The curve $b(a)$ is worth picturing. For $p = 1000$ it is a smooth, gently falling arc, and the search is just a hunt along it for the lone lattice point that is simultaneously an integer and correctly ordered:

The number-theoretic shortcut: Euclid's parametrisation

The linear scan is already fast, but it still walks past ~333 values of $a$ to find one hit. We can jump nearly straight to the answer by parametrising all Pythagorean triples instead of testing candidates.

Euclid's formula says every primitive triplet — one where $a, b, c$ share no common factor — is generated by a pair of integers $m > n > 0$ that are coprime and of opposite parity (one odd, one even):

Every triplet whatsoever is a primitive one scaled by some positive integer $k$. So the general perimeter is

Setting that equal to $p = 1000$ turns the whole problem into a factoring question:

Now enumerate structurally. The factor $m(m+n)$ is strictly between $m^2$ and $2m^2$, so $m$ only ranges up to $\sqrt{p/2} \approx 22$ — that's $O(\sqrt p)$. For each $m$ dividing 500, the co-factor $k(m+n)$ is fixed; sweep the divisor $s = m + n$ (which must lie in $(m, 2m)$ so that $0 < n < m$), recover $n = s - m$ and $k$, and keep the pair that is coprime and opposite-parity. That's a handful of iterations rather than hundreds.

The code

All three live in one file, pe9.py, with a shared benchmark and a tracer that prints exactly which $(m, n, k)$ the Euclid search visits.

#!/usr/bin/env python3
"""Project Euler #9 — the a+b+c=1000 Pythagorean triplet, three ways."""
import math
import time

P = 1000

def brute_quadratic(perim):
    """O(p^2): pick a, pick b, let c fall out of the perimeter, test."""
    for a in range(1, perim // 3):          # a < b < c forces a < p/3
        for b in range(a + 1, (perim - a) // 2):
            c = perim - a - b
            if a * a + b * b == c * c:
                return a, b, c
    return None

def algebraic_linear(perim):
    """O(p): one loop over a, solve for b in closed form.

        c = p - a - b, and a^2 + b^2 = c^2  =>
        b = p*(p - 2a) / (2*(p - a))
    Only accept integer b with a < b < c.
    """
    for a in range(1, perim // 3):
        num = perim * (perim - 2 * a)
        den = 2 * (perim - a)
        if num % den:
            continue
        b = num // den
        c = perim - a - b
        if a < b < c:
            return a, b, c
    return None

def euclid(perim):
    """O(sqrt(p)) via Euclid's parametrisation.

    Every triple is k*(m^2-n^2, 2mn, m^2+n^2), m>n>0 coprime, opposite
    parity. Its perimeter is 2*k*m*(m+n). Solve 2*k*m*(m+n) = p by
    enumerating m up to sqrt(p/2), then the factor s=m+n, then k.
    """
    if perim % 2:
        return None
    half = perim // 2                         # = k*m*(m+n)
    m = 2
    while m * m <= half:                       # m(m+n) < 2m^2 <= half
        if half % m == 0:
            rest = half // m                   # = k*(m+n)
            s = m + 1                           # candidate m+n, s>m
            while s < 2 * m:                    # n = s-m in (0, m)
                if rest % s == 0:
                    n = s - m
                    if math.gcd(m, n) == 1 and (m - n) % 2 == 1:
                        k = rest // s
                        a = k * (m * m - n * n)
                        b = k * (2 * m * n)
                        c = k * (m * m + n * n)
                        return tuple(sorted((a, b, c)))
                s += 1
        m += 1
    return None

def trace_euclid(perim):
    """Print every (m,n,k) the Euclid search actually visits and accepts."""
    half = perim // 2
    print(f"target k*m*(m+n) = p/2 = {half}")
    m = 2
    while m * m <= half:
        if half % m == 0:
            rest = half // m
            s = m + 1
            while s < 2 * m:
                if rest % s == 0:
                    n = s - m
                    g = math.gcd(m, n)
                    parity = (m - n) % 2 == 1
                    k = rest // s
                    ok = g == 1 and parity
                    print(f"  m={m:>2} n={n:>2} k={k:>2}  gcd={g} "
                          f"opp_parity={parity}  {'ACCEPT' if ok else 'skip'}")
                s += 1
        m += 1

def bench(fn, perim, reps=1000):
    t0 = time.perf_counter()
    for _ in range(reps):
        r = fn(perim)
    dt = (time.perf_counter() - t0) / reps
    return r, dt

if __name__ == "__main__":
    print("== Euclid search trace, p =", P, "==")
    trace_euclid(P)
    print()
    for name, fn, reps in [("brute_quadratic", brute_quadratic, 100),
                           ("algebraic_linear", algebraic_linear, 10000),
                           ("euclid", euclid, 10000)]:
        (a, b, c), dt = bench(fn, P, reps)
        assert a * a + b * b == c * c and a + b + c == P
        print(f"{name:>17}: a<b<c found, a+b+c={a+b+c}, "
              f"abc has {len(str(a*b*c))} digits, "
              f"mean {dt*1e6:8.2f} us over {reps} runs")

A few implementation notes that earn their space:

• if num % den: in algebraic_linear is the integer test doing double duty — a non-zero remainder means $b$ would be fractional, i.e. no integer triangle for this $a$, so we skip without ever forming a float.
• In euclid, the loop guard while m * m <= half is the $O(\sqrt p)$ bound made literal: since $m(m+n) < 2m^2$ and $m(m+n)$ divides half, we know $m^2 < \text{half}$, so $m$ never exceeds $\sqrt{500} \approx 22.4$.
• The inner while s < 2 * m keeps $n = s - m$ strictly between 0 and $m$, which is exactly the $m > n > 0$ Euclid requires — we never generate an out-of-range or duplicate parametrisation.
• The assert in the driver is not decoration: it re-checks $a^2+b^2=c^2$ and $a+b+c=p$ on whatever each function returns, so a wrong triple would crash the run rather than print a plausible lie.

Running it

$ python3 pe9.py
== Euclid search trace, p = 1000 ==
target k*m*(m+n) = p/2 = 500
  m= 4 n= 1 k=25  gcd=1 opp_parity=True  ACCEPT
  m=20 n= 5 k= 1  gcd=5 opp_parity=True  skip

  brute_quadratic: a<b<c found, a+b+c=1000, abc has 8 digits, mean 12276.80 us over 100 runs
 algebraic_linear: a<b<c found, a+b+c=1000, abc has 8 digits, mean    32.22 us over 10000 runs
           euclid: a<b<c found, a+b+c=1000, abc has 8 digits, mean     1.56 us over 10000 runs

All three agree on the same triplet (the driver's assert guarantees it satisfies both equations), and the product has 8 digits — which is all I'll say about the number. The performance story is the interesting part: the quadratic scan averages 12.3 ms, the algebraic single loop 32 µs — a ~380× cut for one page of algebra — and Euclid's parametrisation 1.6 µs, another ~20× on top for thinking in terms of triples instead of candidates. Every one of them clears Project Euler's one-minute rule by roughly six orders of magnitude, but the gap between them is the whole lesson.

A worked trace

The tracer prints precisely what the Euclid search does, and it's short enough to follow by hand. We need $k\,m\,(m+n) = 500$.

Two candidates survive the divisibility sieve, and only one survives the number theory. Take the accepted row: $m = 4$, $n = 1$, $k = 25$. Euclid's formula gives the primitive triple

which is the well-known $8, 15, 17$ right triangle (perimeter 40). Scale by $k = 25$ and the perimeter scales to $25 \times 40 = 1000$ — exactly our target. The rejected row, $m = 20, n = 5$, would have produced the same triangle a second time in un-reduced form ($\gcd = 5$ soaks up the very $k = 25$ we already applied); the coprimality test is what stops us from double-counting it. So the search finds the answer, proves it's the only answer for $p = 1000$, and explains structurally why: there is a single scaled copy of the $8\!-\!15\!-\!17$ triangle whose perimeter hits 1000.

What this costs

The one-minute rule is not in danger here — even the deliberately naive quadratic finishes in about twelve milliseconds. The point of walking down the three versions isn't speed for its own sake; it's that each rewrite comes from spending a constraint you already have instead of checking it after the fact. The perimeter kills the third loop. The Pythagorean equation, run through the cancellation of $a^2 + b^2$, kills the second. And recognising that "integer right triangles" is a parametrised family rather than a haystack turns the last loop into a short divisor walk over $p/2$. That progression — filter, then solve, then parametrise — is the same ladder that scales a lot of Project Euler problems from "runs in a minute" to "runs before your finger leaves Enter."

Reproduce it

Save pe9.py from the block above and run python3 pe9.py. The tracer prints the accepted $(m, n, k)$, each solver prints the digit count of $abc$ and its mean runtime, and the built-in assert refuses to let a wrong triple through. The eight-digit product is yours to read off in one more line — Project Euler would rather you earn it.

References

• Project Euler, Problem #9 — Special Pythagorean Triplet: https://projecteuler.net/problem=9
• Euclid's formula for Pythagorean triples (primitive generation via coprime, opposite-parity $m > n$): a standard result; see e.g. the Wikipedia entry on Pythagorean triple. Timings measured in a Kali sandbox with CPython 3, time.perf_counter, averaged over the run counts shown.