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

Let $x$ be a real number and $d$ a positive integer. The statement calls a fraction $r/s$ in lowest terms with $s\le d$ the best approximation to $x$ for the denominator bound $d$ if every rational number closer to $x$ than $r/s$ has a denominator larger than $d$:

$$\left|\frac pq-x\right|\lt\left|\frac rs-x\right|\ \Longrightarrow\ q\gt d.$$

In other words, $r/s$ is the fraction nearest to $x$ among all fractions whose denominator is at most $d$. For the irrational numbers of this problem that fraction is unique: two different fractions at the same distance from $x$ would lie on opposite sides of $x$, so $x$ would be their midpoint, which is rational. The statement gives two checkpoints: the best approximation to $\sqrt{13}$ is $18/5$ for the bound $20$ and $101/28$ for the bound $30$. We must add up the denominators of the best approximations to $\sqrt n$ for the bound $10^{12}$, over all integers $n$ with $1\lt n\le100000$ that are not perfect squares. There are $99684$ such $n$, since $315$ of the numbers $2,3,\dots,100000$ are squares. ^[1]

The whole solution in one sentence: we expand $\sqrt n$ into its continued fraction, using only small whole numbers, until the denominators of its convergents pass the bound; the best approximation is then one of just two fractions built from the last two convergents, and comparing two integers decides which one.

Mathematical Approach

We watch the best approximations to $\sqrt{13}$ as the bound grows, recall continued fractions, prove that only two fractions can win, and turn the choice between them into a comparison of integers.

1. Small cases: the best approximations to $\sqrt{13}$

We have $\sqrt{13}=3.605551\ldots$. Trying every denominator $q\le D$ with the two nearest numerators $\lfloor q\sqrt{13}\rfloor$ and $\lfloor q\sqrt{13}\rfloor+1$ finds the best approximation for the bound $D$; for $D=1,2,\dots,80$ it gives:

$4/1$ for $D=1$, $7/2$ for $D=2$ and $11/3$ for $D=3,4$;
$18/5$ for $5\le D\le22$;
$83/23$ for $23\le D\le27$ and $101/28$ for $28\le D\le32$;
$119/33$ for $33\le D\le37$, $137/38$ for $38\le D\le70$ and $256/71$ for $71\le D\le80$.

Both checkpoints appear in this list. The winners $4/1$, $7/2$, $11/3$, $18/5$, $119/33$, $137/38$ and $256/71$ will turn out to be the convergents of $\sqrt{13}$. The other two, $83/23$ and $101/28$, belong to the family $(18t+11)/(5t+3)$, which runs from $11/3$ at $t=0$ to $119/33$ at $t=6$. Its members $29/8$, $47/13$ and $65/18$, for $t=1,2,3$, never win. The closest call is $D=20$: the fraction $65/18$ is allowed, but its error $0.0055598\ldots$ is slightly larger than the error $0.0055513\ldots$ of $18/5$.

2. Continued fractions and convergents

Every irrational number $x$ has a unique expansion

$$x=[a_0;a_1,a_2,\dots]=a_0+\cfrac{1}{a_1+\cfrac{1}{a_2+\cdots}},$$

produced by $x_0=x$, $a_k=\lfloor x_k\rfloor$ and $x_{k+1}=1/(x_k-a_k)$. The integers $a_k$ are the partial quotients, with $a_k\ge1$ for $k\ge1$, and the numbers $x_k$ are the complete quotients, with $x_k\gt1$ for $k\ge1$. Cutting the expansion after $a_k$ gives the convergent $p_k/q_k=[a_0;a_1,\dots,a_k]$. For example, $\sqrt{13}=[3;1,1,1,1,6,\dots]$, and its convergents $[3;1]=3+\frac11=4$ and $[3;1,1]=3+\frac1{1+1/1}=\frac72$ are two winners of Section 1. Numerator and denominator follow the recurrence ^[2]

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

with the starting values $p_{-2}=0$, $p_{-1}=1$, $q_{-2}=1$ and $q_{-1}=0$. We use three facts about it. First, the recurrence evaluates $[a_0;a_1,\dots,a_k]$ for any positive real last entry, and $x=[a_0;a_1,\dots,a_{k-1},x_k]$, so replacing $a_k$ by $x_k$ gives $x$ itself:

$$x=\frac{x_kp_{k-1}+p_{k-2}}{x_kq_{k-1}+q_{k-2}}.$$

Second, each step of the recurrence changes the sign of $p_{k-1}q_{k-2}-p_{k-2}q_{k-1}$, which equals $1$ for $k=0$, so

$$p_{k-1}q_{k-2}-p_{k-2}q_{k-1}=(-1)^k,$$

and every convergent is in lowest terms. Third, $q_0=1\le q_1\lt q_2\lt\cdots$, and since $a_k\ge1$ we have $q_k\ge q_{k-1}+q_{k-2}$; hence $q_k$ is at least the Fibonacci number $F_{k+1}$, where $F_1=F_2=1$. For $\sqrt{13}$ the convergents begin $3/1$, $4/1$, $7/2$, $11/3$, $18/5$, $119/33$, $137/38$, $256/71$, in agreement with the OEIS. ^[6]

3. The continued fraction of $\sqrt n$ in whole numbers

For $x=\sqrt n$ every complete quotient has the form $x_k=(\sqrt n+m_k)/d_k$ with integers $m_k$ and $d_k$. Starting from $m_0=0$, $d_0=1$ and $a_0=\lfloor\sqrt n\rfloor$, they follow the rule ^[3]

$$m_{k+1}=d_ka_k-m_k,\qquad 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.$$

Indeed $x_k-a_k=(\sqrt n-m_{k+1})/d_k$, and multiplying the reciprocal by $\sqrt n+m_{k+1}$ above and below gives

$$\begin{aligned}x_{k+1}&=\frac{d_k\,(\sqrt n+m_{k+1})}{n-m_{k+1}^2}\\&=\frac{\sqrt n+m_{k+1}}{d_{k+1}}.\end{aligned}$$

The number $d_{k+1}$ is an integer because $d_k$ divides $n-m_k^2$ at every step. This holds for $k=0$; since $m_{k+1}\equiv-m_k\pmod{d_k}$, also $n-m_{k+1}^2\equiv n-m_k^2\equiv0\pmod{d_k}$; and $d_kd_{k+1}=n-m_{k+1}^2$ shows that $d_{k+1}$ divides $n-m_{k+1}^2$, which continues the induction.

Changing the sign of $\sqrt n$ turns $x_k$ into its conjugate $\bar x_k=(m_k-\sqrt n)/d_k$, and the conjugates obey the same rule $\bar x_{k+1}=1/(\bar x_k-a_k)$. They start with $\bar x_1=-1/(\sqrt n+a_0)$, which lies between $-1$ and $0$, and they stay there, because $\bar x_k-a_k\lt-1$ whenever $-1\lt\bar x_k\lt0$ and $a_k\ge1$. Together with $x_k\gt1$ this gives, for $k\ge1$, $d_k=2\sqrt n/(x_k-\bar x_k)\gt0$, $m_k=d_k(x_k+\bar x_k)/2\gt0$ and $m_k\lt\sqrt n$; then $a_k\ge1$ forces $d_k\le a_0+m_k$. In short, for every $k\ge1$, ^[3]

$$1\le m_k\le a_0,\qquad 1\le d_k\le 2a_0,\qquad 1\le a_k\le 2a_0.$$

Because $d_{k+1}\ge1$, the rule may use $a_0$ instead of $\sqrt n$: $\lfloor y/d\rfloor=\lfloor\lfloor y\rfloor/d\rfloor$ for every real $y$ and every integer $d\ge1$, and $\lfloor\sqrt n+m\rfloor=a_0+m$. For $n\le100000$ all three numbers are at most $632$, so the expansion needs only small integers and no square root beyond $a_0$. Numbers with $x_k\gt1$ and $-1\lt\bar x_k\lt0$ are called reduced; they make the expansion of $\sqrt n$ periodic, with each period ending in $2a_0$, ^[3] but the method does not need this. For $\sqrt{13}$ the triples $(m_k,d_k,a_k)$ for $k=0,\dots,5$ are $(0,1,3)$, $(3,4,1)$, $(1,3,1)$, $(2,3,1)$, $(1,4,1)$ and $(3,1,6)$; after that $(3,4,1)$ returns, and the partial quotients repeat with period $1,1,1,1,6$. For instance $x_1=(\sqrt{13}+3)/4=1.651\ldots$ and $x_5=\sqrt{13}+3=6.605\ldots$.

4. Intermediate fractions and their errors

Fix $k\ge1$. The intermediate fractions, also called semiconvergents, that lead up to the convergent $p_k/q_k$ are ^[2]

$$L_t=\frac{t\,p_{k-1}+p_{k-2}}{s_t},\qquad s_t=t\,q_{k-1}+q_{k-2},$$

for $t=0,1,\dots,a_k$, except that $t=0$ is left out when $k=1$, because then $s_0=q_{-1}=0$. The last one, $L_{a_k}$, is the convergent $p_k/q_k$, and for $k\ge2$ the first one, $L_0$, is $p_{k-2}/q_{k-2}$. For $\sqrt{13}$ and $k=5$ they are $11/3$, $29/8$, $47/13$, $65/18$, $83/23$, $101/28$ and $119/33$, the family of Section 1. When $x-L_t$ is brought to a common denominator, with $x$ written as in Section 2, its numerator collapses to $x_k-t$ times the left side of the determinant identity, so

$$x-L_t=\frac{(-1)^k\,(x_k-t)}{s_t\,(x_kq_{k-1}+q_{k-2})}.$$

The same computation with $p_{k-1}/q_{k-1}$ gives

$$x-\frac{p_{k-1}}{q_{k-1}}=\frac{(-1)^{k-1}}{q_{k-1}\,(x_kq_{k-1}+q_{k-2})}.$$

Since $t\le a_k\lt x_k$, the factor $x_k-t$ is positive: all the fractions $L_t$ lie on one side of $x$, and $p_{k-1}/q_{k-1}$ lies on the other. As $t$ grows, $x_k-t$ shrinks and $s_t$ grows, so $L_t$ moves steadily towards $x$ without reaching it. In the example, $18/5=3.6$ lies below $\sqrt{13}$, and the fractions $(18t+11)/(5t+3)$ lie above it and come closer as $t$ grows.

5. No small denominators between neighbours

Call two fractions $u/v\lt u'/v'$ with positive denominators neighbours if $vu'-uv'=1$. Every fraction $w/z$ strictly between two neighbours has $z\ge v+v'$. Indeed, $vw-uz$ and $u'z-v'w$ are positive integers, so

$$\begin{aligned}\frac1{vv'}&=\Bigl(\frac{u'}{v'}-\frac wz\Bigr)+\Bigl(\frac wz-\frac uv\Bigr)\\&\ge\frac1{v'z}+\frac1{vz},\end{aligned}$$

and the right side equals $(v+v')/(vv'z)$, which gives $z\ge v+v'$. This is the familiar property of neighbours in a Farey sequence, between which the first fraction to appear is the mediant $(u+u')/(v+v')$. ^[4]

The fractions $p_{k-1}/q_{k-1}$ and $L_t$ are neighbours, in one order or the other: by the determinant identity,

$$p_{k-1}s_t-q_{k-1}\,(t\,p_{k-1}+p_{k-2})=(-1)^k.$$

Therefore every $L_t$ is in lowest terms, and every fraction strictly between $p_{k-1}/q_{k-1}$ and $L_t$ has a denominator at least $q_{k-1}+s_t=s_{t+1}$. For example $18\cdot18-5\cdot65=-1$, and the only fractions between $18/5$ and $65/18$ with denominators up to $30$ are $83/23$ and $101/28$, both with denominators at least $5+18=23$.

6. Only two candidates

Theorem. Let $x$ be irrational and $D\ge1$ an integer, let $k$ be the first index with $q_k\gt D$, and put

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

Then the best approximation to $x$ for the bound $D$ is whichever of $p_{k-1}/q_{k-1}$ and $L_{t^*}$ is closer to $x$.

Proof. Since $q_0=1\le D$, we have $k\ge1$ and $q_{k-1}\le D$. By the choice of $t^*$, $s_{t^*}\le D\lt s_{t^*+1}$, and $t^*\ge0$ because $q_{k-2}\le q_{k-1}\le D$. From $q_k=s_{a_k}\gt D$ we get $t^*\lt a_k$, so $L_{t^*}$ is one of the fractions of Section 4. For $k=1$, the value $q_{-1}=0$ gives $t^*=D\ge1$, so $L_{t^*}=a_0+1/D$ is a genuine fraction. By Section 4, $x$ lies strictly between $L_{t^*}$ and $p_{k-1}/q_{k-1}$, and by Section 5 every fraction strictly between them has a denominator at least $s_{t^*+1}\gt D$. Hence any other fraction with denominator at most $D$ lies outside the closed interval with these two endpoints, beyond one of them, and that endpoint is strictly closer to $x$. Both endpoints are in lowest terms with denominators at most $D$, so the closer endpoint is the best approximation. $\square$

This is a sharp form of the classical fact that every best approximation in the sense of this problem, a best approximation of the first kind, is a convergent or an intermediate fraction. ^[2] ^[5] For $\sqrt{13}$ and $D=20$, the first denominator above $20$ is $q_5=33$, so $k=5$, $q_4=5$, $q_3=3$, $t^*=\lfloor17/5\rfloor=3$, and the candidates are $18/5$ and $65/18$. For $D=30$ we get $t^*=\lfloor27/5\rfloor=5$ and the candidates $18/5$ and $101/28$.

7. Which of the two wins

Dividing the two error formulas of Section 4 gives

$$\frac{|x-L_t|}{|x-p_{k-1}/q_{k-1}|}=\frac{q_{k-1}\,(x_k-t)}{s_t},$$

so $L_t$ is closer to $x$ exactly when $q_{k-1}(x_k-t)\lt s_t$, that is,

$$q_{k-1}\,x_k\lt s_t+t\,q_{k-1}=2t\,q_{k-1}+q_{k-2}.$$

This criterion has a simple reading, which we call the half rule, because it compares $t$ with $a_k/2$. ^[2] Write $x_k=a_k+1/x_{k+1}$ with $0\lt1/x_{k+1}\lt1$, and note that $0\le q_{k-2}/q_{k-1}\le1$. Divided by $q_{k-1}$, the criterion becomes $2t\gt a_k+1/x_{k+1}-q_{k-2}/q_{k-1}$, and the right side lies strictly between $a_k-1$ and $a_k+1$. So $L_t$ beats $p_{k-1}/q_{k-1}$ when $t\gt a_k/2$ and loses when $t\lt a_k/2$, in particular for $t=0$. Only for $t=a_k/2$ is a closer look needed: then $L_t$ wins exactly when $1/x_{k+1}\lt q_{k-2}/q_{k-1}$. For $k\ge2$ we have $q_{k-1}/q_{k-2}=[a_{k-1};a_{k-2},\dots,a_1]$, so this is the comparison of $[a_k;a_{k-1},\dots,a_1]$ with $[a_k;a_{k+1},\dots]$ in the usual statement of the rule; for $k=1$ the quotient $q_{-1}/q_0$ is $0$, and the half value always loses.

The checkpoint $D=20$ is exactly this borderline case, with $a_5=6$ and $t^*=3$. Here $x_6=(\sqrt{13}+3)/4$, so $1/x_6=\sqrt{13}-3=0.60555\ldots$, just above $q_3/q_4=3/5$. Therefore $65/18$ loses, and its error exceeds that of $18/5$ by only about $8.56\cdot10^{-6}$.

8. An exact test with integers

For $x=\sqrt n$, insert $x_k=(\sqrt n+m_k)/d_k$ into the criterion and multiply by $d_k\ge1$. The intermediate fraction $L_t$ wins exactly when

$$q_{k-1}\sqrt n\lt R,\qquad R=d_k\,(s_t+t\,q_{k-1})-q_{k-1}m_k.$$

Here $R$ is an integer, while $q_{k-1}\sqrt n$ is irrational, so the two sides are never equal, and three cases settle the comparison without any rounding:

if $R\le0$, then $L_t$ loses;
if $R\gt q_{k-1}(a_0+1)$, then $L_t$ wins, because $\sqrt n\lt a_0+1$;
otherwise both sides are positive, and $L_t$ wins exactly when $R^2\gt n\,q_{k-1}^2$.

9. The checkpoints, step by step

For $\sqrt{13}$ and $D=20$, Section 6 found $k=5$ and $t^*=3$. When the denominator $q_5=33$ appears, the state is $(m_5,d_5,a_5)=(3,1,6)$, so with $s_3+3q_4=18+15=33$ we get

$$R=1\cdot33-5\cdot3=18.$$

Since $0\lt18\le5\cdot4=20$, the square test decides: $R^2=324$ is smaller than $13\cdot5^2=325$. Thus $18/5$ wins, as the statement says, and the whole decision rests on $324$ against $325$. For $D=30$ we get $t^*=5$, $s_5+5q_4=28+25=53$ and $R=53-15=38\gt20$, so $101/28$ wins at once, as stated.

How the Code Works

The C++ program ^[9] fixes kDenominatorBound $=10^{12}$ and kMaxN $=100000$. isqrt_u64 corrects a long double square root with exact integer tests, and is_square uses it to skip the perfect squares.

best_denominator_cf(n, d_bound) carries out Sections 2 to 8 in one loop. It keeps the state m, den, a of Section 3, starting from $(0,1,a_0)$, and the last two convergents p_prev1/q_prev1 and p_prev2/q_prev2, starting from $1/0$ and $0/1$, the values for the indices $-1$ and $-2$. Each pass forms the next convergent p/q by the recurrence of Section 2. While q fits, the convergents shift by one place and the state advances through m_next, den_next and a_next. When q first exceeds the bound, q_prev1 and q_prev2 are $q_{k-1}$ and $q_{k-2}$, and m, den, a are $m_k$, $d_k$, $a_k$, exactly the quantities of Section 6.

The function then sets best_q to $q_{k-1}$ and replaces it only if the intermediate fraction wins. If d_bound exceeds q_prev1, the integer division (d_bound - q_prev2) / q_prev1 gives t $=t^*$. If t is positive, the function forms q_semiconv $=s_{t^*}$ and rhs $=R$, computed as $d_k(2tq_{k-1}+q_{k-2})-q_{k-1}m_k$, and applies the three cases of Section 8 with upper $=q_{k-1}(a_0+1)$, lhs_sq $=nq_{k-1}^2$ and rhs_sq $=R^2$. Both shortcuts that keep $q_{k-1}$ without a test are correct. For $t^*=0$ the other candidate, $p_{k-2}/q_{k-2}$, loses by Section 7. If the bound equals $q_{k-1}$, then $t^*=0$ unless $k=1$, and in that case $D=1$ and both candidates have the denominator $1$. Only the denominator is returned, which is all the sum needs.

The state stays below $633$ by Section 3, so 64-bit words are ample for it. The convergents are unsigned 128-bit integers, a GCC and Clang extension, ^[8] although the largest values that occur, about $1.9\cdot10^{17}$ for $p_k$ and $6.2\cdot10^{14}$ for $q_k$, would fit in 64 bits. The test of Section 8 does need more: $|R|$ stays below $3\cdot10^{14}$, but $R^2$ and $nq_{k-1}^2$ reach $96$ and $97$ bits, so they use the signed 128-bit type.

run_validations checks best_denominator_cf(13, 20) $=5$ and best_denominator_cf(13, 30) $=28$, and compares the function with best_denominator_bruteforce for every non-square $n\le200$ and every bound from $2$ to $200$. The brute force tries every denominator $q$ with the numerators $\lfloor q\sqrt n\rfloor$ and $\lfloor q\sqrt n\rfloor+1$, skips fractions that are not in lowest terms and compares errors in long double. main stops if a check fails, adds best_denominator_cf(n, kDenominatorBound) over the non-square $n\le$ kMaxN in unsigned 128-bit integers, and prints the total with to_string_u128, since the standard streams cannot print that type.

The Python and Java versions skip the validations and make the same decisions under slightly different variable names. Python's solve runs the same loop with math.isqrt and unlimited-precision integers. Java computes $a_0$ with isqrt, a double square root corrected by the same integer tests as isqrt_u64, keeps m, d and a in long, and uses BigInteger for the convergents, t, rhs, the squares and the total. ^[8]

Complexity and Verification

For one $n$, each pass of the loop costs a fixed number of integer operations, and there are at most $60$ passes: by Section 2, $q_k\ge F_{k+1}$, and $F_{60}=1548008755920$ exceeds $10^{12}$, so the first denominator above the bound has an index $k\le59$. Over all $99684$ values of $n$ the loop makes $2418093$ passes, $24.26$ per value on average and at most $44$, reached only for $n=3$ with $\sqrt3=[1;1,2,1,2,\dots]$. The final test costs a constant number of operations, so for $N$ values of $n$ and the bound $D$ the method uses $O(N\log D)$ arithmetic operations on integers of $O(\log N+\log D)$ bits and a constant number of variables; on our test machine each of the three programs, the C++ one including its checks, finishes in well under a second. Trying all $10^{12}$ denominators for each $n$, by contrast, would take about $10^{17}$ steps.

Before the main computation the C++ program runs the checks described above, and they pass. For this article we also verified, with separate programs:

For every non-square $n\le1000$ and every bound $D\le1000$, $969000$ pairs in all, the nearest fraction found by trying all denominators, with errors compared to $50$ significant digits, is the candidate chosen by the method; the intermediate fraction wins in $202309$ of these pairs.
The bounds of Section 3 hold over a full period of the expansion of every non-square $n\le100000$, and $d_k$ always divides $n-m_{k+1}^2$.
At full size the intermediate fraction wins for $19815$ values of $n$ and the convergent for $79869$. As the half rule predicts, the intermediate fraction always wins when $t^*\gt a_k/2$, in $15321$ cases, and always loses when $t^*\lt a_k/2$, in $74500$ cases. The borderline case $t^*=a_k/2$ occurs $9863$ times, with $4494$ wins, and the square test is needed $23572$ times.
Python's Fraction.limit_denominator, which returns the nearest fraction with a denominator at most a given bound, ^[7] applied to a rational number within $10^{-80}$ of $\sqrt n$, gives the same fraction as the method for all $99684$ values of $n$. This check is valid because the distances from $\sqrt n$ of two different fractions with denominators at most $10^{12}$ differ by more than $7\cdot10^{-52}$. On the same side of $\sqrt n$, the difference is the distance between the fractions, at least $10^{-24}$. On opposite sides, it is twice the distance from $\sqrt n$ to their midpoint $u/v$, where $v\le2\cdot10^{24}$; that distance is at least $1/(v^2(2\sqrt n+1))$, because $u^2-nv^2$ is a nonzero integer and, when the distance is below $1$, $u/v+\sqrt n\lt2\sqrt n+1$. Moving $\sqrt n$ by less than $10^{-80}$ changes each such difference by less than $2\cdot10^{-80}$, so the nearest fraction stays the same.
The C++ program built with GCC and with Clang, the Java program, the Python program and our own transcription of the method print the same total.

Footnotes and References

The references below supply the background facts used above. The two-candidate theorem in the form used here, the integer test and all the counts are derived or computed explicitly in this article.

Project Euler 192 — Best Approximations. The official statement defines the best approximation to x for the denominator bound d as a reduced fraction r/s with s ≤ d such that every rational number closer to x has a denominator larger than d, gives 18/5 and 101/28 for √13 with the bounds 20 and 30, and asks for the sum of the denominators of the best approximations to √n for the bound 10¹² over the non-square n with 1 < n ≤ 100000.
Wikipedia — Simple continued fraction. The convergents satisfy hₙ = aₙhₙ₋₁ + hₙ₋₂ and kₙ = aₙkₙ₋₁ + kₙ₋₂ with h₋₁ = 1, h₋₂ = 0, k₋₁ = 0, k₋₂ = 1, and kₙhₙ₋₁ − kₙ₋₁hₙ = (−1)ⁿ; the semiconvergents (hₙ₋₁ + m·hₙ)/(kₙ₋₁ + m·kₙ), 0 ≤ m ≤ aₙ₊₁, form a monotonic sequence between the convergents. A best rational approximation is defined as a fraction closer to x than any approximation with a smaller or equal denominator, and all of them arise by truncating the continued fraction and reducing its last term to no less than half its value; when aₖ is even, the halved term aₖ/2 is admissible if and only if [aₖ; aₖ₋₁, …, a₁] > [aₖ; aₖ₊₁, …].
Wikipedia — Periodic continued fraction. For non-square S, the expansion of √S follows from m₀ = 0, d₀ = 1, a₀ = ⌊√S⌋ by mₙ₊₁ = dₙaₙ − mₙ, dₙ₊₁ = (S − mₙ₊₁²)/dₙ and aₙ₊₁ = ⌊(a₀ + mₙ₊₁)/dₙ₊₁⌋; it is periodic, and each period ends with 2a₀. A quadratic surd is reduced if it exceeds 1 and its conjugate lies between −1 and 0; by a theorem of Galois, exactly the reduced surds have purely periodic expansions.
Wikipedia — Farey sequence. If a/b < c/d are neighbours in a Farey sequence, then bc − ad = 1, and the first term that appears between them as the order of the sequence grows is the mediant (a + c)/(b + d), which first appears in the Farey sequence of order b + d.
Wikipedia — Diophantine approximation. A best approximation in the first sense, the notion of this problem, is closer to α than every other fraction with a smaller or equal denominator. The article explains that continued fractions yield these best approximations once the semiconvergents are considered along with the convergents.
OEIS A010122 — Continued fraction for sqrt(13), with A041018 and A041019, the numerators and denominators of its convergents. The partial quotients are 3, 1, 1, 1, 1, 6, 1, 1, 1, 1, 6, …, the numerators 3, 4, 7, 11, 18, 119, 137, 256, … and the denominators 1, 1, 2, 3, 5, 33, 38, 71, …
Python — fractions. Fraction.limit_denominator(max_denominator) returns the closest Fraction with a denominator at most max_denominator.
GCC manual — 128-bit Integers: __int128 and unsigned __int128 exist on targets with a 128-bit integer mode. Java — BigInteger: immutable arbitrary-precision integers. Python — numeric types: integers have unlimited precision.
C++ — Euler192.cpp. The linked immutable revision contains best_denominator_cf, best_denominator_bruteforce and run_validations; it is the implementation record for this article.

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);
    }
}