Integer GCD

The greatest common divisor of two integers is computed by the Euclidean algorithm or one of its accelerations. For n-limb operands the Euclidean algorithm takes O(n \mathtt{FLINT\_BITS}) division steps of decreasing size and O(n^2) limb operations in total; Lehmer's algorithm keeps this complexity but replaces almost all multiprecision steps by word operations; and the half-GCD algorithm of Knuth and Schönhage reduces the cost to O(M(n) \log n). The same algorithms compute continued fractions, Bézout coefficients and rational reconstructions, which are all byproducts of the sequence of quotients. The polynomial analogue is described in Half-GCD.

The Euclidean algorithm as a product of matrices#

Starting from r_0 = a > r_1 = b > 0, the Euclidean algorithm computes r_{i+1} = r_{i-1} - q_i r_i with q_i = \lfloor r_{i-1}/r_i \rfloor until r_{k+1} = 0, and then r_k = \gcd(a, b). In matrix form,

The matrix M_i has nonnegative entries and determinant (-1)^i, so its inverse is again an integer matrix, and its entries are the numerators and denominators of the convergents of the continued fraction a/b = [q_1; q_2, \ldots, q_k]. The second row of M_i^{-1} gives the cofactors s_i, t_i with s_i a + t_i b = r_i; at i = k these are Bézout coefficients. Each accelerated algorithm computes a product of many quotient matrices at once from leading parts of the operands, and applies it to the full operands.

The quotients that leading parts determine are characterised by a statement about continued fractions of real numbers. For a matrix M as above, write M \cdot y = (m_{11} y + m_{12})/(m_{21} y + m_{22}).

Thus quotients can be computed from the top bits of a and b, as long as the computation checks that the quotients found are valid for the whole interval; by Lemma 1, it suffices to check that the transformed endpoints stay above 1.

Binary GCD#

Stein's binary algorithm [Ste1967] avoids divisions altogether: it removes common factors of 2, then repeatedly replaces the larger of two odd numbers by their difference with its factors of 2 removed. For single words this is the fastest method on current processors, since a subtraction, a count of trailing zeros and a shift are cheap. Stehlé and Zimmermann [SZ2004] gave a subquadratic "binary recursive" variant that works from the least significant bits; it is not used in FLINT.

Lehmer's algorithm#

Lehmer [Leh1938] observed that the first quotients of a/b are almost always determined by the leading word of each operand. His algorithm runs the Euclidean algorithm on the leading words \hat a, \hat b (taken from the same bit position of a and b), accumulating the quotient matrix in single words, stops as soon as a quotient might differ from that of the full numbers, and then applies the matrix to a and b with a few multiply-and-add passes. Each pass reduces the operands by about one word at a cost of O(n), so the total is O(n^2) with a small constant.

Collins and Jebelean [Jeb1993b] showed that, with the remainders \hat r_i and cofactors u_i, v_i of the single-word computation, the quotients q_1, \ldots, q_{i} are those of a/b if

(with u in place of v for the other parity of i). This is the interval condition of Corollary 2 written in terms of the cofactors. Using two words of each operand ("double-digit" Lehmer), as GMP does in mpn_hgcd2, roughly doubles the progress per pass.

The half-GCD algorithm#

The half-GCD of a > b of n limbs is a matrix M such that (a, b)^T = M (a', b')^T, where a', b' are consecutive remainders with about n/2 limbs; the entries of M then have about n/2 limbs too. Knuth [Knu1970] and Schönhage [Sch1971] observed that it can be computed recursively from the top halves of a and b.

Each level does two recursive calls of half size and a constant number of products of n/2-limb numbers and matrices, so the cost is H(n) = 2H(n/2) + O(M(n)), which is O(M(n) \log n). The GCD is obtained by applying the half-GCD repeatedly, each time halving the size, at the same total cost. The correctness of the truncations in steps 2 and 4, which is delicate in Schönhage's original formulation, was simplified by Möller [Mol2008], whose version is the one in GMP. FLINT's own version formulates it with intervals as in Corollary 2, which makes the check of step 2 a comparison of the transformed endpoints with 1.

Use in FLINT#

fmpz_gcd, fmpz_gcd_ui, fmpz_gcd3. An inline operand reduces the problem to GMP's mpn_gcd_1, which divides the large operand by the small one once and finishes with a binary GCD on words. Two large operands go to GMP's mpz_gcd, that is mpn_gcd: Lehmer's algorithm with double-word steps (mpn_hgcd2) below GMP's threshold GCD_DC_THRESHOLD (several hundred limbs) and Möller's half-GCD above [GMP].

fmpz_xgcd, fmpz_xgcd_canonical_bezout, fmpz_gcdinv, fmpz_invmod. GMP's mpz_gcdext and mpz_invert, built on the same hgcd machinery with cofactor tracking; for word-size operands, GMP's mpn_gcdext_1 (fmpz_xgcd_canonical_bezout) or FLINT's n_gcdinv (fmpz_gcdinv, fmpz_invmod) (GCD, LCM and Bézout coefficients).

n_gcd, n_xgcd, n_gcdinv. Word GCDs: n_gcd calls GMP's binary mpn_gcd_11; n_xgcd and n_gcdinv are FLINT's own extended Euclidean algorithms with special handling of small quotients.

fmpz_xgcd_partial. FLINT's own Lehmer algorithm with the Jebelean condition of Lehmer's algorithm, which stops when the remainder falls below a bound L and returns the last two remainders and cofactors. It is the partial Euclidean algorithm of Shanks's NUCOMP and NUDUPL for binary quadratic forms in qfb, and is also used by Cornacchia's algorithm in ecpp.

_fmpq_hgcd and _fmpq_ball_get_cfrac. FLINT's own integer half-GCD in the interval formulation of Lemma 1: it maintains a matrix M and a pair (x_a, x_b) representing the interval M \cdot (x_a/(x_b + 1), (x_a + 1)/x_b), and generates quotients valid for the whole interval. It uses single Euclidean steps when the remaining size is below 4 limbs, a Lehmer algorithm with a two-limb half-GCD (_uiui_hgcd) up to 500 limbs, and above that recursion on the top halves as in Algorithm 1. It computes continued fractions (fmpq_get_cfrac), the simplest rational between two numbers (fmpq_simplest_between), Dedekind sums (fmpq_dedekind_sum) and rational reconstruction (fmpq_reconstruct_fmpz_2, which also has its own Lehmer loop for the final steps). Continued fractions of rational numbers and Rational reconstruction describe these applications.