Modular GCD

The modular algorithm computes the gcd of two polynomials over the integers without ever handling numbers much larger than those in the answer^[1][2]. Remainder-sequence methods, even the subresultant PRS, pass through intermediate coefficients whose size grows linearly with the degree, although the gcd itself usually has small coefficients. The modular algorithm reduces both inputs modulo a prime $p$, where the coefficients are single machine words and the gcd can be computed by Euclid's algorithm in the field ${\mathbb{F}}_p$. The images for several primes are combined by the Chinese remainder theorem into the gcd modulo their product. A bound on the size of the gcd's coefficients, from the Landau–Mignotte inequality, says when the product is large enough to determine them; a trial division then confirms the answer. W. S. Brown published the algorithm in 1971^[1].

Let $a$ and $b$ be primitive, with gcd $h$, and let $\gamma =\gcd (\mathrm{lc}(a),\mathrm{lc}(b))$. The leading coefficient of $h$ divides $\gamma$, so $\bar{h}=\frac{\gamma }{\mathrm{lc}(h)}h$ is a multiple of $h$ with the known leading coefficient $\gamma$, and its coefficients are bounded by

by the Landau–Mignotte bound^[6]. For each prime $p$ not dividing $\gamma$, compute the monic gcd of $a\mathrm{mod}p$ and $b\mathrm{mod}p$ in ${\mathbb{F}}_p[x]$ and scale it to leading coefficient $\gamma \mathrm{mod}p$. Combine the images with the Chinese remainder theorem, coefficient by coefficient. Once the product $M$ of the primes used exceeds $2B$, the representatives of the coefficients in $(-M/2,M/2]$ are the true integer coefficients of $\bar{h}$, and its primitive part is $h$^[1][2].

(* The modular gcd of two polynomials over Z: gcds modulo many primes,
   combined by the Chinese remainder theorem until the modulus exceeds
   twice a bound on the answer's coefficients, then checked by division. *)
open Zpoly
 
let primes_used = ref 0
 
(* Combine x mod m with y mod p into the residue mod m p. *)
let crt m (x : t) p (y : int array) =
  let minv = Modp.inv p (Z.to_int (Z.erem m (Z.of_int p))) in
  let n = max (Array.length x) (Array.length y) in
  normalize (Array.init n (fun i ->
      let xi = coeff x i and yi = Modp.coeff y i in
      let d = Modp.md p (yi - Z.to_int (Z.erem xi (Z.of_int p))) * minv mod p in
      Z.add xi (Z.mul m (Z.of_int d))))
 
let gcd a b =
  let c = Z.gcd (content a) (content b) in
  let a = primitive a and b = primitive b in
  let g = Z.gcd (lc a) (lc b) in
  (* the gcd, scaled to have leading coefficient g, has coefficients
     below this bound *)
  let bound = Z.mul g (Z.min (mignotte a (min (deg a) (deg b))) (mignotte b (min (deg a) (deg b)))) in
  let primes = Modp.primes_below (1 lsl 30) in
  let rec loop primes m h =
    match primes () with
    | Seq.Nil -> failwith "out of primes"
    | Seq.Cons (p, rest) ->
      if Z.equal (Z.erem g (Z.of_int p)) Z.zero then loop rest m h        (* the degree would drop *)
      else begin
        incr primes_used;
        let hp = Modp.gcd p (Modp.of_z p a) (Modp.of_z p b) in
        let hp = Modp.scale p hp (Z.to_int (Z.erem g (Z.of_int p))) in
        if Modp.deg hp = 0 then [| c |]                                    (* coprime *)
        else if Z.equal m Z.zero || Modp.deg hp < deg h then
          check rest (Z.of_int p) (Array.map Z.of_int hp)                  (* start again *)
        else if Modp.deg hp > deg h then loop rest m h                     (* an unlucky prime *)
        else check rest (Z.mul m (Z.of_int p)) (crt m h p hp)
      end
  and check primes m h =
    if Z.gt m (Z.shift_left bound 1) then
      let cand = primitive (symmetric m h) in
      if divides cand a && divides cand b then scale cand c else loop primes m h
    else loop primes m h
  in
  primes_used := 0;
  loop primes Z.zero [||]

For all but finitely many primes, the gcd modulo $p$ is the image of the true gcd. For the others, the images of $a$ and $b$ have an extra common factor, and the image gcd has too high a degree. A prime is unlucky exactly when it divides a certain nonzero subresultant of $a/h$ and $b/h$; there are only finitely many, and they cannot be recognized in advance^[1][2].

(* Images of the gcd of two coprime polynomials modulo small primes. *)

a = x^2 + 7x + 5, b = x^2 + 3x + 2, resultant 5
  mod  2: gcd = 1
  mod  3: gcd = 1
  mod  5: gcd = x + 2
  mod  7: gcd = 1
  mod 11: gcd = 1
  mod 13: gcd = 1

The only unlucky prime here is 5, which divides the resultant. An image of too high a degree is simply discarded, and one of lower degree than those collected so far shows that all earlier primes were unlucky, so the accumulation starts again from it. Since the images of too high a degree can agree with each other, a run could in principle converge on a wrong candidate; the final trial division rules that out, and the algorithm continues with further primes if it fails.

open Zpoly
 
(* Random inputs a = g f1, b = g f2 with a planted common factor g. *)
let random n = normalize (Array.init (n + 1) (fun i ->
    if i = n then Z.of_int (1 + Random.int 9) else Z.of_int (Random.int 19 - 9)))

deg g  deg a  answer bits  primes  largest subresultant bits
    5     10            3       1                         44
   10     20            4       2                        126
   20     40            4       2                        239
   30     60            4       3                        401
   40     80            4       4                        562

The number of primes depends only on the size of the answer and the bound, and each prime costs a gcd over ${\mathbb{F}}_p$ with word-size arithmetic, quadratic in the degree. The subresultant PRS, on the same inputs, works with coefficients of hundreds of bits at every step. In practice the bound is pessimistic and implementations stop early, as soon as the reconstructed candidate stops changing and divides the inputs^[9].

For polynomials in several variables, Brown's dense algorithm evaluates all but one variable at many points, computes univariate gcds, and interpolates, recursively; it needs a number of evaluations exponential in the number of variables^[1]. Zippel's sparse interpolation of 1979 learns the monomial structure from one image and determines the coefficients from few further evaluations, a probabilistic method suited to sparse inputs^[3]. The EEZ-GCD algorithm of Wang lifts a univariate image to the full gcd by Hensel lifting^[5], and the heuristic GCDHEU of Char, Geddes and Gonnet evaluates the polynomials at a large integer, takes an integer gcd and reads the polynomial gcd off its digits, which is fast when it succeeds^[4].

Brown's 1971 paper introduced the modular method for polynomial gcds, together with the analysis of unlucky primes and evaluation points, and compared it with the remainder-sequence methods of Collins, showing it to be faster by large factors^[1]. Collins had used modular arithmetic for resultants in the same period. Mignotte's bound of 1974 put the termination condition on a firm footing^[6]. Zippel's sparse algorithm^[3], Wang's EEZ-GCD^[5] and GCDHEU^[4] followed, and textbook accounts appeared in Knuth^[7], Geddes, Czapor and Labahn^[8] and von zur Gathen and Gerhard^[2]. Computer algebra systems use modular gcds by default.

• Subresultant PRSThe subresultant polynomial remainder sequence is a way of running Euclid's algorithm on polynomials over the integers, or any integral domain, without fractions and without the exponential coefficient growth of plain pseudo-division. After each pseudo-division, the remainder is divided by a factor computed from the previous leading coefficients and degree differences; the division is always exact, and the remainders obtained are the subresultants of the input, determinants of submatrices of the Sylvester matrix, whose size grows only linearly with the degree. The sequence yields the gcd, up to a constant, and the resultant. It was introduced by Collins in 1967 and developed by Brown and Traub.
• Landau–Mignotte boundThe Landau–Mignotte bound limits the coefficients of any factor of an integer polynomial in terms of the polynomial itself: if g of degree k divides f, then each coefficient of g is at most binom(k, j) times the Euclidean norm of f, up to the ratio of leading coefficients, so all are below 2^k ‖f‖₂. It follows from Landau's inequality for the Mahler measure. Computer algebra algorithms use it to decide how large a modulus is needed to recover a factor or a gcd from its image modulo a prime power or a product of primes, as in Hensel lifting and the modular gcd. Factors can have larger coefficients than the polynomial they divide, so some such bound is necessary.
• Hensel liftingHensel lifting turns a factorization of a polynomial modulo a prime p into a factorization modulo p², p⁴ and so on, as long as the factors modulo p are coprime. Each step solves a linear equation using the Bézout relation between the factors, and in quadratic form it doubles the precision, the p-adic analogue of Newton's method. It is the core of Zassenhaus's algorithm for factoring integer polynomials: factor modulo a small prime, lift until the modulus exceeds twice the Landau–Mignotte bound, and recombine the lifted factors into true factors over the integers.
• Pseudo-divisionPseudo-division is division of polynomials over an integral domain that is not a field, such as the integers or a ring of polynomials in other variables. To avoid fractions, the dividend a is first multiplied by lc(b)^(deg a − deg b + 1), the leading coefficient of the divisor raised to one more than the degree difference; then ordinary long division succeeds and yields a pseudo-quotient q and pseudo-remainder r with lc(b)^e a = q b + r and deg r < deg b. It is the basic step of polynomial remainder sequences over the integers, of resultant and gcd computation, and of Wu's method in geometry theorem proving. Iterated naively, it makes coefficients grow exponentially.
• Fast Fourier transformThe fast Fourier transform (FFT) computes the discrete Fourier transform of n values, the evaluation of a polynomial at the n complex n-th roots of unity, in O(n log n) operations instead of O(n^2). The Cooley-Tukey algorithm splits the input into even and odd positions, transforms each half recursively, and combines them with n/2 butterflies. Because the transform turns convolution into pointwise multiplication, the FFT multiplies polynomials and large integers in near-linear time; over a finite field with suitable roots of unity the same algorithm is exact.
• Knuth's Algorithm DKnuth's Algorithm D divides one multiple-precision integer by another: schoolbook long division in base 2^k, where each quotient digit is guessed from the leading limbs of the remainder and divisor. After normalizing the divisor so its top limb is at least half the base, the guess is never too small and at most two too large, a cheap test with one more limb removes almost all errors, and a rare add-back step fixes the rest. It takes O(mn) operations for an (m+n)-limb dividend and an n-limb divisor.

• Building a CAS in OCaml Part 3Sparse multivariate polynomials, pseudo-division, the subresultant PRS against a modular GCD that is a thousand times faster, resultants, and a parser to drive it all.

1. [1]W. S. Brown, “On Euclid’s algorithm and the computation of polynomial greatest common divisors”, Journal of the ACM 18 (1971).
2. [2]J. von zur Gathen, J. Gerhard, Modern Computer Algebra, ch. 6, Cambridge University Press (3rd ed., 2013).
3. [3]R. Zippel, “Probabilistic algorithms for sparse polynomials”, EUROSAM, Lecture Notes in Computer Science 72 (1979).
4. [4]B. W. Char, K. O. Geddes, G. H. Gonnet, “GCDHEU: heuristic polynomial GCD algorithm based on integer GCD computation”, Journal of Symbolic Computation 7 (1989).
5. [5]P. S. Wang, “The EEZ-GCD algorithm”, ACM SIGSAM Bulletin 14 (1980).
6. [6]M. Mignotte, “An inequality about factors of polynomials”, Mathematics of Computation 28 (1974).
7. [7]D. E. Knuth, The Art of Computer Programming, Vol. 2: Seminumerical Algorithms, §4.6.1, Addison-Wesley (3rd ed., 1997).
8. [8]K. O. Geddes, S. R. Czapor, G. Labahn, Algorithms for Computer Algebra, ch. 7, Kluwer (1992).
9. [9]M. B. Monagan, A. D. Wittkopf, “On the design and implementation of Brown’s algorithm over the integers and number fields”, ISSAC (2000).

last updated 28 September 2026