xi. Continued Fractions and Pell Equations

Every real $x$ has a continued fraction $[a_0;a_1,a_2,\dots ]$ ($a_0=\lfloor x\rfloor$, then recurse on $1/(x-a_0)$); for rationals it is Euclid’s quotient sequence. The convergents $p_k/q_k$ ($p_k=a_k{p}_{k-1}+p_{k-2}$, same for $q$) are the best rational approximations, and they are exactly the turning points of the Stern–Brocot tree, the binary search tree containing every positive rational once. Quadratic irrationals have periodic expansions, which solve Pell’s equation $x^2-D{y}^2=1$.

Task	Tool	Cost
convergents, best approximation with denominator $\le N$	continued fraction / Stern–Brocot	$O(\log N)$
smallest fraction with a monotone property, $p,q\le N$	binary search on the Stern–Brocot tree with run-length jumps	$O({\log }^{2}N)$ predicate calls
simplest fraction strictly between two fractions	Stern–Brocot LCA / continued-fraction comparison	$O(\log )$
$x^2-D{y}^2=1$	period of $\sqrt{D}$; fundamental solution generates all	$O(\text{period})$ (period $O(\sqrt{D}\log D)$)
$x^2-D{y}^2=-1$	solvable iff the period length is odd	same

The library (all with exact i64 arithmetic):

vector<i64> a = continuedFraction(p, q);        // p / q = a0 + 1 / (a1 + 1 / (a2 + ...)), q > 0
vector<Frac> c = convergents(a);               // c[i] = [a0; a1, ..., ai] in lowest terms
auto [lo, hi] = boundingFractions(p, q, N);    // lo <= p / q <= hi, denominators <= N, closest on each side
Frac b = bestApproximation(p, q, N);           // closest to p / q with denominator <= N (lower on ties)
auto [no, yes] = fracBinarySearch([&](Frac x) { return x.p * x.p >= 2 * x.q * x.q; }, N);
                                               // yes = smallest fraction with p, q <= N and f true
auto path = sternBrocotEncode(p, q);          // runs {('L' | 'R'), count} from the root 1/1
Frac f = sternBrocotDecode(path);
Frac l = sternBrocotLca(p1, q1, p2, q2);
optional<Frac> a = sternBrocotAncestor(p, q, depth);   // depth of 1/1 is 0
auto [lo, hi] = sternBrocotRange(p, q);       // descendants = fractions strictly between lo and hi

struct Frac {
    i64 p, q;
};
 
i64 floorDiv(i64 a, i64 b) { return a / b - ((a % b != 0) && ((a < 0) != (b < 0))); }
 
vector<i64> continuedFraction(i64 p, i64 q) {
    vector<i64> a;
    while (q != 0) {
        i64 t = floorDiv(p, q);
        a.push_back(t);
        p -= t * q;
        swap(p, q);
    }
    return a;
}
 
vector<Frac> convergents(const vector<i64>& a) {
    vector<Frac> res;
    Frac prev{1, 0}, prev2{0, 1};
    for (i64 x : a) {
        Frac cur{x * prev.p + prev2.p, x * prev.q + prev2.q};
        res.push_back(cur);
        prev2 = prev;
        prev = cur;
    }
    return res;
}
 
pair<Frac, Frac> boundingFractions(i64 p, i64 q, i64 N) {
    Frac prev{1, 0}, prev2{0, 1};
    while (true) {
        i64 a = floorDiv(p, q);
        if (prev.q > 0 && a > (N - prev2.q) / prev.q) {
            i64 t = (N - prev2.q) / prev.q;
            Frac c = prev, s{t * prev.p + prev2.p, t * prev.q + prev2.q};
            if ((i128)c.p * s.q > (i128)s.p * c.q) swap(c, s);
            return {c, s};
        }
        Frac cur{a * prev.p + prev2.p, a * prev.q + prev2.q};
        prev2 = prev;
        prev = cur;
        p -= a * q;
        swap(p, q);
        if (q == 0) return {cur, cur};
    }
}
 
Frac bestApproximation(i64 p, i64 q, i64 N) {
    auto [lo, hi] = boundingFractions(p, q, N);
    // lo and hi are adjacent, so below * hi.q + above * lo.q = q: no overflow.
    i128 below = (i128)p * lo.q - (i128)lo.p * q, above = (i128)hi.p * q - (i128)p * hi.q;
    return below * hi.q <= above * lo.q ? lo : hi;
}
 
template <class F>
pair<Frac, Frac> fracBinarySearch(F f, i64 N) {
    Frac lo{0, 1}, hi{1, 0};
    for (int dir = 0, stuck = 0; stuck < 2; dir ^= 1) {
        Frac& a = dir ? hi : lo;   // dir 0 moves lo towards hi, dir 1 moves hi towards lo
        Frac b = dir ? lo : hi;
        i64 tmax = LLONG_MAX;
        if (b.p > 0) tmax = min(tmax, (N - a.p) / b.p);
        if (b.q > 0) tmax = min(tmax, (N - a.q) / b.q);
        auto ok = [&](i64 t) { return t <= tmax && f(Frac{a.p + t * b.p, a.q + t * b.q}) == (dir == 1); };
        i64 t = 0, step = 1;
        while (ok(t + step)) {
            t += step;
            step *= 2;
        }
        while (step > 1) {
            step /= 2;
            if (ok(t + step)) t += step;
        }
        a = Frac{a.p + t * b.p, a.q + t * b.q};
        stuck = t > 0 ? 0 : stuck + 1;
    }
    return {lo, hi};
}
 
// Stern-Brocot tree (root 1/1): the path to p/q as runs of 'L' (smaller) / 'R' (larger) moves.
vector<pair<char, i64>> sternBrocotEncode(i64 p, i64 q) {   // p, q >= 1, gcd 1
    vector<i64> a = continuedFraction(p, q);
    a.back()--;
    vector<pair<char, i64>> path;
    for (int i = 0; i < (int)a.size(); i++) {
        if (a[i] > 0) path.push_back({i % 2 ? 'L' : 'R', a[i]});
    }
    return path;
}
 
// walks a path from the root; returns {node, left bound, right bound} (descendants lie strictly between).
array<Frac, 3> sternBrocotWalk(const vector<pair<char, i64>>& path) {
    Frac lo{0, 1}, hi{1, 0};
    for (auto [c, k] : path) {
        if (c == 'R') lo = Frac{lo.p + k * hi.p, lo.q + k * hi.q};
        else hi = Frac{hi.p + k * lo.p, hi.q + k * lo.q};
    }
    return {Frac{lo.p + hi.p, lo.q + hi.q}, lo, hi};
}
Frac sternBrocotDecode(const vector<pair<char, i64>>& path) { return sternBrocotWalk(path)[0]; }
 
i64 sternBrocotDepth(i64 p, i64 q) {
    i64 d = 0;
    for (auto [c, k] : sternBrocotEncode(p, q)) d += k;
    return d;
}
 
Frac sternBrocotLca(i64 p1, i64 q1, i64 p2, i64 q2) {
    auto a = sternBrocotEncode(p1, q1), b = sternBrocotEncode(p2, q2);
    vector<pair<char, i64>> common;
    for (int i = 0; i < (int)min(a.size(), b.size()); i++) {
        if (a[i].first != b[i].first) break;
        common.push_back({a[i].first, min(a[i].second, b[i].second)});
        if (a[i].second != b[i].second) break;
    }
    return sternBrocotDecode(common);
}
 
optional<Frac> sternBrocotAncestor(i64 p, i64 q, i64 depth) {   // the ancestor at depth (root 1/1 = 0)
    vector<pair<char, i64>> prefix;
    for (auto [c, k] : sternBrocotEncode(p, q)) {
        if (depth == 0) break;
        prefix.push_back({c, min(k, depth)});
        depth -= min(k, depth);
    }
    if (depth > 0) return nullopt;
    return sternBrocotDecode(prefix);
}
 
pair<Frac, Frac> sternBrocotRange(i64 p, i64 q) {   // subtree of p/q = fractions strictly in (first, second)
    auto w = sternBrocotWalk(sternBrocotEncode(p, q));
    return {w[1], w[2]};
}

Pell’s equation

vector<ll> cfSqrt(ll D) {                                  // [a0; period...] of sqrt(D), empty if D is a square
    ll a0 = (ll)sqrtl((long double)D); while (a0 * a0 > D) a0--; while ((a0 + 1) * (a0 + 1) <= D) a0++;
    if (a0 * a0 == D) return {};
    vector<ll> cf{a0}; ll m = 0, d = 1, a = a0;
    while (a != 2 * a0) { m = d * a - m; d = (D - m * m) / d; a = (a0 + m) / d; cf.push_back(a); }
    return cf;                                             // cf[1..] is one full period (it ends with 2 a0)
}
pair<i128, i128> pellFundamental(ll D) {                   // least x, y > 0 with x^2 - D y^2 = 1 (values may need i128)
    auto cf = cfSqrt(D); int r = sz(cf) - 1;               // period length
    i128 p0 = 1, q0 = 0, p1 = cf[0], q1 = 1;
    for (int i = 1; i < 2 * r; i++) {                      // the (r-1)-th convergent (or (2r-1)-th if r is odd) solves it
        if (p1 * p1 - (i128)D * q1 * q1 == 1) return {p1, q1};
        ll a = cf[(i - 1) % r + 1];
        i128 p2 = a * p1 + p0, q2 = a * q1 + q0; p0 = p1, q0 = q1, p1 = p2, q1 = q2;
    }
    return {p1, q1};
}
// all solutions: x_k + y_k sqrt(D) = (x_1 + y_1 sqrt(D))^k, i.e. x_{k+1} = x_1 x_k + D y_1 y_k, y_{k+1} = x_1 y_k + y_1 x_k

Fact	Statement
convergent recurrences	$p_k=a_k{p}_{k-1}+p_{k-2}$, $q_k=a_k{q}_{k-1}+q_{k-2}$, $p_k{q}_{k-1}-p_{k-1}{q}_k=(-1{)}^{k-1}$
approximation quality	$\left
best approximations	every best approximation of the second kind is a convergent; of the first kind, a convergent or an intermediate fraction $\frac{p_{k-1}+t{p}_k}{q_{k-1}+t{q}_k}$
Stern–Brocot	mediant insertion; the path to $p/q$ is $R^{a_0}{L}^{a_1}{R}^{a_2}\cdots$ (last run shortened by one)
quadratic irrationals	eventually periodic expansions (Lagrange); $\sqrt{D}=[a_0;\overline{a_1,\dots ,a_{r-1},2a_0}]$ with a palindromic body
Pell	fundamental solution is the convergent $p_{r-1}/q_{r-1}$ ($r$ even) or $p_{2r-1}/q_{2r-1}$ ($r$ odd); $x^2-D{y}^2=-1$ solvable iff $r$ odd
Farey neighbours	$\frac{a}{b}<\frac{c}{d}$ adjacent in a Farey sequence iff $bc-ad=1$; the mediant $\frac{a+c}{b+d}$ is the simplest fraction between them

Proof: convergents, the determinant identity, and approximation bounds

The recurrences follow by induction from $[a_0;\dots ,a_k,t]$ being a Möbius transformation $\frac{t{p}_k+p_{k-1}}{t{q}_k+q_{k-1}}$ in its last argument. The determinant identity: $p_k{q}_{k-1}-p_{k-1}{q}_k=-(p_{k-1}{q}_{k-2}-p_{k-2}{q}_{k-1})$, with initial value $1$. Hence consecutive convergents differ by $\frac{(-1{)}^{k-1}}{q_k{q}_{k-1}}$, they alternate around $x$, and $x$ lies between $p_k/q_k$ and $p_{k+1}/q_{k+1}$, giving $|x-p_k/q_k|<1/(q_k{q}_{k+1})$. Legendre’s criterion: if $|x-p/q|<1/(2q^2)$ but $p/q$ is not a convergent, place $q$ between $q_k\le q<q_{k+1}$; then $|p/q-p_k/q_k|\ge 1/(q{q}_k)$ while both are within $1/(2q^2)$ and $1/(q_k{q}_{k+1})$ of $x$, which leads to a contradiction (standard algebra). Stern–Brocot: the tree inserts mediants between neighbours with determinant 1, which preserves determinant 1 and lowest terms and reaches every positive rational; descending toward $p/q$ moves right while $p/q$ exceeds the current node, and the runs of equal moves are the continued-fraction digits.

Proof sketch: Pell's equation

The convergents of $\sqrt{D}$ satisfy $p_k^2-D{q}_k^2=(-1{)}^{k+1}{Q}_{k+1}$ where $Q$ are the denominators $d$ of the expansion algorithm, which equal 1 exactly at the ends of periods. So $p_{r-1}^2-D{q}_{r-1}^2=(-1{)}^r$: a solution of $x^2-D{y}^2=1$ for even $r$, of $-1$ for odd $r$ (then squaring, the convergent at index $2r-1$ solves $+1$). Every solution with $x,y>0$ satisfies $|x/y-\sqrt{D}|<1/(2y^2)$, so it is a convergent (Legendre), which is why the first period end gives the smallest one. The solutions form a cyclic group under $(x+y\sqrt{D})(x^{\prime }+y^{\prime }\sqrt{D})$: if $u$ were a solution between two consecutive powers of the fundamental $\epsilon$, then $u{\epsilon }^{-k}$ would be a solution with $1<u{\epsilon }^{-k}<\epsilon$, contradicting minimality. Solutions grow exponentially, so they overflow 64 bits quickly (use i128 or big integers).