Problem Analysis
I am asked to find the perimeter p ≤ 1000 shared by the largest number of distinct integer-sided right triangles. A right triangle with integer sides {a, b, c} obeying a² + b² = c² is called a Pythagorean triple. The naive method - for every candidate perimeter, loop over all side pairs (a, b), set c = p - a - b, and test the squared relation - costs O(p²) per perimeter and O(P²) across all perimeters up to a limit P. That is fine at P = 1000 but becomes hopeless at P = 1,000,000, which the test suite exercises.
The efficient route rests on Euclid's formula. A primitive Pythagorean triple - one whose sides share no common factor, i.e. gcd(a, b, c) = 1 - is generated uniquely by a = m² - n², b = 2mn, c = m² + n² with m > n > 0, gcd(m, n) = 1, and m, n of opposite parity (exactly one even). Here the greatest common divisor being 1 is precisely what makes the triple irreducible. The primitive perimeter simplifies beautifully: p0 = a + b + c = (m² - n²) + 2mn + (m² + n²) = 2m² + 2mn = 2m(m + n). Every non-primitive triple is an integer scaling k of a primitive one, so its perimeter is k · p0. This collapses the problem to enumerating primitive (m, n) pairs and tallying p0 together with all its multiples in a frequency table. Because p0 ≈ 2m², the constraint p0 ≤ P forces m ≲ √(P/2), so the number of primitive pairs is O(P); distributing each pair's multiples adds a harmonic factor, yielding roughly O(P log P) overall - far below the quadratic baseline.
Solution Approach
Euclid parametrisation with multiple distribution (Python and C)
Both p0039_s0.py and p0039_s0.c implement the identical algorithm, so I describe them together and flag the language-specific points as I go.
Step 1 - tight analytic loop bounds. Rather than iterate m and n over a generous range and test the perimeter inside the loop, I invert the perimeter formula algebraically into closed-form limits. For the outer loop, the smallest valid m at a given n is m = n + 1; substituting into p0 = 2m(m + n) ≤ P and solving for n gives the bound n < (int(8√P) - 6) / 8. For the inner loop I fix n and solve the quadratic 2m² + 2mn ≤ P for m via the quadratic formula, giving m < (int(√(4 + 8P)) - 2n) / 4. With these bounds almost every iteration produces an in-range perimeter, so the loop body does real work instead of failing a range test.
Step 2 - parity and coprimality filter. Euclid's formula demands opposite parity and coprimality. I enforce opposite parity structurally: the inner loop starts at m = n + 1 and steps by 2, so m and n always differ in parity. I then keep only pairs with gcd(m, n) == 1, discarding non-primitive generators. The C version uses an iterative inline Euclidean algorithm; the Python version calls math.gcd, which is itself implemented in C.
Step 3 - multiple distribution. For each surviving pair I compute p0 = 2m(m + n), count it once, then count every multiple k · p0 for k = 2, 3, ... up to P. In C I increment a value-indexed array counts[k * p0]; in Python I append each perimeter to a list. The C version also keeps a defensive if (p0 > max_perimeter) break; in case the floating-point bound drifts slightly loose.
Step 4 - find the winner. C performs a single linear scan over counts[], tracking the running maximum and its index. Python feeds the perimeter list to a collections.Counter and reads most_common()[0][0]. Both steps are O(P) and negligible beside the enumeration.