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