Problem 192: Best Approximations

Mathematical Foundation

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

For each non-square d <= 100000, find best rational approximation p/q to sqrt(d) with q <= 10^12. Sum all q. We generate continued fraction and convergents. We then actually use standard initialization. Finally, loop.

Pseudocode

Generate continued fraction and convergents
Actually use standard initialization:
loop
Check semiconvergent
Compare |ks*sqrt(d) - hs| vs |k1*sqrt(d) - h1|
for d from 2 to 100000

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 <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef __int128 lll;
typedef unsigned long long ull;

// Integer square root
ll isqrt(ll n) {
    ll x = (ll)sqrt((double)n);
    while (x * x > n) x--;
    while ((x+1)*(x+1) <= n) x++;
    return x;
}

ll best_denominator(ll d, ll Q) {
    ll a0 = isqrt(d);
    if (a0 * a0 == d) return 0; // perfect square, skip

    // Continued fraction expansion of sqrt(d)
    // We track convergents h2/k2 (two steps back), h1/k1 (one step back)
    ll m = 0, dd = 1, a = a0;
    ll h2 = 0, h1 = 1, k2 = 1, k1 = 0;

    // Process a0
    ll h = a0 * h1 + h2;
    ll k = a0 * k1 + k2;
    h2 = h1; h1 = h;
    k2 = k1; k1 = k;

    while (true) {
        m = a * dd - m;
        dd = (d - m * m) / dd;
        a = (a0 + m) / dd;

        ll knext = a * k1 + k2;
        if (knext > Q) {
            // The full convergent exceeds Q.
            // Check semiconvergent with max valid m_sc
            ll m_sc = (Q - k2) / k1;
            if (m_sc < 1) return k1; // no valid semiconvergent

            ll qs = m_sc * k1 + k2;
            ll ps = m_sc * h1 + h2;

            // Compare |qs*sqrt(d) - ps| vs |k1*sqrt(d) - h1|
            // |qs*sqrt(d) - ps|^2 = qs^2*d - 2*ps*qs*sqrt(d) + ps^2
            // But we can use: for convergent p_{n-1}/q_{n-1}, the error is
            // 1/(q_{n-1}*(a_n*q_{n-1}+q_{n-2})) approximately.
            // More precisely:
            // |q*sqrt(d)-p| compared by squaring is tricky due to sqrt.
            //
            // Instead use: semiconvergent is better iff
            //   |qs*sqrt(d) - ps| < |k1*sqrt(d) - h1|
            // Both sides positive (after abs). Square both sides:
            //   (qs^2*d + ps^2 - 2*ps*qs*sqrt(d)) < (k1^2*d + h1^2 - 2*h1*k1*sqrt(d))
            // Let A = qs^2*d + ps^2 - k1^2*d - h1^2
            // Let B = 2*(ps*qs - h1*k1)
            // We need A < B*sqrt(d)
            // If B > 0: A < B*sqrt(d) iff A < 0 or A^2 < B^2*d
            // If B <= 0: A < B*sqrt(d) iff A < 0 and A^2 > B^2*d

            lll A = (lll)qs*qs*d + (lll)ps*ps - (lll)k1*k1*d - (lll)h1*h1;
            lll B = 2*((lll)ps*qs - (lll)h1*k1);

            bool semi_better;
            if (B > 0) {
                if (A <= 0) semi_better = true;
                else semi_better = (A * A < B * B * d);
            } else {
                if (A >= 0) semi_better = false;
                else semi_better = (A * A > B * B * d);
            }

            return semi_better ? qs : k1;
        }

        ll hnext = a * h1 + h2;
        h2 = h1; h1 = hnext;
        k2 = k1; k1 = knext;

        if (a == 2 * a0 && k1 <= Q) {
            // Completed one period; if still within bound, continue
            // (the CF is periodic, it will repeat)
        }
    }
}

int main() {
    const ll Q = 1000000000000LL; // 10^12
    const int N = 100000;
    ll total = 0;
    for (int d = 2; d <= N; d++) {
        ll sq = isqrt(d);
        if (sq * sq == d) continue;
        total += best_denominator(d, Q);
    }
    cout << total << endl;
    return 0;
}

"""
Problem 192: Best Approximations

For each non-square d <= 100000, find best rational approximation p/q
to sqrt(d) with q <= 10^12. Sum all q.
"""

import math

def best_denom(d, Q):
    a0 = int(math.isqrt(d))
    if a0 * a0 == d:
        return 0

    # Continued fraction of sqrt(d)
    m, dd, a = 0, 1, a0

    # Convergents: h2/k2 = p_{n-2}/q_{n-2}, h1/k1 = p_{n-1}/q_{n-1}
    h2, h1 = 0, a0
    k2, k1 = 1, 1

    while True:
        m = a * dd - m
        dd = (d - m * m) // dd
        a = (a0 + m) // dd

        knext = a * k1 + k2
        if knext > Q:
            # Check semiconvergent
            m_sc = (Q - k2) // k1
            if m_sc < 1:
                return k1

            qs = m_sc * k1 + k2
            ps = m_sc * h1 + h2

            # Compare |qs*sqrt(d) - ps| vs |k1*sqrt(d) - h1|
            # Using exact integer arithmetic:
            # Need: (qs^2*d + ps^2 - k1^2*d - h1^2) < 2*(ps*qs - h1*k1)*sqrt(d)
            A = qs * qs * d + ps * ps - k1 * k1 * d - h1 * h1
            B = 2 * (ps * qs - h1 * k1)

            if B > 0:
                if A <= 0:
                    semi_better = True
                else:
                    semi_better = (A * A < B * B * d)
            else:
                if A >= 0:
                    semi_better = False
                else:
                    semi_better = (A * A > B * B * d)

            return qs if semi_better else k1

        hnext = a * h1 + h2
        h2, h1 = h1, hnext
        k2, k1 = k1, knext

def solve():
    Q = 10**12
    N = 100000
    total = 0
    for d in range(2, N + 1):
        sq = int(math.isqrt(d))
        if sq * sq == d:
            continue
        total += best_denom(d, Q)
    return total

answer = solve()
assert answer == 43289071131872028
print(answer)