Problem 192: Best Approximations

Mathematical Development

Theorem 1. (Continued Fraction Expansion of Quadratic Irrationals.) For every non-square positive integer $d$, $\sqrt{d}$ has an eventually periodic continued fraction expansion

where $a_0 = \lfloor \sqrt{d} \rfloor$ and the period $T$ satisfies $a_T = 2a_0$. The partial quotients are computed via:

Proof. This is a classical result due to Lagrange. The quadratic irrational $\frac{a_0 + m_n}{d_n}$ undergoes a purely periodic transformation under the continued fraction algorithm. Since $d_n \mid (d - m_n^2)$ is maintained as an invariant and $0 < m_n < \sqrt{d}$, there are finitely many possible states, forcing periodicity. $\square$

Theorem 2. (Convergent Recurrence.) The convergents $h_k/k_k$ of a continued fraction $[a_0; a_1, a_2, \ldots]$ satisfy

with $h_{-1} = 1, h_{-2} = 0, k_{-1} = 0, k_{-2} = 1$.

Proof. By induction on $k$, using the matrix identity

$\square$

Theorem 3. (Best Approximation Theorem.) Let $\alpha$ be irrational with convergents $h_k/k_k$. For a given bound $Q$ on the denominator, the best approximation to $\alpha$ with denominator $\leq Q$ is either:

1. A convergent $h_k/k_k$ with $k_k \leq Q$, or
2. A semiconvergent $\frac{m \cdot h_{k} + h_{k-1}}{m \cdot k_{k} + k_{k-1}}$ for some integer $1 \leq m < a_{k+1}$, where $k_{k+1} > Q$.

Proof. This follows from the theory of best rational approximations (see Hardy & Wright, Chapter 10). Every best approximation of the second kind to $\alpha$ is either a convergent or a semiconvergent. The key identity is that convergents alternate in being above and below $\alpha$, with

for all $n$, and semiconvergents interpolate monotonically between consecutive convergents. $\square$

Lemma 1. (Semiconvergent Selection Criterion.) Let $\alpha = \sqrt{d}$, and suppose the next convergent $h_n/k_n$ has $k_n > Q$. The largest valid semiconvergent uses $m = \lfloor (Q - k_{n-2}) / k_{n-1} \rfloor$. This semiconvergent $h_s/k_s$ is better than the previous convergent $h_{n-1}/k_{n-1}$ if and only if

This can be checked with exact integer arithmetic by squaring both sides and using $\alpha^2 = d$:

Proof. Both $|k_s \alpha - h_s|$ and $|k_{n-1} \alpha - h_{n-1}|$ are positive (since $\alpha$ is irrational). Squaring preserves the inequality. Substituting $\alpha^2 = d$ converts the comparison to integer arithmetic. The signs of $k_i \alpha - h_i$ alternate with the parity of the convergent index, but this does not affect the absolute value comparison. $\square$

Editorial

Continued fractions are the right tool here because every best approximation to sqrt(d) with a denominator cap comes from the convergents or the semiconvergents between two consecutive convergents. For each non-square d, we expand the continued fraction only until the next convergent denominator would exceed 10^12. At that moment there are only two serious candidates left: the previous convergent, and the last semiconvergent that still respects the bound.

The implementation keeps the continued-fraction state and the last two convergent denominators. Once the next denominator is too large, it computes the largest admissible semiconvergent multiplier and uses the exact inequality derived from the complete quotient to decide whether that semiconvergent is closer than the previous convergent. That gives the correct denominator for one d, and the global answer is just the sum over all 2 <= d <= 100000 that are not perfect squares.

Pseudocode

Set the denominator bound Q to 10^12 and the running total to 0.

For each d from 2 to 100000:
    Skip d if it is a perfect square.
    Initialize the continued-fraction state for sqrt(d).
    Track the last two convergent denominators.

    Repeatedly generate the next partial quotient and next convergent denominator.
    When the next denominator would exceed Q:
        let t be the largest multiplier that keeps a semiconvergent within Q;
        compare that semiconvergent against the previous convergent
        with the exact integer test coming from the complete quotient;
        add the winning denominator to the total;
        stop processing this d.

Return the total.

Complexity Analysis

• Time: For each non-square $d$, the continued fraction period is $O(\sqrt{d})$ in the worst case. Summing over $d \leq N$: $O(N^{3/2})$ total. For $N = 10^5$, this is approximately $3 \times 10^7$ operations.
• Space: $O(1)$ per value of $d$ (streaming computation), plus $O(\sqrt{N})$ for the perfect-square sieve.

Answer

#include <cassert>
#include <cmath>
#include <iostream>

using namespace std;

const unsigned long long BASE = 1000000000ULL;

struct Big4 {
    unsigned long long limb[4] = {0, 0, 0, 0};
};

Big4 square_u64(unsigned long long value) {
    unsigned long long low = value % BASE;
    unsigned long long high = value / BASE;
    Big4 result;

    unsigned long long term0 = low * low;
    result.limb[0] = term0 % BASE;
    unsigned long long carry = term0 / BASE;

    unsigned long long term1 = 2ULL * low * high + carry;
    result.limb[1] = term1 % BASE;
    carry = term1 / BASE;

    unsigned long long term2 = high * high + carry;
    result.limb[2] = term2 % BASE;
    result.limb[3] = term2 / BASE;
    return result;
}

void multiply_small(Big4& value, unsigned long long factor) {
    unsigned long long carry = 0;
    for (int i = 0; i < 4; ++i) {
        unsigned long long current = value.limb[i] * factor + carry;
        value.limb[i] = current % BASE;
        carry = current / BASE;
    }
}

int compare_big(const Big4& left, const Big4& right) {
    for (int i = 3; i >= 0; --i) {
        if (left.limb[i] < right.limb[i]) {
            return -1;
        }
        if (left.limb[i] > right.limb[i]) {
            return 1;
        }
    }
    return 0;
}

long long best_denominator(int d, long long bound) {
    long long a0 = static_cast<long long>(sqrt(static_cast<long double>(d)));
    while ((a0 + 1) * (a0 + 1) <= d) {
        ++a0;
    }
    while (a0 * a0 > d) {
        --a0;
    }
    if (a0 * a0 == d) {
        return 0;
    }

    long long m = 0;
    long long den = 1;
    long long a = a0;
    long long q_prev2 = 1;
    long long q_prev1 = 0;

    while (true) {
        long long q_curr = a * q_prev1 + q_prev2;
        if (q_curr > bound) {
            long long t = (bound - q_prev2) / q_prev1;
            if (t <= 0) {
                return q_prev1;
            }

            long long rhs_value = den * (2LL * t * q_prev1 + q_prev2) - m * q_prev1;
            unsigned long long rhs = static_cast<unsigned long long>(rhs_value < 0 ? -rhs_value : rhs_value);
            Big4 left = square_u64(rhs);
            Big4 right = square_u64(static_cast<unsigned long long>(q_prev1));
            multiply_small(right, static_cast<unsigned long long>(d));
            if (compare_big(left, right) > 0) {
                return t * q_prev1 + q_prev2;
            }
            return q_prev1;
        }

        q_prev2 = q_prev1;
        q_prev1 = q_curr;
        m = a * den - m;
        den = (d - m * m) / den;
        a = (a0 + m) / den;
    }
}

int main() {
    const long long bound = 1000000000000LL;
    long long total = 0;

    for (int d = 2; d <= 100000; ++d) {
        total += best_denominator(d, bound);
    }

    assert(total == 57060635927998347LL);
    cout << total << '\n';
    return 0;
}

"""
Project Euler Problem 192: Best Approximations

For each non-square d <= 100000, find the denominator of the best
approximation to sqrt(d) with denominator bound 10^12, then sum them.
"""

import math


def best_denominator(d, bound):
    a0 = math.isqrt(d)
    if a0 * a0 == d:
        return 0

    m = 0
    den = 1
    a = a0
    q_prev2, q_prev1 = 1, 0

    while True:
        q_curr = a * q_prev1 + q_prev2
        if q_curr > bound:
            t = (bound - q_prev2) // q_prev1
            if t <= 0:
                return q_prev1

            # Compare the last admissible semiconvergent with the previous convergent
            # using the complete quotient alpha = (sqrt(d) + m) / den.
            rhs = den * (2 * t * q_prev1 + q_prev2) - m * q_prev1
            if rhs * rhs > d * q_prev1 * q_prev1:
                return t * q_prev1 + q_prev2
            return q_prev1

        q_prev2, q_prev1 = q_prev1, q_curr
        m = a * den - m
        den = (d - m * m) // den
        a = (a0 + m) // den


def solve(limit=100000, bound=10**12):
    total = 0
    for d in range(2, limit + 1):
        total += best_denominator(d, bound)
    return total


if __name__ == "__main__":
    answer = solve()
    assert answer == 57060635927998347
    print(answer)