Continued Fractions & Pell's Equation

Intuition, Vocabulary and Foundations

A continued fraction records a real number as a sequence of integer quotients, much as the Euclidean algorithm records successive divisions. For example, $\sqrt{2}=[1;\overline{2}]$ produces the rational approximations $1/1$, $3/2$, $7/5$, and $17/12$. Each new convergent is obtained from the previous two, so the method behaves like a small exact state machine rather than a floating-point search.

For coefficients $a_0,a_1,\ldots$, define $p_{-2}=0$, $p_{-1}=1$, $q_{-2}=1$, and $q_{-1}=0$. The recurrence

$$p_k=a_kp_{k-1}+p_{k-2},\qquad q_k=a_kq_{k-1}+q_{k-2}$$

produces the convergent $p_k/q_k$. The determinant identity $p_kq_{k-1}-p_{k-1}q_k=(-1)^{k-1}$ proves that adjacent convergents are reduced and extremely close. For nonsquare $D>0$, the continued fraction of $\sqrt D$ is periodic. Its convergents eventually supply the least positive solution of Pell's equation $x^2-Dy^2=1$; for $D=2$, the first such solution is $(x,y)=(3,2)$.

Step-by-Step Derivations and Algorithms

Write a quadratic-surd state as $(m_k,d_k,a_k)$, initially $(0,1,\lfloor\sqrt D\rfloor)$. Exact rearrangement of the reciprocal remainder gives

$$m_{k+1}=d_ka_k-m_k,$$ $$d_{k+1}=\frac{D-m_{k+1}^2}{d_k},$$ $$a_{k+1}=\left\lfloor\frac{a_0+m_{k+1}}{d_{k+1}}\right\rfloor.$$

All divisions are exact consequences of the preceding state. For a nonsquare $D$, the state returns to the period boundary when $a_k=2a_0$. If the period length is $L$, the fundamental solution occurs at convergent index $L-1$ when $L$ is even and at index $2L-1$ when $L$ is odd. This parity rule follows from the alternating determinant identity and the norm of the convergent reached after one period.

To generate convergents, process coefficients once while retaining only the previous two numerator-denominator pairs. Test the exact certificate $p_k^2-Dq_k^2=1$ rather than relying on a decimal approximation. Once the fundamental solution $(x_1,y_1)$ is known, later solutions satisfy

$$(x_{t+1},y_{t+1})=(x_1x_t+Dy_1y_t,\;x_1y_t+y_1x_t),$$

which is the coefficient expansion of $(x_1+y_1\sqrt D)^{t+1}$.

flowchart TD
    A[Validate positive nonsquare D] --> B[Expand exact surd states]
    B --> C[Detect period at coefficient 2a0]
    C --> D{Period length even?}
    D -->|Yes| E[Use one period of coefficients]
    D -->|No| F[Use two periods of coefficients]
    E --> G[Build convergent and verify Pell identity]
    F --> G

Implementation, Complexity and Engineering

The reference implementation uses isqrt, so square detection and every transition are exact. It rejects $D\le0$ and square $D$ because the positive nonsquare Pell workflow does not apply. Python integers grow automatically; fixed-width implementations must either use checked big integers or reject results that exceed their documented range.

The period is stored for clarity and reuse. A streaming implementation can reduce auxiliary memory when only one solution is needed, but storing the period makes parity handling, traces, and tests easier to audit. The strongest oracle is always the defining identity, supplemented by known minimal solutions and comparisons with brute force on small values.

from math import isqrt


def sqrt_continued_fraction(d_value):
    if not isinstance(d_value, int) or isinstance(d_value, bool):
        raise TypeError("D must be an integer")
    if d_value <= 0:
        raise ValueError("D must be positive")
    a0 = isqrt(d_value)
    if a0 * a0 == d_value:
        raise ValueError("D must be nonsquare")

    period = []
    m_value, denominator, coefficient = 0, 1, a0
    while True:
        m_value = denominator * coefficient - m_value
        denominator = (d_value - m_value * m_value) // denominator
        coefficient = (a0 + m_value) // denominator
        period.append(coefficient)
        if coefficient == 2 * a0:
            return a0, period


def convergent(coefficients):
    p_prev2, p_prev1 = 0, 1
    q_prev2, q_prev1 = 1, 0
    for coefficient in coefficients:
        p_value = coefficient * p_prev1 + p_prev2
        q_value = coefficient * q_prev1 + q_prev2
        p_prev2, p_prev1 = p_prev1, p_value
        q_prev2, q_prev1 = q_prev1, q_value
    return p_prev1, q_prev1


def pell_fundamental(d_value):
    a0, period = sqrt_continued_fraction(d_value)
    length = len(period)
    coefficient_count = length if length % 2 == 0 else 2 * length
    coefficients = [a0]
    for index in range(coefficient_count - 1):
        coefficients.append(period[index % length])
    x_value, y_value = convergent(coefficients)
    if x_value * x_value - d_value * y_value * y_value != 1:
        raise ArithmeticError("Pell certificate failed")
    return x_value, y_value


def pell_solutions(d_value, count):
    if not isinstance(count, int) or isinstance(count, bool):
        raise TypeError("count must be an integer")
    if count < 0:
        raise ValueError("count must be nonnegative")
    if count == 0:
        return []
    x1, y1 = pell_fundamental(d_value)
    results = []
    x_value, y_value = 1, 0
    for _ in range(count):
        x_value, y_value = (
            x1 * x_value + d_value * y1 * y_value,
            x1 * y_value + y1 * x_value,
        )
        results.append((x_value, y_value))
    return results


known = {
    2: (3, 2),
    3: (2, 1),
    5: (9, 4),
    6: (5, 2),
    7: (8, 3),
    13: (649, 180),
    61: (1766319049, 226153980),
}
for d_value, expected in known.items():
    actual = pell_fundamental(d_value)
    assert actual == expected
    assert actual[0] ** 2 - d_value * actual[1] ** 2 == 1

assert sqrt_continued_fraction(2) == (1, [2])
assert pell_solutions(2, 4) == [(3, 2), (17, 12), (99, 70), (577, 408)]
assert pell_solutions(2, 0) == []

for invalid in (0, -2, 1, 4, 9):
    try:
        pell_fundamental(invalid)
    except ValueError:
        pass
    else:
        raise AssertionError("invalid D was accepted")

print("all tests passed")

Applications, Case Studies and Exercises

Continued fractions solve two related tasks: they produce provably strong rational approximations, and their periodic quadratic-surd states expose exact Pell solutions. The approachable case is approximation: successive convergents to $\sqrt2$ improve rapidly while their determinant certificate proves reducedness. The professionally realistic case is an exact solver that processes many $D$ values, rejects squares, records period lengths, returns minimal solutions, and emits a certificate that a separate system can verify with one multiplication-subtraction check.

For analysis, distinguish recurrence steps from bit complexity: the number of states may be modest even when the fundamental solution contains many digits. For testing, enumerate small $D$, compare minimality against bounded brute force where feasible, and always verify $x^2-Dy^2=1$. Then add known stress cases such as $D=61$ to expose accidental fixed-width arithmetic.