Subresultant PRS

A polynomial remainder sequence (PRS) is the sequence of remainders produced by Euclid's algorithm applied to two polynomials. Over the integers, the remainders are computed by pseudo-division to avoid fractions, and then something must be done about coefficient growth: the plain pseudo-remainders grow exponentially, and dividing each by its content, as the primitive PRS does, requires a gcd of all its coefficients at every step^[4]. The subresultant PRS divides each pseudo-remainder by a factor that can be predicted from the leading coefficients and degrees already seen, without any gcd computation. George Collins proved in 1967 that this division is always exact^[1]. The resulting remainders are, up to sign, the subresultants of the two inputs, determinants built from their coefficients, and Hadamard's inequality then bounds their size linearly in the degree. The sequence gives the gcd of the inputs up to a constant factor and, from its last element, their resultant^[2]^[5].

§ 01

Resultants and subresultants

For $a = \sum_{i = 0}^{m} a_{i} x^{i}$ and $b = \sum_{i = 0}^{n} b_{i} x^{i}$ , the Sylvester matrix is the $(m + n) \times (m + n)$ matrix whose rows are the coefficients of $x^{n - 1} a, \dots, a$ and $x^{m - 1} b, \dots, b$ . Its determinant is the resultant $res (a, b)$ , which vanishes exactly when $a$ and $b$ have a common factor^[6]^[10]. The $j$ -th subresultant $S_{j}$ is the polynomial of degree at most $j$ obtained from the submatrix that uses only $n - j$ shifts of $a$ and $m - j$ shifts of $b$ . Every subresultant is an integer combination $S_{j} = u a + v b$ with $de g u < n - j$ and $de g v < m - j$ , and its coefficients are determinants of integer matrices of size at most $m + n$ .

The fundamental theorem of subresultants says that the remainders of Euclid's algorithm, computed over the fraction field, are proportional to subresultants, and it gives the proportionality constants explicitly in terms of leading coefficients^[1]^[2]. By Hadamard's inequality, the determinant of a $k \times k$ matrix with entries bounded by $B$ is at most $(k B)^{k}$ , so the coefficients of subresultants have $O (k lo g (k B))$ bits: linear in the degree.

§ 02

Algorithm

Let $F_{1} = a$ , $F_{2} = b$ with $de g a \geq de g b$ , and $g_{1} = h_{1} = 1$ . For $i \geq 2$ , with $δ_{i} = de g F_{i - 1} - de g F_{i}$ , the sequence continues as^[3]^[4]:

F_{i + 1} = \frac{prem ( F _{i - 1} , F _{i} )}{g _{i - 1} h _{i - 1}^{δ_{i}}} g_{i} = lc (F_{i}), h_{i} = \frac{g _{i}^{δ_{i}}}{h _{i - 1}^{δ_{i} - 1}}

until the pseudo-remainder is zero. Both divisions are exact; this is the content of Collins's theorem^[1]. The $h_{i}$ are themselves leading coefficients of subresultants. When the sequence ends at a constant, the inputs are coprime and the resultant can be read off the last element with a sign correction; when it ends at a polynomial of positive degree, that polynomial is a multiple of the gcd and its primitive part is the gcd.

The subresultant PRS and the resultant computed from it, and a determinant of the Sylvester matrix by fraction-free elimination for comparison.

(* The subresultant polynomial remainder sequence: each pseudo-remainder
   is divided by g h^delta, which is exact, and the remainders are then
   the subresultants of the inputs (up to sign). *)
let subresultant_prs a b =
  let a, b = if deg a >= deg b then (a, b) else (b, a) in
  let rec go u v g h acc =
    if is_zero v then List.rev acc
    else
      let delta = deg u - deg v in
      let r = prem u v in
      if is_zero r then List.rev (v :: acc)
      else
        let r = divexact_const r (Z.mul g (Z.pow h delta)) in
        let g' = lc v in
        let h' = if delta = 0 then h else Z.divexact (Z.pow g' delta) (Z.pow h (delta - 1)) in
        go v r g' h' (v :: acc)
  in
  go a b Z.one Z.one [ a ]
 
(* The resultant, following the same sequence (Cohen, Algorithm 3.3.7). *)
let resultant a b =
  if is_zero a || is_zero b then Z.zero
  else begin
    let signed_content p = let c = content p in if Z.sign (lc p) < 0 then Z.neg c else c in
    let t = Z.mul (Z.pow (signed_content a) (deg b)) (Z.pow (signed_content b) (deg a)) in
    let a = primitive a and b = primitive b in
    let s = ref (if deg a < deg b && deg a mod 2 = 1 && deg b mod 2 = 1 then -1 else 1) in
    let a, b = if deg a >= deg b then (a, b) else (b, a) in
    let rec go a b g h =
      if is_zero b then Z.zero
      else if deg b = 0 then
        if deg a = 0 then Z.mul (Z.of_int !s) t
        else Z.mul (Z.of_int !s) (Z.mul t (Z.divexact (Z.pow (lc b) (deg a)) (Z.pow h (deg a - 1))))
      else begin
        let delta = deg a - deg b in
        if deg a mod 2 = 1 && deg b mod 2 = 1 then s := - !s;
        let r = prem a b in
        if is_zero r then Z.zero
        else
          let r = divexact_const r (Z.mul g (Z.pow h delta)) in
          let g' = lc b in
          let h' = if delta = 0 then h else Z.divexact (Z.pow g' delta) (Z.pow h (delta - 1)) in
          go b r g' h'
      end
    in
    go a b Z.one Z.one
  end

Knuth's example.

open Zpoly

Running it.

  x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x - 5
  3x^6 + 5x^4 - 4x^2 - 9x + 21
  -15x^4 + 3x^2 - 9
  65x^2 + 125x - 245
  -9326x + 12300
  260708
resultant from the sequence: 260708
Sylvester determinant:       260708

The coefficients 15, 65, 9326 and the resultant 260708 are the subresultants of the two polynomials; they match the determinant of the $14 \times 14$ Sylvester matrix computed independently.

§ 03

Coefficient growth

Bits in the largest coefficient of each remainder, for a random pair of degree 20 and 19 with one-digit coefficients.

open Zpoly
 
(* Remainder sequences of two random polynomials of degree 20 and 19 with
   coefficients between -9 and 9: bits in the largest coefficient of each
   remainder, for three ways of keeping the sequence in the integers. *)
let rec seq step a b = if is_zero b then [] else b :: seq step b (step a b)

Running it; the pseudo-remainder sequence was stopped after nine steps.

degree  pseudo-remainders  primitive PRS  subresultant PRS
    19                  4              4                 4
    18                  8              8                 8
    17                 16             16                16
    16                 36             22                22
    15                 84             29                29
    14                199             36                38
    13                478             44                45
    12               1152             53                53
    11               2779             59                59
    10                                64                69
     9                                75                78
     8                                86                87
     7                                97                97
     6                               107               107
     5                               116               116
     4                               125               125
     3                               133               133
     2                               136               142
     1                               147               149
     0                                 1               156

The same data on a logarithmic scale. Plain pseudo-remainders grow exponentially and pass 2,000 bits after nine steps. The primitive and subresultant sequences grow linearly and stay within about 150 bits; the primitive sequence ends at 1 because the inputs are coprime, the subresultant sequence at the resultant.

The subresultants are sometimes slightly larger than the primitive remainders, since they are not reduced by their content, but they are obtained without a single gcd of coefficients, which on large inputs is the more expensive operation^[3]. For computing gcds of polynomials over the integers, both are now superseded by the modular algorithm, whose intermediate numbers are as small as the answer; subresultants remain the method of choice when the resultant itself or the subresultants are wanted.

§ 04

Applications

Resultants eliminate variables: the resultant of two bivariate polynomials with respect to $y$ is a polynomial in $x$ whose roots are the $x$ -coordinates of their common zeros. Discriminants, which detect repeated roots, are resultants of a polynomial and its derivative. Collins developed the subresultant theory in part for his cylindrical algebraic decomposition, which decides the first-order theory of the real numbers and computes subresultant coefficients of many pairs of polynomials^[9]. Sturm–Habicht sequences, built from subresultants, count real roots without fractions^[7]^[11].

§ 05

History

Sylvester introduced the matrix of the resultant in 1840 and subresultants in his memoir of 1853^[6]. Habicht used them in 1948 for a fraction-free version of Sturm's theorem^[7]. With the arrival of computer algebra, Collins proved the fundamental theorem and gave the reduced and subresultant PRS algorithms in 1967^[1], and Brown and Traub reformulated the theory in 1971^[2]. Brown's 1978 paper gave the algorithm in its standard form^[3], and Knuth's presentation as Algorithm C made it widely known^[4]. Ducos gave an improved version in 2000 that reduces the cost of the exact divisions^[8].