Problem 192: Best Approximations

Project Euler Problem 192 Solution

EulerSolve provides an optimized solution for Project Euler Problem 192, Best Approximations, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every non-square integer $2 \le n \le 100000$, the task is to find the reduced fraction $\frac{p}{q}$ with $1 \le q \le 10^{12}$ that minimizes the error $ \left| \sqrt{n} - \frac{p}{q} \right| $, and then sum all winning denominators $q$. Because $\sqrt{n}$ is irrational whenever $n$ is not a square, there is no tie: exactly one reduced fraction is closest for each such $n$. The difficulty is the denominator limit $10^{12}$. A direct search over all denominators is impossible, so the solution uses the continued fraction of $\sqrt{n}$ and turns the choice into an exact integer comparison between two final candidates. Mathematical Approach The whole argument rests on one classical fact: once a denominator bound is imposed, the best rational approximation to an irrational number must come from the continued-fraction ladder, namely from a convergent or from an intermediate fraction between two consecutive convergents. For $\sqrt{n}$, the code never leaves exact integer arithmetic while it walks that ladder. Continued fractions of $\sqrt{n}$ Let $a_0 = \lfloor \sqrt{n} \rfloor$....

Detailed mathematical approach

Problem Summary

For every non-square integer $2 \le n \le 100000$, the task is to find the reduced fraction $\frac{p}{q}$ with $1 \le q \le 10^{12}$ that minimizes the error $ \left| \sqrt{n} - \frac{p}{q} \right| $, and then sum all winning denominators $q$.

Because $\sqrt{n}$ is irrational whenever $n$ is not a square, there is no tie: exactly one reduced fraction is closest for each such $n$. The difficulty is the denominator limit $10^{12}$. A direct search over all denominators is impossible, so the solution uses the continued fraction of $\sqrt{n}$ and turns the choice into an exact integer comparison between two final candidates.

Mathematical Approach

The whole argument rests on one classical fact: once a denominator bound is imposed, the best rational approximation to an irrational number must come from the continued-fraction ladder, namely from a convergent or from an intermediate fraction between two consecutive convergents. For $\sqrt{n}$, the code never leaves exact integer arithmetic while it walks that ladder.

Continued fractions of $\sqrt{n}$

Let $a_0 = \lfloor \sqrt{n} \rfloor$. For a non-square $n$, the continued fraction of $\sqrt{n}$ is periodic after the first term, and the standard quadratic-irrational state can be written as

$$x_k = \frac{\sqrt{n} + m_k}{d_k}, \qquad a_k = \lfloor x_k \rfloor.$$

The next state is obtained entirely with integers:

$$m_{k+1} = d_k a_k - m_k,$$

$$d_{k+1} = \frac{n - m_{k+1}^2}{d_k},$$

$$a_{k+1} = \left\lfloor \frac{a_0 + m_{k+1}}{d_{k+1}} \right\rfloor.$$

This recurrence is exactly what the implementations maintain. It generates the partial quotients of $\sqrt{n}$ without ever approximating $\sqrt{n}$ numerically.

Convergents and the only fractions that matter

$$\sqrt{n} = [a_0; a_1, a_2, a_3, \dots],$$

then the convergents satisfy

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

with initial values $p_{-2}=0$, $p_{-1}=1$, $q_{-2}=1$, $q_{-1}=0$. These convergents are the principal best approximations.

Between two consecutive convergents there are also semiconvergents, sometimes called intermediate fractions:

$$\frac{p_t}{q_t} = \frac{t p_{k-1} + p_{k-2}}{t q_{k-1} + q_{k-2}}, \qquad 1 \le t < a_k.$$

When a denominator cap $D$ is present, the optimum fraction under $q \le D$ must be either the last convergent whose denominator still fits, or one of these semiconvergents for the next continued-fraction digit. Nothing else can beat them.

Stopping at the first denominator that exceeds the bound

Let $k$ be the first index for which $q_k > D$. Then $q_{k-1} \le D$, so the previous convergent is admissible, while the next full convergent is already too large.

Every admissible semiconvergent on that last edge has denominator

$$q_t = t q_{k-1} + q_{k-2}.$$

The largest allowed choice is therefore

$$t = \left\lfloor \frac{D - q_{k-2}}{q_{k-1}} \right\rfloor,$$

which gives the boundary semiconvergent $q_t \le D < q_{t+1}$. Because $q_k = a_k q_{k-1} + q_{k-2} > D$, this $t$ automatically satisfies $0 \le t < a_k$.

Among semiconvergents on the same edge, larger $t$ means a larger denominator and also a smaller error. In other words, once the bound has cut the edge, only the boundary semiconvergent needs to be compared with the last admissible convergent.

An exact comparison with no floating point

The key identity is

$$\sqrt{n} = \frac{p_{k-2} + x_k p_{k-1}}{q_{k-2} + x_k q_{k-1}}, \qquad x_k = \frac{\sqrt{n} + m_k}{d_k}.$$

From this one obtains the exact errors

$$\left| \sqrt{n} - \frac{p_{k-1}}{q_{k-1}} \right| = \frac{1}{q_{k-1}(x_k q_{k-1} + q_{k-2})},$$

and for the semiconvergent with parameter $t$,

$$\left| \sqrt{n} - \frac{p_t}{q_t} \right| = \frac{x_k - t}{(t q_{k-1} + q_{k-2})(x_k q_{k-1} + q_{k-2})}.$$

So the semiconvergent is better exactly when

$$q_{k-1} x_k < 2 t q_{k-1} + q_{k-2}.$$

Substituting $x_k = \frac{\sqrt{n} + m_k}{d_k}$ turns that into

$$q_{k-1} \sqrt{n} < R,$$

where

$$R = d_k (2 t q_{k-1} + q_{k-2}) - q_{k-1} m_k.$$

This is the integer test used by the implementations. If $R \le 0$, the semiconvergent cannot win. If $R > q_{k-1}(a_0+1)$, then it definitely wins because $\sqrt{n} < a_0 + 1$. Otherwise both sides are positive, so the comparison becomes the exact square test

$$R^2 > n q_{k-1}^2.$$

No floating-point approximation is needed at any stage.

Worked example: $\sqrt{13}$

The continued fraction is

$$\sqrt{13} = [3; 1,1,1,1,6,\dots].$$

The convergents begin

$$\frac{3}{1}, \frac{4}{1}, \frac{7}{2}, \frac{11}{3}, \frac{18}{5}, \frac{119}{33}, \dots$$

Suppose first that $D=20$. The last admissible convergent is $\frac{18}{5}$, because the next denominator $33$ is too large. At this cutoff stage the quadratic-irrational state is $m_k=3$, $d_k=1$, and $q_{k-1}=5$, $q_{k-2}=3$. The boundary parameter is

$$t = \left\lfloor \frac{20-3}{5} \right\rfloor = 3,$$

so the boundary semiconvergent has denominator $q_t = 3 \cdot 5 + 3 = 18$, namely $\frac{65}{18}$. The comparison value is

$$R = 1 \cdot (2 \cdot 3 \cdot 5 + 3) - 5 \cdot 3 = 18.$$

Since $5 \sqrt{13} \approx 18.0278$, we have $R < 5\sqrt{13}$, so $\frac{65}{18}$ is worse than $\frac{18}{5}$. The winning denominator is therefore $5$.

Now change the bound to $D=30$. Then

$$t = \left\lfloor \frac{30-3}{5} \right\rfloor = 5,$$

so the boundary semiconvergent is $\frac{101}{28}$. This time

$$R = 1 \cdot (2 \cdot 5 \cdot 5 + 3) - 5 \cdot 3 = 38,$$

and $38 > 5\sqrt{13}$, so $\frac{101}{28}$ beats $\frac{18}{5}$. The winning denominator becomes $28$. This is exactly the behavior checked by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline: skip perfect squares, generate the continued fraction of each $\sqrt{n}$ with the integer recurrence above, stop at the first convergent whose denominator exceeds $10^{12}$, and then decide between the previous convergent and the single boundary semiconvergent.

Maintaining the continued-fraction state

For each non-square $n$, the implementation starts from $a_0=\lfloor \sqrt{n} \rfloor$ and updates the quadratic-irrational state $(m_k,d_k,a_k)$ together with the convergent numerators and denominators. The recurrence for $(p_k,q_k)$ is synchronized with the recurrence for $(m_k,d_k,a_k)$, so when the denominator first crosses the bound, all data needed for the final comparison is already available.

Choosing the winning denominator

Once the first oversized convergent appears, the implementation sets the provisional answer to the previous denominator $q_{k-1}$. If the denominator bound still leaves room beyond $q_{k-1}$, it computes the boundary value $t = \left\lfloor \frac{D-q_{k-2}}{q_{k-1}} \right\rfloor$ and forms the corresponding semiconvergent denominator $t q_{k-1} + q_{k-2}$.

The final decision is made with the exact criterion above. The C++ and Java implementations use wider integer arithmetic for intermediate values, while Python relies on arbitrary-precision integers automatically. The C++ implementation also performs small exact validations against brute force before launching the full computation, which confirms the continued-fraction logic on tiny cases.

Accumulating the global sum

After the winning denominator for one $n$ is known, it is added to the running total. The C++ implementation distributes different $n$-values across several worker threads; the Python and Java implementations execute the same logic serially. The mathematical decision rule is identical in all three languages.

Complexity Analysis

For one fixed $n$, the loop continues until a convergent denominator exceeds $D$. Since every partial quotient satisfies $a_k \ge 1$, the denominators obey $q_k \ge q_{k-1} + q_{k-2}$, so they grow at least as fast as Fibonacci numbers. Therefore the number of continued-fraction steps is $O(\log D)$.

Each step performs only a constant amount of integer arithmetic, so over all non-squares $n \le N$ the total work is $O(N \log D)$, with $N=100000$ and $D=10^{12}$ in this problem. Memory usage is $O(1)$ per worker aside from the running total. The practical speed comes from replacing a hopeless scan over $10^{12}$ denominators by a very short continued-fraction walk for each $n$.

Footnotes and References

Problem page: https://projecteuler.net/problem=192
Continued fractions: Wikipedia - Continued fraction
Convergents and semiconvergents: Wikipedia - Convergents of continued fractions
Quadratic irrationals and periodic expansions: Wikipedia - Periodic continued fraction
Diophantine approximation: Wikipedia - Diophantine approximation

Problem 192 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128;

constexpr u64 kDenominatorBound = 1'000'000'000'000ULL;
constexpr u64 kMaxN = 100'000ULL;

u64 isqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1) <= x / (r + 1)) {
        ++r;
    }
    while (r > x / r) {
        --r;
    }
    return r;
}

bool is_square(u64 n) {
    const u64 r = isqrt_u64(n);
    return r * r == n;
}

u64 best_denominator_cf(u64 n, u64 d_bound) {
    const u64 a0 = isqrt_u64(n);
    if (a0 * a0 == n) {
        return 0;
    }

    u64 m = 0;
    u64 den = 1;
    u64 a = a0;

    u128 p_prev2 = 0;
    u128 q_prev2 = 1;
    u128 p_prev1 = 1;
    u128 q_prev1 = 0;

    while (true) {
        const u128 p = static_cast<u128>(a) * p_prev1 + p_prev2;
        const u128 q = static_cast<u128>(a) * q_prev1 + q_prev2;

        if (q > d_bound) {
            u128 best_q = q_prev1;

            if (d_bound > q_prev1) {
                const u128 t = (static_cast<u128>(d_bound) - q_prev2) / q_prev1;
                if (t > 0) {
                    const u128 q_semiconv = t * q_prev1 + q_prev2;
                    const i128 rhs = static_cast<i128>(den) *
                                         (2 * static_cast<i128>(t) * static_cast<i128>(q_prev1) +
                                          static_cast<i128>(q_prev2)) -
                                     static_cast<i128>(q_prev1) * static_cast<i128>(m);

                    if (rhs > 0) {
                        const i128 upper =
                            static_cast<i128>(q_prev1) * static_cast<i128>(a0 + 1);

                        if (rhs > upper) {
                            best_q = q_semiconv;
                        } else {
                            const i128 lhs_sq = static_cast<i128>(n) *
                                                static_cast<i128>(q_prev1) *
                                                static_cast<i128>(q_prev1);
                            const i128 rhs_sq = rhs * rhs;
                            if (rhs_sq > lhs_sq) {
                                best_q = q_semiconv;
                            }
                        }
                    }
                }
            }

            return static_cast<u64>(best_q);
        }

        p_prev2 = p_prev1;
        q_prev2 = q_prev1;
        p_prev1 = p;
        q_prev1 = q;

        const u64 m_next = den * a - m;
        const u64 den_next = (n - m_next * m_next) / den;
        const u64 a_next = (a0 + m_next) / den_next;

        m = m_next;
        den = den_next;
        a = a_next;
    }
}

u64 best_denominator_bruteforce(u64 n, u64 d_bound) {
    const long double sq = std::sqrt(static_cast<long double>(n));

    long double best_err = 0.0L;
    u64 best_den = 0;
    bool found = false;

    for (u64 q = 1; q <= d_bound; ++q) {
        const long double v = sq * static_cast<long double>(q);
        const u64 p0 = static_cast<u64>(std::floor(v));

        for (u64 p : {p0, p0 + 1}) {
            if (std::gcd(p, q) != 1) {
                continue;
            }

            const long double err = std::fabsl(static_cast<long double>(p) / static_cast<long double>(q) - sq);
            if (!found || err < best_err) {
                found = true;
                best_err = err;
                best_den = q;
            }
        }
    }

    return best_den;
}

bool run_validations() {
    if (best_denominator_cf(13, 20) != 5) {
        std::cerr << "Validation failed: sqrt(13), d=20 should yield denominator 5.\n";
        return false;
    }

    if (best_denominator_cf(13, 30) != 28) {
        std::cerr << "Validation failed: sqrt(13), d=30 should yield denominator 28.\n";
        return false;
    }

    for (u64 n = 2; n <= 200; ++n) {
        if (is_square(n)) {
            continue;
        }

        for (u64 d = 2; d <= 200; ++d) {
            const u64 fast = best_denominator_cf(n, d);
            const u64 brute = best_denominator_bruteforce(n, d);
            if (fast != brute) {
                std::cerr << "Validation failed: n=" << n << ", d=" << d
                          << ", fast=" << fast << ", brute=" << brute << "\n";
                return false;
            }
        }
    }

    return true;
}

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }

    std::string s;
    while (value > 0) {
        const u64 digit = static_cast<u64>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u128 solve_parallel(u64 max_n, u64 d_bound, unsigned thread_count) {
    if (thread_count == 0) {
        thread_count = 1;
    }

    std::atomic<u64> next_n{2};
    std::vector<u128> partial(thread_count, 0);
    std::vector<std::thread> workers;
    workers.reserve(thread_count);

    for (unsigned tid = 0; tid < thread_count; ++tid) {
        workers.emplace_back([&, tid]() {
            u128 local_sum = 0;

            while (true) {
                const u64 n = next_n.fetch_add(1, std::memory_order_relaxed);
                if (n > max_n) {
                    break;
                }

                if (is_square(n)) {
                    continue;
                }

                local_sum += best_denominator_cf(n, d_bound);
            }

            partial[tid] = local_sum;
        });
    }

    for (auto& worker : workers) {
        worker.join();
    }

    u128 total = 0;
    for (u128 part : partial) {
        total += part;
    }
    return total;
}

}  // namespace

int main() {
    if (!run_validations()) {
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 4;
    }

    const u128 answer = solve_parallel(kMaxN, kDenominatorBound, threads);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

# Problem 192: Best Approximations
# Sum of denominators of best rational approximations to sqrt(n) for non-square n<=100000, with denominator<=10^12.

from math import isqrt

def solve():
    d_bound = 10**12
    total = 0
    for n in range(2, 100001):
        a0 = isqrt(n)
        if a0*a0 == n: continue
        # Continued fraction expansion of sqrt(n)
        m, d, a = 0, 1, a0
        p_prev2, q_prev2 = 0, 1
        p_prev1, q_prev1 = 1, 0
        best_q = 1
        while True:
            p = a * p_prev1 + p_prev2
            q = a * q_prev1 + q_prev2
            if q > d_bound:
                # Check semi-convergent
                best_q = q_prev1
                if d_bound > q_prev1:
                    t = (d_bound - q_prev2) // q_prev1
                    if t > 0:
                        q_sc = t * q_prev1 + q_prev2
                        # Compare errors: |p_sc/q_sc - sqrt(n)| vs |p_prev1/q_prev1 - sqrt(n)|
                        # Use the Stern-Brocot criterion:
                        # q_sc is better if 2*t*q_prev1*d + q_prev2*d > q_prev1*(m + a0 + 1)
                        # where we need exact arithmetic
                        rhs = d * (2*t*q_prev1 + q_prev2) - q_prev1*m
                        if rhs > 0:
                            upper = q_prev1 * (a0 + 1)
                            if rhs > upper:
                                best_q = q_sc
                            else:
                                # Compare rhs^2 vs n * q_prev1^2
                                if rhs*rhs > n * q_prev1 * q_prev1:
                                    best_q = q_sc
                break
            p_prev2, q_prev2 = p_prev1, q_prev1
            p_prev1, q_prev1 = p, q
            m_next = d*a - m
            d_next = (n - m_next*m_next) // d
            a_next = (a0 + m_next) // d_next
            m, d, a = m_next, d_next, a_next
        total += best_q
    print(total)

solve()

Java

import java.math.BigInteger;

public class Euler192 {
    static long isqrt(long x) {
        long r = (long) Math.sqrt((double) x);
        while ((r + 1) <= x / (r + 1))
            r++;
        while (r > x / r)
            r--;
        return r;
    }

    public static void main(String[] args) {
        long dBound = 1000000000000L;
        BigInteger total = BigInteger.ZERO;
        for (long n = 2; n <= 100000; n++) {
            long a0 = isqrt(n);
            if (a0 * a0 == n)
                continue;
            long m = 0, d = 1, a = a0;
            BigInteger pp2 = BigInteger.ZERO, qp2 = BigInteger.ONE, pp1 = BigInteger.ONE, qp1 = BigInteger.ZERO;
            BigInteger bestQ = BigInteger.ONE;
            while (true) {
                BigInteger ba = BigInteger.valueOf(a);
                BigInteger p = ba.multiply(pp1).add(pp2);
                BigInteger q = ba.multiply(qp1).add(qp2);
                if (q.compareTo(BigInteger.valueOf(dBound)) > 0) {
                    bestQ = qp1;
                    BigInteger bdBound = BigInteger.valueOf(dBound);
                    if (bdBound.compareTo(qp1) > 0) {
                        BigInteger t = bdBound.subtract(qp2).divide(qp1);
                        if (t.signum() > 0) {
                            BigInteger qsc = t.multiply(qp1).add(qp2);
                            BigInteger rhs = BigInteger.valueOf(d)
                                    .multiply(BigInteger.TWO.multiply(t).multiply(qp1).add(qp2))
                                    .subtract(qp1.multiply(BigInteger.valueOf(m)));
                            if (rhs.signum() > 0) {
                                BigInteger upper = qp1.multiply(BigInteger.valueOf(a0 + 1));
                                if (rhs.compareTo(upper) > 0)
                                    bestQ = qsc;
                                else if (rhs.multiply(rhs)
                                        .compareTo(BigInteger.valueOf(n).multiply(qp1).multiply(qp1)) > 0)
                                    bestQ = qsc;
                            }
                        }
                    }
                    break;
                }
                pp2 = pp1;
                qp2 = qp1;
                pp1 = p;
                qp1 = q;
                long mn = d * a - m;
                long dn = (n - mn * mn) / d;
                long an = (a0 + mn) / dn;
                m = mn;
                d = dn;
                a = an;
            }
            total = total.add(bestQ);
        }
        System.out.println(total);
    }
}