Continued fractions of rational numbers

The continued fraction of a rational number is the sequence of quotients of the Euclidean algorithm. FLINT computes it in quasi-linear time by the half-GCD method [Sch1971]. For the basic theory see [HarWri2008, Chapter X] and [GKP1994, §6.7].

Continued fractions and matrices#

For x = p/q with q > 0, the Euclidean algorithm r_{-1} = p, r_0 = q, r_{i-1} = c_i r_i + r_{i+1} with c_i = \lfloor r_{i-1}/r_i \rfloor gives the expansion x = [c_0; c_1, \ldots, c_{k-1}], with c_0 = \lfloor x \rfloor and c_i \ge 1 for i \ge 1, and c_{k-1} \ge 2 if k \ge 2 (the shorter of the two expansions). Writing

the fractions p_i/q_i are the convergents of x, and for every y > 0

So the first i terms of x are a matrix M of determinant (-1)^i with nonnegative entries (apart from the first row when c_0 < 0), and the rest of the expansion is that of y = M^{-1}(x), where M acts as a Möbius transformation. The entries of M have about as many bits in total as the terms consumed, and y has correspondingly fewer bits than x.

Terms from leading bits#

The key observation, due to Lehmer [Leh1938], is that the first terms of x depend only on its leading bits.

If x = a/b with a > b > 0 of n bits and a', b' are the leading k bits of a and b (so a = 2^{n-k} a' + \alpha, 0 \le \alpha < 2^{n-k}, and similarly for b), then x lies in the interval (a'/(b' + 1), (a'+1)/b'), whose endpoints have O(k) bits. Expanding both endpoints in parallel until they disagree gives about k/2 bits worth of terms of x at the cost of a computation on k-bit numbers. Applying the resulting matrix M^{-1} to the full (a, b) then reduces it to a pair of about n - k/2 bits.

The half-GCD algorithm#

Each level costs a constant number of multiplications of n-bit numbers, so the cost is O(M(n) \log n), against O(n^2) for the Euclidean algorithm. The half-GCD algorithm for integer GCDs, which this one parallels, is described on Integer GCD. As with all half-GCD algorithms, the delicate part is the correction when the leading bits produce a term that is not a term of the full number; FLINT avoids it by working with intervals throughout (Lemma 1), so that every term written is guaranteed.

Evaluation#

Evaluating [c_0; \ldots, c_{n-1}] amounts to computing the matrix product E_{c_0} \cdots E_{c_{n-1}}. Multiplying the factors in order costs O(nb) for a b-bit result, since the partial products grow linearly. Splitting the product in halves recursively balances the sizes of the operands, and the cost becomes O(M(b) \log n), as in Binary splitting and product trees.

The simplest fraction in an interval#

The Stern–Brocot tree contains every positive rational once: its root is 1/1, and the children of a node are obtained as mediants (a+c)/(b+d) of the node and its nearest ancestors a/b < c/d [GKP1994, §4.5]. A path in the tree is a continued fraction expansion, read as run lengths of left and right moves.

In terms of continued fractions, the common terms of l and r are the common part of the two paths. FLINT computes them with the interval version of the half-GCD algorithm, which handles a ball [l, r] directly (_fmpq_ball_get_cfrac), and then completes the expansion with the smallest integer in the transformed interval. Negative and mixed-sign intervals are reduced to this case by subtracting an integer or noting that 0 or the smallest integer in the interval is the answer.

Farey neighbours#

Two fractions a/b < c/d with bc - ad = 1 are adjacent in every Farey sequence of order Q with \max(b, d) \le Q < b + d. Given x = p/q and Q \ge q, the left neighbour p_l/q_l is therefore determined by p q_l \equiv 1 \pmod q and Q - q < q_l \le Q, which fmpq_farey_neighbors solves with one modular inverse; the right neighbour follows from the mediant relation p_l/q_l + p_r/q_r = v \cdot p/q with v = \lfloor (Q + q_l)/q \rfloor [HarWri2008, Chapter III].

Dedekind sums#

By the reciprocity law of Dedekind sums [RadGro1972], s(h, k) is a simple function of the continued fraction of k/h: with the partial quotients c_1, \ldots, c_n and the convergent matrix M = (m_{ij}) of k/(h \bmod k),

[Hickerson1977]. FLINT evaluates it with word-size quotients when k fits in a word and otherwise with the half-GCD code, which accumulates the alternating sum of the terms instead of storing them. The cost is that of a GCD.

Use in FLINT#

• fmpq_get_cfrac computes the terms with _fmpq_ball_get_cfrac, the interval half-GCD (with _fmpq_hgcd as its recursive kernel and Lehmer's method on one and two words below), then finishes with division steps; fmpq_get_cfrac_naive is the Euclidean algorithm.
• fmpq_set_cfrac multiplies the elementary matrices in a balanced tree, with plain products below 32 terms.
• fmpq_simplest_between, fmpq_dedekind_sum and the large-modulus path of rational reconstruction (_fmpq_reconstruct_fmpz_2, see Rational reconstruction) use the same half-GCD code.
• fmpq_farey_neighbors uses the modular inverse formula above.