notes.html

problem 66 · 27 lines

<article>
<h3>Problem Analysis</h3>
<p>
    This problem asks me to solve <a class="reference-source" href="https://en.wikipedia.org/wiki/Pell%27s_equation" rel="noopener noreferrer" target="_blank">Pell's equation</a>, x² - Dy² = 1, for every non-square integer D up to 1000, and to report the D whose minimal (fundamental) solution has the largest x. The phrase "minimal solution in x" means the smallest positive-integer pair (x, y) satisfying the equation; every other solution is generated from it. The naive idea - iterate over candidate values of x or y and test - is hopeless, because the fundamental solution can be astronomically large. For D = 661 the value of x exceeds 10^17, so any search bounded by the magnitude of the answer would never terminate. The problem only becomes tractable through a classical theorem of number theory: the fundamental solution of Pell's equation always appears among the <a class="reference-source" href="https://en.wikipedia.org/wiki/Continued_fraction" rel="noopener noreferrer" target="_blank">continued fraction</a> convergents of sqrt(D).
  </p>
<p>
    A <a class="reference-source" href="https://en.wikipedia.org/wiki/Continued_fraction" rel="noopener noreferrer" target="_blank">continued fraction</a> expresses a number as a0 + 1/(a1 + 1/(a2 + ...)), and a convergent is the rational p/q obtained by truncating that expansion at a chosen depth. For a quadratic irrational such as sqrt(D), the expansion is eventually periodic - the integer terms repeat in a fixed-length cycle - which is exactly what lets me compute the whole thing in finite time. The crucial refinement is that one specific convergent yields the fundamental solution, and which one depends on the parity of the period length L. If L is odd, the solution sits at convergent index L - 1; if L is even, at index 2L - 1. The period length is O(sqrt(D)), and the numerators of the convergents grow exponentially with the index, so the digits of x grow as O(sqrt(D)) too - which forces me to use <a class="reference-source" href="https://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic" rel="noopener noreferrer" target="_blank">arbitrary-precision arithmetic</a>. The overall cost is O(D * sqrt(D)) big-integer operations, fast enough to finish in milliseconds for D up to 1000.
  </p>
<h3>Solution Approach</h3>
<h4>Continued-Fraction Convergents for Pell's Equation (Python and C)</h4>
<p>
    Both files implement the same algorithm; Python defers big-integer arithmetic to its runtime, while C carries its own big-integer type. I chose this method because it is correct by theorem and runs in time proportional to the period length rather than the magnitude of the answer.
  </p>
<p>
<strong>Step 1 - Generate the continued fraction.</strong> Starting from a0 = floor(sqrt(D)), I produce each subsequent term with the standard recurrence over three state variables m, n, a: <code>m = n*a - m</code>, <code>n = (D - m*m)/n</code>, <code>a = floor((sqrt(D) + m)/n)</code>. These three integers fully encode the "reduced" form of the quadratic surd at each step. I stop when a term equals 2*a0, which is the classical signal that one full period of the <a class="reference-source" href="https://en.wikipedia.org/wiki/Continued_fraction" rel="noopener noreferrer" target="_blank">continued fraction</a> has been emitted. In C this lives in <code>compute_cf</code>; in Python it is the <code>while</code> loop inside <code>find_fundamental_solution_to_pell_equation</code>.
  </p>
<p>
<strong>Step 2 - Pick the right convergent index.</strong> With the stored array length including a0, the period length is <code>cf_len - 1</code>. If that length is even I evaluate the convergent at index <code>2*cf_len - 3</code>; if odd, at index <code>cf_len - 2</code>. This parity rule comes directly from <a class="reference-source" href="https://en.wikipedia.org/wiki/Pell%27s_equation" rel="noopener noreferrer" target="_blank">Pell theory</a> and means I never have to generate or test more than the one convergent that matters.
  </p>
<p>
<strong>Step 3 - Evaluate the convergent backward.</strong> I evaluate the n-th convergent from the innermost term outward, repeatedly applying <code>term + 1/running_value</code>. To reach indices beyond a single period I wrap the term index with <code>(i - 1) % period_length + 1</code>, holding a0 fixed at position 0 and cycling through the periodic part. The backward evaluation keeps only two big integers (numerator and denominator) alive at any moment, which matters because the convergent index can be large.
  </p>
<p>
<strong>Language-specific differences.</strong> In Python I use <a class="reference-source" href="https://docs.python.org/3/library/fractions.html" rel="noopener noreferrer" target="_blank">fractions.Fraction</a>, which carries an exact <a class="reference-source" href="https://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic" rel="noopener noreferrer" target="_blank">arbitrary-precision</a> numerator and denominator and reduces to lowest terms automatically. In C, where there is no native big integer, I implement a <code>BigInt</code> struct holding up to 512 limbs in base 10^9, little-endian. The only operations I need are addition (<code>bi_add</code>) and multiplication by a small scalar (<code>bi_mul_small</code>); both are simple schoolbook loops. I chose base 10^9 because a limb fits comfortably in a 64-bit word and the product of two limbs fits in 128 bits, which I obtain through GCC's <code>unsigned __int128</code> extension.
  </p>
</article>