Algorithms for Polynomial GCD Computation over Algebraic Function Fields

Mark van Hoeij^{∗}

Department of Mathematics

Florida State University Tallahassee, FL 32306-4510, USA.

Michael Monagan^{†}

Department of Mathematics Simon Fraser University Burnaby, B.C. Canada. V5A 1S6.

ABSTRACT

Let L be an algebraic function field in k ≥ 0 parameters t_{1}, …, t_{k}. Let f_{1}, f_{2} be non-zero polynomials in L[x]. We give two algorithms for computing their gcd. The first, a modular GCD algorithm, is an extension of the modular GCD algorithm of Brown for ℤ[x_{1}, …, x_{n}] and Encarnacion for ℚ(α)[x] to function fields. The second, a fractionfree algorithm, is a modification of the Moreno Maza and Rioboo algorithm for computing gcds over triangular sets. The modification reduces coefficient growth in L to be linear. We give an empirical comparison of the two algorithms using implementations in Maple.

1. INTRODUCTION

Let D = ℤ[t_{1}, …, t_{k}] and F = ℚ(t_{1}, …, t_{k}). Let m(z) ∈ F [z] be monic and irreducible of degree d and let L = F [z]/⟨m(z)⟩. L is an algebraic function field of degree d in k parameters t_{1}, …, t_{k}. In our examples, if k = 1 we use t without subscript to denote the parameter. Our problem is to compute the gcd of two non-zero polynomials in L[x]. We denote the input polynomials by f_{1} and f_{2} and their monic gcd by g.

Our first algorithm, presented in section 2, is a modular GCD algorithm. It uses rational number and rational function reconstruction to recover the coefficients of the gcd and trial division to prove the correctness of the result. Like Encarnacion’s algorithm [3], our algorithm is output sensitive, that is, the number of primes and evaluation points it uses depends on the size of the gcd and not on bounds based on the size of the inputs. As in [5], our algorithm does not compute discriminants. We show also how to use our algorithm to compute gcd’s in L[x_{1}, …, x_{n}] by moving the variables x_{2}, …, x_{n} into F.

∗Supported by NSF grant 0098034. †Supported by NSERC of Canada and the MITACS NCE of Canada.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for pro�t or commercial advantage and that copies bear this notice and the full citation on the �rst page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior speci�c permission and/or a fee. ISSAC’04, July 4�7, 2004, Santander, Spain. Copyright 2004 ACM 1-58113-827-X/04/0007 ...$5.00.

We may also compute g using the Euclidean algorithm. But, there is a linear growth in the degrees and heights of the coefficients in F in the Euclidean algorithm which causes it to bog down doing arithmetic in F. Arithmetic with fractions in F and hence gcd computation in D can be eliminated using a fraction-free approach. In [6] Moreno Maza and Rioboo show how to do this for univariate gcd computation modulo a triangular set for which our setting, L[x], is a special case. In section 3 we demonstrate that their algorithm has a serious coefficient swell in D. We modify their algorithm to have linear growth by doing O(δ) gcd computations in D where δ = max(deg_{x} f_{1}, deg_{x} f_{2}) − deg_{x} g + 1. In section 4 we compare Maple implementations for the two algorithms to demonstrate their effectiveness.

Both of our algorithms work with associates (scalar multiples) of f_{1}, f_{2} and g. We make some definitions. A non-zero polynomial in D[z, x] is said to be primitive wrt (z, x) if the gcd of its coefficients in D is 1. Let f be nonzero in L[x]. The denominator of f is the polynomial den(f) ∈ D of least (total) degree in (t_{1}, …, t_{k}) and with smallest integer content such that den(f)f is in D[z, x]. The primitive associate f̌ of f is the associate of den(f)f which is primitive in D[z, x]. These definitions for den(f) and f̌ are unique up to sign; we impose uniqueness by requiring them to have positive leading coefficient in a term ordering. The monic associate f̃ of f is defined as f̃ = ȟ where h = monic(f). Here monic(f) is defined as lc_{x}(f)^{−1}f where lc_{x}(f) ∈ L is the leading coefficient of f wrt x. Notice that deg_{z}(lc_{x} f̃) = 0.

Examples: Let f = 2tx − t/(1 − t) ∈ ℚ(t)[x]. Then den(f) = t − 1, monic(f) = x − 1/(2 − 2t) and f̌ = f̃ = (2t − 2)x + 1. Let f = 6zx − 3t ∈ ℚ(t)[z][x] where z = \sqrt{t}. Then f̌ = 2zx − t, monic(f) = x − z/2 and f̃ = 2x − z.

To provide the reader with an overview of the two algorithms, we first work through an example. Let z = \sqrt{t}. Consider the input polynomials

f_{1} = x^{2} + \frac{−2 t + 3}{3} zx + \frac{5}{t} x + \frac{5}{t} z + \frac{−2t^{2}}{3} ,
f_{2} = zx^{2} + \frac{5}{t} zx + \frac{−3 + 2t^{2}}{3} x + \frac{−2t}{3} z + \frac{5}{t} .

Since f_{1} = (x + z)(x − 2tz/3 + 5/t) and f_{2} = (zx + 1)(x − 2tz/3 + 5/t) the monic gcd g of f_{1} and f_{2} is the polynomial x−2tz/3+5/t. Both algorithms output ǧ = 3tx−2t^{2}z +15.

On input of f_{1} and f_{2} they first compute

f̌_{1} = 3tx^{2} + (−2t^{2} + 3t)zx + 15x + 15z − 2t^{3},
f̌_{2} = 3tzx^{2} + 15zx + (−2t^{3} + 3t)x − 2t^{2}z + 15.

The Fraction-Free Algorithm

Let r_{1}, r_{2}, r_{3}, …, r_{n}, r_{n+1} = 0 be the remainder sequence for the Euclidean algorithm with input r_{1} := f̌_{1}, r_{2} := f̌_{2}. Our fraction-free algorithm computes the remainder sequence ř_{1}, r̃_{2}, r̃_{3}, …, r̃_{n} = ǧ, r̃_{n+1} = 0 without introducing fractions in ℚ(t). We multiply f̌_{2} by z, a quasi-inverse (see section 3) of lc_{x}f̌_{2}, to eliminate z from the leading coefficient of f̌_{2}.

p_{2} := zf̌_{2} = 3t^{2}x^{2} + (−2t^{3} + 3t)zx + 15tx + 15z − 2t^{3}.

Since p_{2} is primitive, p_{2} = f̃_{2} = r̃_{2}. Using pseudo-division to avoid fractions in ℚ(t), we compute the remainder of μf̌_{1} divided by p_{2} in x where μ = lc_{x}p_{2} = 3t^{2} is an element of D, all working modulo m(z) = z^{2} − t. We obtain

p_{3} = 3t(t − 1)zx + 15(t − 1)z − 2t^{3}(t − 1).

To minimize coefficient growth in ℤ[t] we make p_{3} primitive wrt (z, x); we compute and divide out by gcd(3t(t−1), 15(t− 1), −2t^{3}(t − 1)) = t − 1. We obtain

p_{3} := p_{3}/(t − 1) = 3tzx + 15z − 2t^{3}.

Now we compute p_{3} := r̃_{3} by first multiplying p_{3} by z and then making the result primitive. We obtain

p_{3} := zp_{3}/t = 3tx − 2t^{2}z + 15.

Now we divide p_{2} by p_{3} using pseudo-division modulo m(z). The pseudo-remainder is 0, thus, ǧ = p_{3} and we are done.

The Modular GCD Algorithm

We compute the gcd of f̌_{1} and f̌_{2} modulo a sequence of primes. Unlike the fraction-free algorithm, we do not compute f̃_{2}, because inverting lc_{x}(f_{2}) may lead to a blowup. Suppose we start with the prime p = 11. We will obtain

g_{p} = tx + 3t^{2}z + 5 mod 11.

We apply Wang’s rational reconstruction (see [7, 3]) to the coefficients of g_{p} modulo p. It fails so we compute the gcd modulo a second prime, q = 13. We will obtain

g_{q} = tx − 5t^{2}z + 5 mod 13.

Note that we normalized the leading coefficient (of tx) to be 1 for both images. We apply the Chinese remainder theorem to obtain g_{m} mod m = pq such that g_{m} ≡ g_{p} mod p and g_{m} ≡ g_{q} mod q.

g_{m} = tx + 47t^{2}z + 5 mod 143.

We apply rational reconstruction to the coefficients of g_{m} modulo m = 143. This time we succeed. We obtain

h = tx − \frac{2}{3} t^{2}z + 5.

Now we clear denominators; we set

h := 3h = 3tx − 2t^{2}z + 15

and test if h|f̌_{1} and h|f̌_{2}. It does, therefore h = ǧ and we are done. It remains to describe how we compute g_{p} and g_{q}. We do so by computing the gcd of f̌_{1}(t, x) and f̌_{2}(t, x) modulo p at a sequence of values for t in ℤ_{p}. The evaluation point t = 0 causes lc_{x}(f̌_{2}) to vanish, so it is lc-bad (defined in section 2.1). Our algorithm avoids such t. Now consider t = 1 (our algorithm chooses evaluation points at random, but for this example we will use t = 1, 2, …). We run the Euclidean algorithm modulo p = 11 on

f_{1}(1, x) = 3x^{2} + zx + 4x + 4z − 2, and
f_{2}(1, x) = 3zx^{2} + 4zx + x − 2z + 4.

modulo m(1, z) = z^{2} − 1. This is a univariate computation over ℤ_{p}[z]/⟨m(1, z)⟩, a finite ring, so no coefficient growth can occur. It succeeds because it did not try to divide by zero divisor in ℤ_{p}[z]/⟨z^{2} − 1⟩. The output is

g_{1} = gcd(f_{1}(1, x), f_{2}(1, x)) = x^{2} + 4zx + 5x + 5z + 3

Now we apply rational function reconstruction to the coefficients w.r.t. (z, x). It succeeds with output h = g_{1}, i.e., with constant rational functions in t. We test if h|f_{1}(t, x) mod 11 and h|f_{2}(t, x) mod 11. It doesn’t so we take a new evaluation point t = 2 and run the Euclidean algorithm over ℤ_{p}[z]/⟨m(2, z)⟩. We obtain

g_{2} = gcd(f_{1}(2, x), f_{2}(2, x)) mod 11 = x − 5z − 3

which has lower degree than g_{1}. This tells us (see section 2.1) that t = 1 is unlucky and thus g_{1} cannot be an associate of g(1, x) mod 11. We discard g_{1}. Continuing, we apply rational function reconstruction to g_{2} and obtain h = g_{2}. We test if h|f_{1}(t, x) mod 11 and h|f_{2}(t, x) mod 11. It doesn’t so we take a new evaluation point t = 3 and compute

g_{3} = gcd(f_{1}(3, x), f_{2}(3, x)) = x − 2z − 2.

We apply the Chinese remainder theorem to obtain the gcd modulo (t − 2)(t − 3). We obtain

c = x + (3t)z + (t − 5).

We apply rational function reconstruction to the coefficients of c to reconstruct numerators of degree 1 in t and denominators of degree 0 in t. This succeeds with h = c. We test if h|f_{1}(t, x) mod 11 and h|f_{2}(t, x) mod 11. It doesn’t so we take a new evaluation point t = 4 and compute

g_{4} = gcd(f_{1}(4, x), f_{2}(4, x)) = x + z + 4.

Apply the Chinese remainder theorem to c and g_{4} to obtain

c = x + 3tz − 3t^{2} + 5t − 1

the gcd modulo (t − 2)(t − 3)(t − 4). We apply rational function reconstruction to the coefficients of c to reconstruct numerators of degree 1 and denominators of degree 1 in t. It succeeds with output

h = x + \frac{3t}{1} z + \frac{5}{t} .

Now we clear the denominators in t. We obtain

h := th = tx + 3t^{2}z + 5

Since h|f_{1}(t, x) mod 11 and h|f_{2}(t, x) mod 11 we are done with p = 11. In this example all prime numbers used were good. However, we did encounter an lc-bad evaluation t = 0 and an unlucky evaluation t = 1. It is also possible to hit a zero divisor in the middle of the Euclidean algorithm but such a fail prime/evaluation did not occur in the example. In section 2.1 we study the different problem cases that can occur before describing the algorithm in section 2.2 which we can prove (using Theorem 1 and 2) is correct.

2. THE MODULAR GCD ALGORITHM

2.1 lc-bad, fail, unlucky, and good primes

Let f_{1}, f_{2} ∈ L[x] where m ∈ F [z], L = F [z]/⟨m⟩ and F = ℚ(t_{1}, …, t_{k}). Let g ∈ L[x] be their monic gcd. Our algorithm in section 2.2 replaces f_{1}, f_{2}, m by their primitive associates, so we can assume that m = m̌ ∈ D[z] and f_{i} = f̌_{i} ∈ D[z, x] (mod m) where D = ℤ[t_{1}, …, t_{k}]. Let l_{m} = lc_{z}(m) ∈ D, and l_{i} = lc_{x}(f_{i}) ∈ D[z]. If I is a prime ideal in D, and if l_{m} ∈ I or l_{1} ∈ ⟨I, m⟩ or l_{2} ∈ ⟨I, m⟩ then I is called an lc-bad prime. Such I will be avoided.

Let I be a maximal ideal in D. Denote \overline{D} = D/I and \overline{R} = D[z]/⟨I, m⟩. For s ∈ D[z, x] define \overline{s} := s mod ⟨I, m⟩. Thus \overline{s} ∈ \overline{R}[x]. Note that \overline{l}_{m}, \overline{l}_{1} and \overline{l}_{2} are not zero because we avoid lc-bad I. The maximal ideals I used in our algorithm are of the form I = ⟨p, t_{1} − α_{1}, …, t_{k} − α_{k}⟩, with p a prime number and α_{i} ∈ ℤ. For such I we can identify \overline{D} resp. \overline{R} with ℤ_{p} resp. ℤ_{p}[z]/⟨m(α_{1}, …, α_{k}, z)⟩.

Definition: A monic gcd of \overline{f}_{1}, \overline{f}_{2} is a monic polynomial g_{I} ∈ \overline{R}[x] for which ⟨\overline{f}_{1}, \overline{f}_{2}⟩ = ⟨g_{I}⟩ as \overline{R}[x]-ideals.

In general, a monic gcd need not exist, but if it does then it is unique. The Euclidean algorithm applied to \overline{f}_{1}, \overline{f}_{2} ∈ \overline{R}[x] fails if it tries to invert a leading coefficient in the polynomial remainder sequence that happens to be a zero-divisor (for more details see section 2 in [5]). In this case I is called a fail prime. This can happen even if the monic gcd of \overline{f_{1}}, \overline{f_{2}} exists. As an example, take f_{1} = f_{2} + 1, f_{2} = (z + 1)x + t, m = z^{2} + 7tz − 1. Then ⟨7, t − α⟩ resp. ⟨p, t − 0⟩ is a fail prime for any integer α resp. prime number p.

We may assume (if not, interchange f_{1}, f_{2}) that the Euclidean algorithm applied to \overline{f}_{1}, \overline{f}_{2} first inverts the leading coefficient of \overline{f}_{2} in order to divide \overline{f}_{1} by \overline{f}_{2}. So if I is not lc-bad nor fail then \overline{l}_{2} is a unit in \overline{R}, which is important in the proof of the theorem below. Note that if the implementation (without first interchanging) always inverts \overline{l}_{2} like in [5] then we may use the asymmetric definition of lc-bad: I is lc-bad when \overline{l}_{m} = 0 or \overline{l}_{2} = 0.

Assume now that the Euclidean algorithm applied to \overline{f_{1}}, \overline{f_{2}} does not fail. Then the monic gcd g_{I} ∈ \overline{R}[x] exists and will be the output. Let d_{I} := deg_{x}g_{I} and d_{g} := deg_{x}g. If d_{I} > d_{g} then I is called an unlucky prime.

If den(g) ∉ I then \overline{g} ∈ \overline{R}[x] is defined, \overline{g} := \overline{den(g)g}/\overline{den(g)}. If this \overline{g} equals g_{I} then I is called a good prime.

Theorem 1. With the above notations, if I is not lc-bad nor fail, then it is either unlucky or good.

We will prove a more general statement:

Theorem 2. Let D be an integral domain, F its field of fractions, let m ∈ D[z] (not necessarily irreducible in F [z]). Let L = F [z]/⟨m⟩. Let f_{1}, f_{2} ∈ D[z, x] represent two nonzero polynomials f_{1}, f_{2} ∈ L[x], and let g ∈ L[x] be any monic common divisor of f_{1}, f_{2} in L[x]. Let I be a maximal ideal in D that is not lc-bad, let \overline{D} = D/I and \overline{R} = D[z]/⟨I, m⟩, and assume that the result g_{I} of the Euclidean algorithm applied to \overline{f}_{1}, \overline{f}_{2} ∈ \overline{R}[x] is not “failed”. Then d_{I} ≥ d_{g}. Furthermore, if d_{I} = d_{g} then \overline{g} is defined and \overline{g} = g_{I}, and the monic gcd of f_{1}, f_{2} exists and equals g.

Proof: If the Euclidean algorithm does not fail then the extended Euclidean algorithm will not fail either, and will produce polynomials s_{I}, t_{I} ∈ \overline{R}[x] such that s_{I} \overline{f}_{1}+t_{I} \overline{f}_{2} = g_{I} and

deg_{x}s_{I} < d_{2} − d_{I}, deg_{x}t_{I} < d_{1} − d_{I}.

Here d_{1} = deg_{x}f_{1}, d_{2} = deg_{x}f_{2}, d_{I} = deg_{x}g_{I}, and d_{g} = deg_{x}g. Note that the degree-bound given for t_{I} is not optimal because deg_{x}\overline{f}_{1} could be smaller than d_{1} if we use the asymmetric definition of lc-bad.

Let D_{I} = {c ∈ F | den(c) ∉ I} = {a/b | a, b ∈ D, b ∉ I}. Then we can identify D_{I}/⟨I⟩ with \overline{D}. An element c ∈ D_{I} is a unit in D_{I} if and only if its image \overline{c} in \overline{D} is not zero. Denote R_{I} = {a/b ∈ L | a ∈ D[z] mod m, b ∈ D, b ∉ I} which we can identify with D_{I}[z]/⟨m⟩ (note that l_{m} is a unit in D_{I} because I is not lc-bad, hence division with remainder by m works as usual in D_{I}[z]). Then we can identify R_{I}/⟨I⟩ = D_{I}[z]/⟨I, m⟩ = \overline{R}.

Consider the free R_{I}-modules: H = {h ∈ R_{I}[x] | deg_{x}h < d_{1} + d_{2} − d_{I}}, H_{0} = {h ∈ R_{I}[x] | deg_{x}h < d_{I}}, G = H/H_{0} and S = {(s, t) | s, t ∈ R_{I}[x], deg_{x}s < d_{2} − d_{I}, deg_{x}t < d_{1} − d_{I}}. Now define a map ϕ : S → G given by ϕ(s, t) = sf_{1} + tf_{2} + H_{0}. Likewise (just replace R_{I} by \overline{R}) define free \overline{R}-modules \overline{H}, \overline{H}_{0}, \overline{S} and \overline{G}, as well as the map \overline{ϕ} : \overline{S} → \overline{G} given by \overline{ϕ}(s, t) = s\overline{f}_{1} + t\overline{f}_{2} + \overline{H}_{0}. A basis for G resp. S as R_{I}-modules is B_{G} := {x^{i} | d_{I} ≤ i < d_{1} + d_{2} − d_{I}} resp. B_{S} := {(x^{i}, 0), (0, x^{j}) | i < d_{2} − d_{I}, j < d_{1} − d_{I}}. The same B_{G} resp. B_{S} is also a basis for \overline{G} resp. \overline{S}.

Counting B_{S} and B_{G} one sees that S and G are free R_{I}-modules of rank n_{0} := d_{1} + d_{2} − 2d_{I}. Viewed as free D_{I}-modules, their ranks are n := deg_{z}(m)n_{0}, so ϕ can be represented by an n by n matrix M with entries in D_{I}. In the same way, we can view \overline{S} and \overline{G} as \overline{D}-vector spaces of dimension n, and \overline{ϕ} is given by an n by n matrix \overline{M} which equals M mod ⟨I⟩.

We will show that \overline{ϕ} is one to one. The extended Euclidean algorithm produces (s_{I}, t_{I}) ∈ \overline{S} with s_{I} \overline{f}_{1} + t_{I} \overline{f}_{2} = g_{I}. Now suppose that \overline{ϕ}(s, t) = 0, i.e. s\overline{f}_{1} + t\overline{f}_{2} ∈ \overline{H}_{0}. Now \overline{f}_{1}, \overline{f}_{2}, and hence s\overline{f}_{1} +t\overline{f}_{2} are divisible by g_{I} (divisions by g_{I} work because g_{I} is monic). But g_{I} has higher degree than any element of \overline{H}_{0}. Hence s\overline{f}_{1} + t\overline{f}_{2} = 0. Then s\overline{f}_{1} and thus also ss_{I} \overline{f}_{1} are divisible by \overline{f}_{2} (divisions by \overline{f}_{2} work because \overline{l}_{2} is a unit). Then 0 ≡ ss_{I} \overline{f}_{1} ≡ s(s_{I} \overline{f}_{1} + t_{I} \overline{f}_{2}) ≡ sg_{I} mod \overline{f}_{2}. So sg_{I} vanishes mod \overline{f}_{2}, yet it has lower degree than \overline{f}_{2}, so sg_{I} must be 0, which implies s = 0 because g_{I} is monic. Now t\overline{f}_{2} = s\overline{f}_{1} + t\overline{f}_{2} = 0, which implies t = 0 because \overline{l}_{2} is a unit. This proves that \overline{ϕ} is one to one, so det(\overline{M}) ∈ \overline{D} is not zero. So the image of det(M) in \overline{D} is not zero. Then det(M) is a unit in D_{I}. This implies that ϕ is invertible.

Now take w = x^{d_{I}} + H_{0} ∈ G and \overline{w} = x^{d_{I}} + \overline{H}_{0} ∈ \overline{G}. Note that \overline{w} = g_{I} + \overline{H}_{0} and hence \overline{ϕ}^{−1}(\overline{w}) = (s_{I}, t_{I}). Now consider ϕ^{−1}(w), which is in S, and hence can be written as (s, t) where s, t ∈ R_{I}[x]. By definition of ϕ, we see that sf_{1} + tf_{2} + H_{0} = x^{d_{I}} + H_{0} so h := sf_{1} + tf_{2} ∈ R_{I}[x] is a monic polynomial of degree d_{I}.

Since \overline{M} is M mod ⟨I⟩ we see that s_{I} = \overline{s}, t_{I} = \overline{t} and g_{I} = \overline{h}. Here \overline{s}, \overline{t}, \overline{h} are the images of s, t, h in \overline{R}[x]. Now g ∈ L[x] was assumed to be a monic common factor of f_{1}, f_{2}, so it must be a monic factor in L[x] of h as well. This shows d_{g} ≤ deg_{x}h = d_{I}. If d_{g} = d_{I} then the monic polynomial h is divisible by the monic polynomial g of the same degree, and hence g = h. Then g ∈ R_{I}[x] so \overline{g} is defined and equals \overline{h} = g_{I}. Since g is assumed to be a monic common factor of f_{1}, f_{2}, and g = sf_{1} + tf_{2} we see that g is the monic gcd of f_{1}, f_{2}. This completes the proof.

2.2 Algorithm MODGCD

Recall that D = ℤ[t_{1}, …, t_{k}], F = ℚ(t_{1}, …, t_{k}), and L = F [z]/⟨m(z)⟩. Our modular GCD algorithm calls subroutine M which calls recursive subroutine P which in turn calls the Euclidean algorithm over a finite ring.

Algorithm MODGCD Input: Non-zero polynomials f_{1}, f_{2} ∈ L[x] and m(z). Output: ǧ, where g is the monic gcd of f_{1} and f_{2} in L[x].

Call subroutine M with input f̌_{1}, f̌_{2} and m̌(z).

Subroutine M. Input: f_{1}, f_{2} ∈ D[z]/⟨m(z)⟩[x] and m(z) ∈ D[z]. Output: ǧ, where g is the monic gcd of f_{1} and f_{2}.

1. Set n = 1.

2. Main Loop: Take a new prime p_{n} that is not lc-bad.

3. Let g_{n} ∈ ℤ_{p_{n}}[t_{1}, …, t_{k}][z, x] be the output of subroutine P applied to f_{1}, f_{2}, m(z) mod p_{n}.

4. If g_{n} =“failed” then go back to step 2.

5. If g_{n} = 1 then return 1.

6. Select from {g_{1}, … g_{n}} those with the same leading term (in pure lexicographic order with x > t_{1} > ⋯ > t_{k}) as g_{n}. Combine those using Chinese remaindering to obtain c mod m_{c}. Set n = n + 1.

7. Apply integer rational reconstruction to obtain h from c mod m_{c}. If this fails, go back to step 2.

8. Clear fractions in ℚ : Set h = ȟ.

9. Trial division: If h|f_{1} and h|f_{2} then return h, otherwise, go back to step 2.

Subroutine P. Input: f_{1}, f_{2} ∈ D_{p}[z]/⟨m(z)⟩[x] and m(z) ∈ D_{p}[z] where D_{p} = ℤ_{p}[t_{1}, …, t_{k}]. Output: Either ǧ, the primitive associate of the monic gcd of f_{1} and f_{2}, or “failed.”

0. If k = 0 then output the result of the Euclidean algorithm applied to f_{1}, f_{2} (which returns “failed” if it tries to divide by a zero-divisor).

1. Set n = 1 and d = 1.

2. Main Loop: Take a random new evaluation point α_{n} in ℤ_{p} that is not lc-bad.

3. Let g_{n} ∈ ℤ_{p}[t_{1}, …, t_{k−1}][z, x] be the output of subroutine P applied to f_{1}, f_{2}, m(z) at t_{k} = α_{n}.

4. If g_{n} =“failed” then Set d = d + 1. If d > n output “failed”, else go back to step 2.

5. If g_{n} = 1 then return 1.

6. Select from {g_{1}, … g_{n}} those with the same leading term in x, t_{1}, …, t_{k−1} as g_{n}. Chinese remainder those to obtain c mod m_{c}(t_{k}). Set n = n + 1.

7. Apply rational function reconstruction to the coefficients of c to obtain h ∈ ℤ_{p}(t_{k})[t_{1}, …, t_{k−1}][z, x] s.t. h ≡ c mod m_{c}(t_{k}). If this fails, go back to step 2.

8. Clear fractions in ℤ_{p}(t_{k}): Set h = ȟ.

9. Trial division: If h|f_{1} and h|f_{2} then return h, otherwise, go back to step 2.

If we often encounter “failed”, or if we run out of evaluation points in step 2 in subroutine P, then we should increase the size of the primes used (this will not be necessary in practice if one uses 30 bit primes). Since lc-bad and fail prime/evaluations are discarded in steps 2 and 4, of subroutines M and P, the prime/evaluations remaining are either unlucky or good by Theorem 2. The g_{i} combined in step 6 have the same leading term, thus the same degree in x, and hence they are either all unlucky or all good.

The images g_{i} in step 6 must be scaled in a consistent way before applying Chinese remaindering, see also Remark 2 below. The rational reconstruction in step 7 in subroutine M is accomplished with Wang’s algorithm [7, 2, 4]. Suppose this succeeds and step 7 finds h. If h|f_{1} and h|f_{2} then ȟ = ǧ and M terminates. If either trial division fails then Theorem 1 implies either m_{c} is not yet large enough to recover the coefficients in ℚ or all primes dividing m_{c} were unlucky.

Zero divisors in subroutine P

Subroutine P attempts to compute (if it exists) the gcd of f_{1}, f_{2} in characteristic p by choosing points α = (α_{1}, …, α_{k}) in ℤ_{p}^{k} and computing gcd(f_{1}(α), f_{2}(α)) with the Euclidean algorithm (in step 0). This is a gcd computation over a finite ring which may fail if a zero divisor is encountered. If we use sufficiently large prime numbers and random evaluation points then the probability of such failure is small, so our modular GCD algorithm will terminate.

We explain the role of the variable d in subroutine P. Suppose m(z) = z^{2} + 7t − 1 and f_{1} = x^{2} + t and f_{2} = (z + 1)x + t. If subroutine M chose the prime p = 7 then when the Euclidean algorithm is called by algorithm P (in step 0) with inputs f_{1}(α) and f_{2}(α), the first division would fail when it attempts to invert lc_{x}(f_{2}) = z + 1 which is a zero-divisor for any choice of α ∈ ℤ_{7}. So when p = 7 then step 0 fails with probability 1. Suppose subroutine M chose a prime p ≠ 7. This time only one evaluation point (α = 0 for t) results in a zero divisor being encountered, so step 0 fails with probability 1/p. When most evaluations fail (in the example when p = 7) then d > n will soon hold in step 4, so subroutine P wastes little time on p and gives up quickly. If most evaluations are good, and if subroutine P has already computed many good g_{n}, then the test d > n prevents, with high probability, that a few unlucky choices in step 2 could cause a lot of useful work to be lost.

Remark 1: It is essential that the evaluation point α_{n} ∈ ℤ_{p} be chosen at random in step 2. Consider m(z) = z^{2} + (t − 1)z −1 and f_{2} = (z +1)x+t. If subroutine P were to always start with α_{1} = 1, then on this example, for any p it would always start with hitting a zero divisor and then give up in step 4. This would cause subroutine M to loop.

Reconstruction of the parameters in subroutine P

We explain the choice of representation used by subroutine P for the coefficients of g and how they are reconstructed. Suppose g = x + \sqrt{s}/(s − t^{2}) − t/(s − t^{2}) is the gcd of the inputs of subroutine P. Three possible choices for canonical representations are:

(1) x + \sqrt{s}/(s − t^{2}) − t/(s − t^{2}),

(2) (s − t^{2})x + \sqrt{s} − t, and

(3) (\sqrt{s} + t)x + 1. Number (1) is the monic gcd g ∈ ℤ_{p}(s, t)[\sqrt{s}, x]. (2) is ǧ ∈ ℤ_{p}[s, t][\sqrt{s}][x]. It is obtained from (1) by clearing fractions. If there remains a variable that does not appear in m(z) then (3) is obtained from (2) by dividing (2) by the gcd of the coefficients in x, i.e., by gcd(s − t^{2}, \sqrt{s} − t) = \sqrt{s} − t.

We prefer (2) to (1) because we avoid arithmetic with fractions in s, t and can compute it using only univariate rational function reconstruction and univariate gcd computation. Although form (3) may be smaller than (2), computing it requires recursive gcd computations with one less variable which could be expensive. A precise definition of form (2) is now given. We require that

(i) g is a gcd of f_{1}, f_{2} (in characteristic p). (ii) g is in D_{p}[z, x] where D_{p} = ℤ_{p}[t_{1}, …, t_{k}]. (iii) g is primitive in D_{p}[z, x]. (iv) lc_{x}(g) is in D_{p}, i.e., deg_{z}(lc_{x} g) = 0, and

(v) (for uniqueness) lc_{t_{1},... ,t_{k}}(lc_{x} g) = 1 under a term ordering.

We show how subroutine P outputs form (2) using univariate rational function reconstruction. For p = 11, t = 1, suppose we computed the gcd at s = 2, 3, 4 with the Euclidean algorithm and have

g(2, 1) mod 11 = x + 1z − 1, g(3, 1) mod 11 = x − 5z + 5, g(4, 1) mod 11 = x + 4z − 4.

In step 7 subroutine P applies rational function reconstruction to the coefficients of z^{i}x^{j} to find rational functions in ℤ_{p}(s) of the form (as + b)/(cs + d). We obtain

h(s, 1) mod 11 = \frac{1}{1} x + \frac{1}{s − 1} z + \frac{−1}{s − 1} .

Now clear denominators to obtain form (2) at t = 1.

g(s, 1) mod 11 := (s − 1)h(s, 1) = (s − 1)x + z − 1.

We repeat this for t = 2, 3, 4 so that we can reconstruct rational functions in t of the form (at^{2} + bt + c)/(dt + e) for the coefficients of s^{i}z^{j}x^{k}. We obtain

h(s, t) mod 11 = ( \frac{1}{1} s + \frac{−t^{2}}{1} )x + \frac{1}{1} zs + \frac{−t}{1} .

There are no denominators in t to clear in our example. Thus we obtain (s − t^{2})x + z − t. Since the leading term, in lexicographical order with x > s > t, is 1sx, condition (v) already holds so no rescaling needs to be done.

Remark 2: In step 6 in subroutines M and P, there are two reasons we only combine images with the same leading term. First, this way we Chinese remainder only images that are scaled by condition (v) in a consistent way. Second, it prevents reconstruction problems arising from an unlucky content. Suppose that p_{1} is a good prime, ǧ = tx + t + p_{1}z and that subroutine M chooses prime p_{1}. Then ǧ mod p_{1} = tx + t has an unlucky content t, and g_{1} will be the monic associate of tx+t which is g_{1} = x+1. Hence g_{1} is not equal to ǧ mod p_{1}. Therefore, we must not Chinese remainder it with images g_{i} that are equal to ǧ mod p_{i}. The easiest way to prevent this is to Chinese remainder only those g_{i} together that have the same leading term. This way g_{1} will not contribute to the reconstruction of ǧ. Note that if lc_{t_{1},... ,t_{k}}(lc_{x}(ǧ)) vanishes mod a good prime p_{2}, then the image mod p_{2} will end up not being used even if there is no unlucky content. So our algorithm may compute images like these g_{1} or g_{2} that end up not being used. The reason we make no effort to avoid (as in Brown’s algorithm) such images before they have been computed is because that would require additional computations and slow down our algorithm (a simple check on lc_{t}(lc_{x}(f_{i})) without computing some resultant (mod p) is not sufficient, take for example f_{1} = f_{2} = (t^{2} + p_{2} z)x + 1 where z = \sqrt{t^{7}}).

Remark 3: The algorithm as stated assumes m(z) is irreducible so L is a field. Our implementation does not check this assumption because that would be expensive, and in many applications where L is known to be a field, unnecessary. There are applications for gcds in L[x] when L is not a field. In this case it is possible that our algorithm as stated will loop. Take for example m(z) = z^{2} − t^{2} and f_{1} = f_{2} = (z − t)x + 1. The algorithm can be modified to handle this case. The basic idea is to reconstruct the zero divisor from some of its images while we reconstruct the gcd, and terminate when trial division confirms either a zero divisor or a gcd.

The Trial Divisions

In step 9 in subroutine P the trial divisions g|f_{1} and g|f_{2} take place in characteristic p. Since g, f_{1} and f_{2} are in D_{p}[z]/⟨m(z)⟩[x] where D_{p} = ℤ_{p}[t_{1}, …, t_{k}] and m(z) ∈ D_{p}[z] we can use pseudo-division in D_{p}[z][x] modulo m(z) where we divide by m(z) also using pseudo-division. Thus we can do the trial divisions without any gcd computation in D_{p}. However, this results in a linear growth in size of the coefficients in D_{p}. In [5] we presented an algorithm for doing trial divisions of polynomials in ℤ[z][x] modulo m(z) ∈ ℤ[z] which uses pseudo-division and some gcds in ℤ to minimize growth of the integer coefficients. The same algorithm can be applied here; the only difference is the coefficient ring, D_{p} instead of ℤ. Our algorithm in [5] can also be used for the trial divisions in subroutine M. There f_{1}, f_{2} and h are in D[z]/⟨m(z)⟩[x] where D = ℤ[t_{1}, …, t_{k}]. Again, we do pseudo-division in D[z][x] modulo m(z) using gcds in D. Because we have effective algorithms for gcd computation in D and D_{p}, this approach is effective in practice. To prevent the algorithm from doing many trial divisions, we build into the rational reconstruction algorithm some redundancy so that if it succeeds with output h then ȟ = ǧ with high probability.

Rational Function Reconstruction

Let n, d ∈ ℤ_{p}[t] with d ≠ 0, lc_{t}(d) = 1, and gcd(n, d) = 1. Let m = m_{c}(t) be the modulus in subroutine P and (assuming gcd(d, m) = 1) let u = n/d mod m. The rational reconstruction problem is: given u and m find n and d. The problem is solved by the Euclidean algorithm as follows (this is Wang’s algorithm for ℚ (see [7]) modified for ℤ_{p}(t)).

Algorithm Rational Reconstruction (RR) Input: Polynomials m, u ∈ ℤ_{p}[t] with deg_{t} m > 0. Output: Either polynomials a, b ∈ ℤ_{p}[t] with deg_{t} a + deg_{t} b < deg_{t} m, lc_{t}(b) = 1, gcd(a, b) = 1, gcd(b, m) = 1 and a/b ≡ u mod m, or FAIL if no such a/b exists.

Set M = deg_{t} m, N = ⌊M/2⌋, D = M − N − 1.
Set (r_{0}, s_{0}) = (m, 0).
Set (r_{1}, s_{1}) = (u, 1).

While deg_{t} r_{1} > N do Set q to be the quotient of r_{0} divided r_{1}. Set (r_{0}, r_{1}) = (r_{1}, r_{0} − q r_{1}). Set (s_{0}, s_{1}) = (s_{1}, s_{0} − q s_{1}).

Set (a, b) = (r_{1}, s_{1}).

If deg_{t} b > D or gcd(b, m) ≠ 1 then output FAIL. Make b monic: Set u = lc_{t} b and (a, b) = (a/u, b/u). Output (a, b).

If N ≥ deg_{t}(n) and D ≥ deg_{t}(d) then algorithm RR outputs (n, d), otherwise, it outputs FAIL or some (a, b) ≠ (n, d).

In each loop of subroutine P, the degree of m_{c}(t_{k}) increases by 1. In step 7 in P, algorithm RR is applied to each coefficient of c until it fails. Thus algorithm RR will be applied many times to many coefficients before it succeeds. Suppose c has T terms in t_{1}, …, t_{k−1}, z, x and suppose subroutine P requires S points to succeed. Then the expected number of calls to algorithm RR is ST. Moreover, the trial divisions in subroutine P will be attempted many times with h ≠ ǧ.

We force algorithm RR to output FAIL with high probability when M = deg_{t}(m) is not large enough for RR to reconstruct a coefficient of h by setting N = ⌊(M − 1)/2⌋ and D = M − N − 2 so that N + D = M − 2, i.e., deg_{t}(m) = N + D + 2. In the Euclidean algorithm, since deg_{t}(q) + deg_{t}(r_{1}) + deg_{t}(s_{1}) = deg_{t}(m) the new settings of N and D mean we are forcing the q in algorithm RR that corresponds to the output to be of degree 2 or more. For u chosen at random in ℤ_{p}[t] of degree less than M, the probability of getting any q of degree 2 or more is O(M/p). Thus, by requiring one evaluation point more than the minimum needed, with high probability, we attempt the trial divisions once, and secondly, we reduce the expected number of calls to algorithm RR from O(ST) to O(S + T).

2.3 Multivariate Inputs

Let f_{1}, f_{2} ∈ L[x_{1}, ..., x_{n}] with n > 1. To compute the gcd g of f_{1} and f_{2} using our modular GCD algorithm do the following: write f_{1} = ∑ a_{i}x_{1}^{i} and f_{2} = ∑ b_{j}x_{1}^{j} where a_{i}, b_{j} ∈ L[x_{2}, …, x_{n}]. Now recall that gcd(f_{1}, f_{2}) = cb where c is the gcd of the contents in x_{1} of f_{1} and f_{2}, i.e., c = gcd(a_{i}, b_{j}), and b is the gcd of h_{1} = f_{1}/c and h_{2} = f_{2}/c. Computing c requires gcd computations in L[x_{2}, …, x_{n}], that is, in one less variable. To compute gcd(h_{1}, h_{2}) we view x_{2}, …, x_{n} as parameters, that is, we write h_{1}, h_{2} as polynomials in K[x_{1}] where K = G[z]/⟨m(z)⟩ and G = ℚ(t_{1}, …, t_{k}, x_{2}, …, x_{n}). Now we compute ǧ, the monic associate of their gcd, using our modular GCD algorithm. When x_{2}, …, x_{n} are regarded as polynomial variables, ǧ may have a non-trivial content in x_{1} which needs to be computed and divided out to obtain b. That ǧ could have a non-trivial content is illustrated by the example (s − t^{2})x + \sqrt{s} − t given in section 2.2.

3. FRACTION-FREE ALGORITHMS

In this section we develop a fraction-free GCD algorithm for L[x] based on the ideas of Moreno Maza and Rioboo in [6]. Their algorithm, a modification of the subresultant GCD algorithm, computes an associate of the gcd of two univariate polynomials modulo a triangular set of polynomials over an integral domain D. Our problem setting is a special case of theirs where we have a triangular set of one polynomial, m(z) = m̌(z), which is irreducible over ℤ[t_{1}, …, t_{k}]. We recall the idea of pseudo-division in D[x].

Definition: Let a(x), b(x) be non-zero polynomials in D[x] with deg a ≥ deg b. Let δ = deg a − deg b + 1. The pseudo remainder r̄ and pseudo quotient q̄ are the remainder and quotient, respectively, of μa divided by b where μ = lc_{x} b^{δ}. Thus μa = bq̄ + r̄ where r̄ = 0 or deg r̄ < deg b.

The point of pseudo-division is that r̄ and q̄ are in D[x] and no fractions appear in the division algorithm when dividing μa by b. The subresultant GCD algorithm on input of non-zero f_{1} and f_{2} in D[x], uses pseudo-division to compute an associate of the gcd of f_{1} and f_{2} in D[x]. It does this without introducing fractions and without gcd computation in D. Moreover, one can show that if D = ℤ[t_{1}, …, t_{k}] the size of the integer coefficients and degree in each parameter of the coefficients of the pseudo-remainder sequence (PRS) grows linearly which is is optimal to within a constant factor. This is achieved by dividing pseudo-remainders by known exact divisors in D.

In [6], Moreno Maza and Rioboo modify the subresultant PRS to work for R[x] where R = D[z]/⟨m(z)⟩. The main idea is to multiply a pseudo-remainder p(x) ∈ R[x] by a scalar i ∈ R, a quasi-inverse of lc_{x}(p), so that after multiplication, lc_{x}(ip) ∈ D and not D[z]. Now pseudo-division does not introduce fractions and the exact divisions (Moreno Maza and Rioboo prove that they remain exact) in the subresultant PRS are by elements of D not R. We recall the definition for quasi-inverse for our commutative ring R.

Definition: Let u ∈ R. Then v ∈ R is a quasi-inverse of u if uv = r for some r ∈ D.

Remark 4: The definition is unique up to multiplication by a non-zero element of D. To reduce coefficient growth we will want to use a quasi-inverse which is smallest.

Example: Let u ∈ R = ℤ[z]/⟨m⟩ where m ∈ ℤ[z] is monic. Then there exist polynomials s, t ∈ ℤ[z] such that sm + tu = r where r ∈ ℤ is the resultant of m and u. Thus v = t is a quasi-inverse of u in R. Let g = gcd(r, t). Now t ∈ ℤ[z] and so g is the gcd of r and the coefficients of t. Then t/g is also a quasi-inverse of u; it is a minimal quasiinverse.

To compute a quasi-inverse of u ∈ R = D[z]/⟨m(z)⟩ we can use either Collin’s reduced PRS, or the subresultant PRS, or the primitive PRS to compute elements s, t ∈ D[z] such that sm + tu = r for some r ∈ D and output t/g where g = gcd(r, t). Here is Maple code to accomplish this using the reduced PRS for D = ℤ[t_{1}, …, t_{k}].

                                                          QuasiInv := proc(x,m,z)
                                                          local u,r0,r1,t0,t1,pr,pq,mu,beta;
                                                            u := primpart(x,z);
                                                            r0,r1,t0,t1 := m,u,0,1;
                                                            beta := 1;
                                                            while degree(r1,z)>0 do
                                                              r0,r1 := r1,prem(r0,r1,z,'mu','pq');
                                                              t0,t1 := t1,expand(mu*t0-pq*t1);
                                                              if r1=0 then ERROR("hit zero divisor",r0) fi;
                                                              divide(r1,beta,'r1');
                                                              divide(t1,beta,'t1');
                                                              beta := mu;
                                                            od;
                                                            divide(t1,gcd(r1,t1),'t1');

    return t1;
  end;

Below is Maple code for the modified subresultant PRS for computing an associate of the gcd of f_{1}, f_{2} in D[z]/⟨m(z)⟩. Moreno Maza and Rioboo prove that the divisions by elements of D remain exact. Note, if m(z) is not monic in D[z] the divisions by m(z) (the Maple command rem(...,m,z);) may be replaced by pseudo-division by m(z) without breaking the exact division properties of the subresultant PRS.

  MMRPRS := proc(f1,f2,x,m,z)
  local p1,p2,p3,del,alp,psi,bet,i;
    p1,p2 := f1,f2;
    del := degree(p1,x)-degree(f2,x);
    if del<0 then return MMRPRS(f2,f1,x,m,z) fi;
    i := QuasiInv(lcoeff(p2,x),m,z);
    p2 := rem(i*p2,m,z);
    psi,bet,alp := -1,(-1)^(del+1),lcoeff(p2,x);
    while degree(p2,x)>0 do
      p3 := rem(prem(p1,p2,x),m,z);
      if p3=0 then return p2 fi;
      divide(p3,bet,'p3');
      i := QuasiInv(lcoeff(p3,x),m,z);
      p3 := rem(i*p3,m,z);
      if del>0 then
         divide((-alp)^del,psi^(del-1),'psi')
      fi;
      p1,p2 := p2,p3;
      if degree(p2,x)>0 then
         del := degree(p1,x)-degree(p2,x);
         bet := -alp*psi^del;
         alp := lcoeff(p2,x);
      fi;
    od;
    p2;
  end:

However, multiplication of p_{3} by a quasi-inverse i ∈ D[z] introduces a new coefficient growth in D into the subresultant PRS. We observe that the growth is approximately cubic instead of linear. We illustrate this. For D = ℤ, m = z^{3} + 80z^{2} + 23z − 20, we created two polynomials f_{1} and f_{2} of degree 21 and 20 in x respectively, with coefficients polynomials of degree 2 in z with integer coefficients chosen at random from (-100,100) in Maple as follows.

  > c := () -> randpoly(z,degree=2,dense):
  > f1 := randpoly(x,coeffs=c,degree=21,dense):
  > f2 := randpoly(x,coeffs=c,degree=20,dense):

The data below shows the degree (column 1) and length in decimal digits of the largest integer coefficient of the polynomials p_{3} appearing in the Moreno Maza and Rioboo PRS (columns 2 and 3) and the primitive PRS (columns 4 and 5). Columns 2 and 3 show, respectively, the length of the largest integer coefficient of p_{3} before and after the quantity bet is divided out. Let r_{3} denote the pseudo remainder of the primitive parts of p_{1} and p_{2}. Columns 4 and 5 show, respectively, the length of the largest integer coefficient of r_{3} and the primitive part of r_{3}.

The data shows that after 20 steps the length of the integers in the Moreno Maza and Rioboo algorithm is over 100 times longer than those in the primitive PRS. We see also that after an initial step, the size of the integers in column 5

2  3  4  5
19  21  21  21 20
18  96  79  62 53
17  308  221 121 95
16  760  497 182 135
15  1572  945 241 174
14  2867 1609 303 216
13  4763 2527 365 257
12  7386 3742 428 299
11 10854 5293 490 342
10 15302 7229 554 384
9 20757 9592 619 427
8 27638 12417 678 465
7 35777 15750 743 509
6 45390 19630 807 553
5 56614 24105 872 596
4 69576 29218 936 639
3 84409 35014 1003 684
2 101248 41538 1072 729
1 120220 48826 1132 766
0 141447 56925 1193 810

Table 1: Integer Coefficient Growth

grows each step by approximately 42 digits and the size of those in column 4 by approximately 62 digits. This is a linear growth. In columns 2 and 3 the growth is approximately cubic. There are no parameters present in our example. If there were a parameter t present in m and f_{1} and f_{2}, what we would see is a similar growth in the degree in t, i.e., there would be a two dimensional growth.

If instead of modifying the subresultant PRS we modify the primitive PRS to compute a gcd in R[x] by making the pseudo remainders primitive in D[x, z] then we get minimal (linear) coefficient growth. To do this we assume we can compute gcds in D effectively. The same idea can be applied to the general setting of the algorithm of Moreno Maza and Rioboo in [6] for gcd computation modulo a triangular set; instead of using the (extended) subresultant PRS use the (extended) primitive PRS throughout. Here is a Maple code to do this for D = ℤ[t_{1}, …, t_{k}].

                                                          PrimitivePRS := proc(f1,f2,x,m,z)
                                                          local i,r1,r2,r3;
                                                             r1 := primpart(f1,[x,z]);
                                                             r2 := primpart(f2,[x,z]);
                                                             while r2 <> 0 do
                                                                i := QuasiInv(lcoeff(r2,x),m,z);
                                                                r2 := primpart(prem(i*r2,m,z), [x,z]);
                                                                r3 := prem(prem(r1,r2,x),m,z);
                                                                r3 := primpart(r3,[x,z]);
                                                                r1,r2 := r2,r3;
                                                             od;
                                                             r1;
                                                          end;

The number of gcds in D = ℤ[t_{1}, …, t_{k}] that are computed when taking the primitive part of a polynomial in D[x, z] of degree n in x and d in z can be effectively reduced from O(nd) to 1 by computing the gcd of the smallest coefficient with a random linear combination of the others. Thus the expected total number of gcds in D that the primitive PRS does is O(δ) where δ = max(deg_{x} f_{1}, deg_{x} f_{2}) − deg_{x} g + 1.

4. IMPLEMENTATION

We compare a Maple implementation of our algorithms on the following problem set. Let

m(z) = z^{3} − (5 − t)z^{2} + (7 − t^{2})z − (9 − t^{3}),

g = (10 − 4 t)x^{2} − 5 xz^{2} + (4 t + 1)xz + (11 − 17 t^{2} + 9 t)x −19 z^{2} + (−7 t + 6)z + (−11 t^{2} + 15 t + 3),

a = (18+10 t)x^{2} +10 xz^{2} +(17 t+2)xz +(2+17 t^{2} +8 t)x +6 z^{2} + (17 t + 6)z + (4 t^{2} − 4 t + 2),

b = (−8−11 t)x^{2}−14 xz^{2}+(8 t−4)xz+(−17−5 t^{2}+19 t)x −11 z^{2} + (17 t − 4)z + (−14 t^{2} − 19 t − 2).

For n = 10, k = 0, 1, …, n let f_{1} = g^{k}a^{n−k} and f_{2} = g^{k}b^{n−k}. Thus there are 11 gcd problems, each of degree 20 in x and 20 in t where the gcd g^{k} has degree 2k in x and degree 2k in t with leading coefficient degree k in t. The reason for this choice of inputs is to compare the algorithms as the size of the gcd increases relative to the size of the cofactors. The modular GCD algorithm will do best for small k and the PRS algorithms will do best for large k.

The timings below (in CPU seconds) were obtained using Maple 9 on a Itanium PC running Linux. Column 1 is the time for the Moreno Maza Rioboo algorithm. Column 2 is the time for our primitive PRS algorithm. Column 3 is the time for our modular GCD algorithm (shown in parentheses is the number of primes used). Column 4 is the time for an improved version of the modular GCD algorithm.

k  1  2  3  4
0  NA  NA 0.07 (1) 0.06
1  NA  NA 0.77 (2) 0.22
2  NA  NA 1.56 (3) 0.54
3  NA  NA 2.86 (3) 0.84
4  NA  NA 4.39 (5) 1.24
5  NA 7459.7 5.43 (6) 1.99
6  NA 1488.7 7.70 (8) 2.61
7  NA 259.6 8.01 (9) 3.34
8 4473.1  42.8 7.70 (10) 4.18
9  31.1  1.2 7.22 (11) 5.07
10  0.0  0.0 6.29 (12) 5.72

Table 2: NA means not attempted

The data demonstrates clearly the superiority of our modular GCD algorithm. Much of the time in our implementation is spent doing trial divisions in subroutine P, that is, in the ring ℤ_{p}(t)[z]/⟨m(z)⟩[x]. Column 4 shows the improvement we obtained by making the following change to our modular GCD algorithm. We modify subroutine P so that if rational function reconstruction succeeds we assume the result is correct and simply return h without trial division. Thus when subroutine P calls itself recursively on input A_{i} = f_{1}(t_{k} = α_{i}), B_{i} = f_{2}(t_{k} = α_{i}), and m_{i} = m_{c}(t_{k} = α_{i}), it is now possible that an output g_{i} from subroutine P is wrong, i.e., g_{i} does not divide A_{i} or g_{i} does not divide B_{i} modulo m_{i}(z). We apply the following strategy to detect an incorrect g_{i} so that subroutine P will eventually terminate. In subroutine P, if k > 1, periodically, we choose β = (β_{1}, …, β_{k−1}) ∈ ℤ_{p}^{k−1} at random s.t. lc_{x}(g_{i})(β) ≠ 0 and lc_{z}(m_{i})(β) ≠ 0. We test whether g_{i}(β)|A_{i}(β) and g_{i}(β)|B_{i}(β) modulo m_{i}(β, z). Thus we replace division over a function field in k − 1 parameters with division over a finite ring. With high probability this check will identify an incorrect g_{i}. Once identified, we discard g_{i} in subroutine P and adjust the c and m_{c} accordingly. We modify subroutine M similarly to detect incorrect output from subroutine P. It remains to say how often we make these checks. At step n, after computing g_{n}, if 2^{j}|n − 1 for some j > 2 we check g_{n−2^{j}} , g_{n−2^{j} +1}, …, g_{n−1} for correctness. Thus the first time we make any checks is at step n = 9 where we check g_{1}, g_{2}, …, g_{8}. The next time is at step 17 when we check g_{1}, g_{2}, …, g_{16}. The next is at step 24 when we check g_{17}, …, g_{24}. Thus we check each g_{i} within 8 steps of being computed and then with decreasing frequency. In practice, because of the way we have designed the rational function reconstruction, the output from subroutine P is hardly ever wrong, hence, the infrequent checks are justified.

5. REFERENCES

[1] W. S. Brown (1971). On Euclid’s Algorithm and the Computation of Polynomial Greatest Common Divisors, J. ACM 18, pp. 476–504.

[2] G. E. Collins and M. J. Encarnacion (1995). Efficient Rational Number Reconstruction. J. Symbolic Computation 20, pp. 287–297.

[3] M. J. Encarnacion (1995). Computing GCDs of Polynomials over Algebraic Number Fields, J. Symbolic Computation 20, pp. 299–313.

[4] J. von zur Gathen and J. Gerhard (1999). Modern Computer Algebra. University of Cambridge Press.

[5] M. van Hoeij, M. B. Monagan, A Modular GCD Algorithm over Number Fields Presented with Multiple Field Extensions. Proceedings of ISSAC ’2002, ACM Press, pp. 109–116, 2002.

[6] M. Moreno Maza, R. Rioboo (1995). Polynomial Gcd Computations over Towers of Algebraic Extensions, Proc. of AAECC-11 Springer-Verlag LNCS 948 (1995), pp. 365–382.

[7] P. Wang (1981). A p-adic Algorithm for Univariate Partial Fractions. Proceedings of SYMSAC ’81, ACM Press, pp 212-217.