Subresultant PRS

also: polynomial remainder sequence, prs, algorithm c

The polynomial remainder sequence that keeps coefficient growth linear rather than doubly exponential. Each remainder is divisible by a factor predictable from the previous two degrees and leading coefficients, and dividing it out at every step holds the intermediates near the true subresultants.

The growth in a naive pseudo-remainder sequence is an artefact. Each remainder is divisible by a factor that can be predicted from the degrees and leading coefficients of the previous two, and dividing it out at every step keeps the entries at the size of the true subresultants, which grow linearly in the degree rather than doubly exponentially.

Knuth's Algorithm C. The division by g*h^delta is exact, which is the content of the theorem.

let g = ref Bigint.one and h = ref Bigint.one in
while not !finished do
  let delta = degree !u - degree !v in
  let r = prem !u !v in
  if is_zero r then (result := !v; finished := true)
  else if degree r = 0 then (result := one; finished := true)
  else begin
    u := !v;
    v := divexact_int r (Bigint.mul !g (Bigint.pow !h delta));
    g := lc !u;
    h := if delta = 0 then !h
         else Bigint.div (Bigint.pow !g delta) (Bigint.pow !h (delta - 1))
  end
done

The bookkeeping in g and h is the whole algorithm: get it wrong and the division is no longer exact, which shows up immediately as a non-integer rather than as a slow path.

It is correct and it is still not fast. The subresultants themselves are large, so a degree-50 input has coefficients hundreds of digits wide and every operation is on bignums. That is what the modular approach avoids.