Pseudo-division
Pseudo-division is the version of polynomial long division used over rings in which not every leading coefficient can be inverted, such as the integers or a polynomial ring [1]. Dividing by over would need the fraction . Pseudo-division avoids it by multiplying the dividend in advance by a power of the divisor's leading coefficient, large enough that every step of the long division is exact. The result is a pseudo-quotient and a pseudo-remainder with coefficients in the ring. Pseudo-remainders are the steps of the polynomial remainder sequences from which gcds and resultants of polynomials over are computed[2][3], and in multivariate form they are the core of Ritt's and Wu's methods for systems of polynomial equations[5][6]. Their drawback is growth: repeated pseudo-division multiplies coefficients so fast that the subresultant and modular algorithms were invented to control it.
Definition
Let be an integral domain and with . Put . There are unique polynomials and in with
called the pseudo-quotient and pseudo-remainder, written and [1]. Existence follows from long division in the fraction field of : each of the steps divides by once, and the factor clears every such denominator. Uniqueness follows from the degree condition, as for ordinary division. Over a field, is times the ordinary remainder, so the two agree up to a unit.
Algorithm
The computation is long division in which, instead of dividing the leading coefficient of the running remainder by , the whole remainder is first multiplied by . With , each of the steps replaces by if , and by otherwise, and updates in the same way.
Pseudo-division over the integers, with Zarith integers as coefficients.
(* Pseudo-division: lc(b)^(deg a - deg b + 1) * a = q * b + r withdeg r < deg b, computed without leaving the integers. *)let pseudo_divmod a b =if is_zero b then invalid_arg "pseudo_divmod";let db = deg b and l = lc b inlet e = deg a - db + 1 inif e <= 0 then ([||], a)else beginlet q = ref [||] and r = ref a infor _ = 1 to e do(* r <- l r - lc(r) x^k b, when deg r >= deg b *)let k = deg !r - db inif k >= 0 then beginlet c = lc !r inq := add (scale !q l) (shift [| c |] k);r := sub (scale !r l) (shift (scale b c) k)end else beginq := scale !q l; r := scale !r lenddone;(!q, !r)endlet prem a b = snd (pseudo_divmod a b)
Two examples, checking the defining identity.
open Zpoly
Running it.
a = x^3 + x + 1, b = 3x^2 + 1q = 3x, r = 6x + 9lc(b)^2 a = q b + r: truea = x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x - 5, b = 3x^6 + 5x^4 - 4x^2 - 9x + 21q = 9x^2 - 6, r = -15x^4 + 3x^2 - 9lc(b)^3 a = q b + r: true
The fixed exponent makes the result unique and easy to reason about, but it can be wasteful: if the remainder's degree drops by more than one in a step, the extra multiplications by are not needed for exactness. The sparse pseudo-remainder omits them and uses the smallest exponent that works[4].
Remainder sequences
Replacing division by pseudo-division in Euclid's algorithm gives a polynomial remainder sequence (PRS) that never leaves . The last nonzero remainder is, up to a factor in , the greatest common divisor, and the last element before zero determines the resultant. The naive sequence, however, is useless in practice.
The example of Knuth: a gcd computation on two polynomials with small coefficients, first with pseudo-remainders alone and then removing the content of each.
open Zpoly(* Knuth's example: the Euclidean algorithm with pseudo-remainders only,against the same with the content removed at each step. *)let rec sequence step a b = if is_zero b then [] else b :: sequence step b (step a b)
Running it.
pseudo-remainders:degree 4, largest coefficient 15degree 2, largest coefficient 59535degree 1, largest coefficient 1654608338437500degree 0, largest coefficient 12593338795500743100931141992187500primitive parts of pseudo-remainders:degree 4, largest coefficient 5degree 2, largest coefficient 49degree 1, largest coefficient 6150degree 0, largest coefficient 1
Each pseudo-remainder is multiplied by a power of the previous leading coefficient, so the number of digits in the coefficients roughly doubles with each step. The polynomials are coprime, and the final result, a constant with 35 digits, only says that the gcd is 1. Dividing each remainder by its content, the gcd of its coefficients, gives the primitive PRS, whose coefficients stay as small as possible, at the price of a gcd of all coefficients at each step[1][2]. The subresultant PRS finds a middle way: it divides by a factor known in advance to be exact, and its coefficients grow only linearly in the degree[2][3].
Multivariate polynomials
A polynomial in several variables can be viewed as a polynomial in one main variable with coefficients in the ring of polynomials in the others, and pseudo-division then eliminates the main variable without introducing rational functions. Ritt used this for systems of differential and algebraic equations, reducing them to triangular characteristic sets[5]. Wu Wen-tsün applied the method to theorems of elementary geometry in 1978: hypotheses and conclusion are written as polynomial equations in the coordinates, the hypotheses are brought into triangular form, and the conclusion follows, under nondegeneracy conditions given by the leading coefficients, if its successive pseudo-remainders with respect to them reduce to zero[6]. Chou proved hundreds of geometry theorems automatically this way[7].
History
Pseudo-division appears in the nineteenth-century theory of resultants and Sturm sequences, where calculations were kept in the integers by clearing denominators. Its systematic study began with computer algebra. Collins introduced reduced and subresultant polynomial remainder sequences in 1967 to control coefficient growth[2], and Brown and Traub developed the theory of subresultants in 1971[3]. Knuth's account in The Art of Computer Programming made the example above the standard illustration of the problem[1]. Ritt's characteristic sets of 1950[5] and Wu's geometry prover of 1978[6] made pseudo-division the central operation of triangular-set methods for polynomial systems.
see also
- 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.
- Modular GCDThe modular GCD algorithm computes the greatest common divisor of two polynomials with integer coefficients by computing it modulo many machine-word primes, where coefficients never grow, and reconstructing the integer coefficients with the Chinese remainder theorem once the product of the primes exceeds twice a bound on them. Primes that divide the leading coefficients are skipped, unlucky primes that give an image of too high a degree are detected and discarded, and a final trial division certifies the result. The intermediate numbers are no larger than the answer, which makes it far faster than remainder-sequence methods; it was introduced by Brown in 1971.
- Gröbner basisA Gröbner basis of an ideal of multivariate polynomials is a generating set with respect to which multivariate division behaves as in one variable: the remainder of any polynomial is unique, and it is zero exactly when the polynomial lies in the ideal. Formally, the leading terms of the basis, for a fixed monomial order, generate the ideal of all leading terms. Buchberger's algorithm computes one from any generating set by adding the remainders of S-polynomials until they all reduce to zero. Gröbner bases decide ideal membership, solve systems of polynomial equations, count their solutions and eliminate variables; Bruno Buchberger introduced them in 1965.
- 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.
- Horner's methodHorner's method evaluates a polynomial a_n x^n + ... + a_0 by nesting the multiplications, (...(a_n x + a_(n-1)) x + ...) x + a_0, which takes n multiplications and n additions, the fewest possible for a general polynomial. The intermediate values are the coefficients of the quotient on division by x - c, so the same loop performs synthetic division, and repeated divisions give the derivatives and the shifted polynomial p(x + c). It is backward stable, but near a multiple root the result can be dominated by rounding errors.
further reading
- [1]D. E. Knuth, The Art of Computer Programming, Vol. 2: Seminumerical Algorithms, §4.6.1, Addison-Wesley (3rd ed., 1997).
- [2]G. E. Collins, “Subresultants and reduced polynomial remainder sequences”, Journal of the ACM 14 (1967).
- [3]W. S. Brown, J. F. Traub, “On Euclid’s algorithm and the theory of subresultants”, Journal of the ACM 18 (1971).
- [4]K. O. Geddes, S. R. Czapor, G. Labahn, Algorithms for Computer Algebra, ch. 7, Kluwer (1992).
- [5]J. F. Ritt, Differential Algebra, American Mathematical Society Colloquium Publications 33 (1950).
- [6]W.-T. Wu, “On the decision problem and the mechanization of theorem-proving in elementary geometry”, Scientia Sinica 21 (1978).
- [7]S.-C. Chou, Mechanical Geometry Theorem Proving, Reidel (1988).
- [8]J. von zur Gathen, J. Gerhard, Modern Computer Algebra, ch. 6, Cambridge University Press (3rd ed., 2013).
- [9]H. Cohen, A Course in Computational Algebraic Number Theory, §3.1, Springer (1993).
last updated