A Modular GCD algorithm over Number Fields presented with Multiple Extensions.

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

We consider the problem of computing the monic gcd of two polynomials over a number field L = ℚ(α_{1}, …, α_{n}). Encarnacion, Langemyr and McCallum have already shown how Brown’s modular GCD algorithm for polynomials over ℚ can be modified to work for ℚ(α).

Our first contribution is an extension of Encarnacion’s modular GCD algorithm to the case n > 1 without converting to a single field extension. Our second contribution is a proof that it is not necessary to test if p divides the discriminant. This simplifies the algorithm; it is correct without this test.

Our third contribution is the design of a data structure for representing multivariate polynomials over number fields with multiple field extensions. We have a complete implementation of the modular GCD algorithm using it. We provide details of some practical improvements.

1. INTRODUCTION

We recall the relevant details of the so called modular GCD algorithm first developed by Brown in [3] for polynomials over ℤ and then by Langemyr and McCallum in [7] and Encarnacion in [4] for polynomials over L = ℚ(α), which we shall generalize to L = ℚ(α_{1}, …, α_{n}). First some notation.

We denote the input polynomials by f_{1} and f_{2}, their monic gcd by g. The cofactors are the polynomials f_{1}/g and f_{2}/g. If f ∈ ℚ[x] then let the denominator den(f) be the smallest positive integer such that den(f)f ∈ ℤ[x]. See section 2 for the definition of den(f) if f ∈ ℚ(α_{1}, …, α_{n})[x]. The associate f̃ of f is defined as f̃ = den(g)g where g = monic(f). Here monic(f) is defined as lc(f)^{−1}f where lc(f) is the leading coefficient of f. Define the semi-associate f̌ as rf

∗Supported by NSF grant 0098034. †Supported by NSERC of Canada and the MITACS NCE of Canada. where r is the smallest positive rational number for which den(rf) = 1.

Computing the associate f̃ is useful for removing denominators, but could be expensive if lc(f) is a complicated algebraic number. So we preprocess the input polynomials in our algorithm by taking the semi-associate instead. If lc(f) ∈ ℚ then the two notions are the same up to a sign:

f̌ = ±f̃ ⟺ lc(f) ∈ ℚ

Examples: If f = 2x−2/3 then f̌ = f̃ = 3x−1. If α = \sqrt{2} and f = −αx + 1 then f̌ = f, monic(f) = x − α/2 and f̃ = 2x − α.

The modular GCD algorithm computes the monic gcd g ∈ L[x] of f_{1} and f_{2}. It does this by reducing f_{1}, f_{2} modulo one or more primes and calling the Euclidean algorithm mod p for each of these primes p. If p is a good prime, the Euclidean algorithm mod p returns g mod p. The modular GCD algorithm reconstructs g from these modular images. Only good primes should be used during the reconstruction for it to be successful. However, not all primes are good. We distinguish the following cases:

Definition 1. Let f_{1}, f_{2} ∈ L[x] and g be their monic gcd. We will distinguish four types of primes.

• lc-bad primes. Let m_{1}, …, m_{n} be the minimal polynomials of the field extensions α_{1}, …, α_{n}. So m_{i}(z) is a monic irreducible polynomial in ℚ(α_{1}, …, α_{i−1})[z] and m_{i}(α_{i}) = 0. If den(f_{1}), den(f_{2}) or any leading coefficient of f̌_{2}, m̌_{1}, …, m̌_{n} vanishes mod p then we call p an lc-bad prime.

• Fail primes. If p is not an lc-bad prime, and the Euclidean algorithm mod p returns “failed”, then p is called a fail prime.

• Unlucky primes. If p is not an lc-bad prime nor a fail prime, and if the output of the Euclidean algorithm mod p has higher degree than g, then p is called an unlucky prime.

• Good primes. A prime p is called a good prime if the Euclidean algorithm mod p returns g mod p. Theorem 1 in section 2 says that all primes that are not lc-bad are either fail, unlucky or good.

Remarks:

1. Our definition of lc-bad prime is not symmetric in f_{1}, f_{2}. It could be that p is lc-bad for f_{1}, f_{2} but not lc-bad for f_{2}, f_{1}. In that case, because of how we set up the algorithm, we should either: not use p, or: interchange f_{1}, f_{2} mod p before calling the Euclidean algorithm mod p.

2. Our definitions are not the same as the standard definitions in [3]. For example, it is possible that the Euclidean algorithm mod p fails even if the monic gcd of f_{1} mod p, f_{2} mod p exists and equals g mod p. We call such p a fail prime and not a good prime. This distinction is not necessary if f_{1}, f_{2} ∈ ℚ[x] where there are no fail primes.

3. If p | den(g) (in the standard definition these primes are called bad primes) then g mod p is not defined and so p can not be a good prime. According to theorem 1, p must then be either lc-bad, fail, or unlucky.

4. Minimal polynomials are monic so the leading coefficients of m̌_{1}, …, m̌_{n} are den(m_{1}), …, den(m_{n}) ∈ ℤ. However, lc(f̌_{2}) is in general not an integer but an algebraic number.

5. It is very easy to tell if a prime p is lc-bad or not, but we can not tell in advance if p is fail, unlucky, or good. So we will end up calling the Euclidean algorithm mod p with fail, unlucky, and good primes but never with lc-bad primes.

1.1 lc-bad primes

If f_{1} = 5x + 1, f_{2} = 5x − 1 and p = 5 then p satisfies our definition of an lc-bad prime as well as the definition of a good prime. However, there are good reasons not to use any lc-bad prime. Take for example f_{1} = f_{2} = 5x + 1. Also, the proof of theorem 1 requires that p not be lc-bad.

Another example is L = ℚ(α), f_{1}, f_{2} ∈ L[x] with gcd g = x + α^{3}, p = 5, and the minimal polynomial of α is m = z^{5} +z^{4} + \frac{1}{5} z^{3} − \frac{1}{5} . Because of preprocessing, in the algorithm we work with m̌ = 5z^{5} + 5z^{4} + z^{3} − 1. Modulo p = 5 this becomes z^{3} + 4. If we used the prime p = 5, it is easy to give an example f_{1}, f_{2} where the Euclidean algorithm mod p returns g mod (5, α^{3} + 4) which is x + 1. But, viewing α as a variable, g ≢ x + 1 mod 5.

For our algorithm, the best solution to the above problems is: never use an lc-bad prime.

1.2 Fail primes

Fail primes are primes for which the Euclidean algorithm mod p tries to divide by a zero divisor, in which case it returns “failed”. Take for example f_{1} = x^{2} − 1, f_{2} = ax − a where a = 2^{1/5} + 5. Denote a mod p as \overline{a}. The Euclidean algorithm mod p will first try to make f_{2} mod p monic by multiplying it with 1/\overline{a}. But if N(a), the norm of a, vanishes mod p then \overline{a} is zero or a zero-divisor, and the computation of 1/\overline{a} fails. In this example N(a) = 53·59 so the fail primes are 53 and 59.

The reason that in our terminology 53 and 59 are called fail primes and not lc-bad primes in the example (after all, the problem was caused by lc(f_{2}) mod p) is to indicate how these primes are discarded: We do not actively avoid these primes, instead, they “discard themselves” when the Euclidean algorithm mod p is called.

One can also construct examples where p is not lc-bad, lc(f_{2}) is a unit mod p, but p still divides den(g) (occasionally such p can be unlucky instead of fail). Take for example α with minimal polynomial m = z^{3} + 3z^{2} − 46z + 1, f_{1} = x^{3}−2x^{2}+(−2α^{2}+8α+2)x−α^{2}+11α−1, f_{2} = x^{3}−2x^{2}−x+1. The monic gcd is g = x − \frac{1}{91} α^{2} − \frac{23}{91} α − \frac{50}{91} . The denominator is den(g) = 91 = 7 · 13. In this example, if p ∈ {7, 13} then p is not lc-bad and the leading coefficient of f_{2} (as well as of f_{1}) is a unit mod p. Nevertheless, p can not be a good prime because p | den(g). In this type of example p must divide the discriminant. For this reason, Encarnacion [4] tests if the discriminant is 0 mod p and avoids such primes. However, even without the discriminant-test, the primes p ∈ {7, 13} would still have been discarded at some point: The Euclidean algorithm mod p will calculate r_{3} = f_{1} mod (p, f_{2}), try to make r_{3} monic and fail because the leading coefficient of r_{3}, namely, −2α^{2} + 8α + 3, is a zero divisor mod p.

Although one can generalize the discriminant-test to L, see section 3 in [6], our algorithm does not use it because it makes no difference for the correctness of the algorithm. For an intuitive explanation see lemma 4 and for a proof see theorem 1.

1.3 Unlucky primes

Unlucky primes are not trivially detectable like lc-bad primes and do not “discard themselves” like fail primes do, but need to be detected and discarded nevertheless. Fortunately, Brown [3] showed how to do this in a way that is efficient and easy to implement: Whenever modular gcd’s do not have the same degree, keep only those of smallest degree and discard the others.

As an example, take f_{1} = x^{2}+(2\sqrt{5}+1)x+3, f_{2} = x^{2}−x−1, g = x + (\sqrt{5} − 1)/2. Then the Euclidean algorithm mod 2 will return x^{2} + x + 1, so p = 2 is an unlucky prime. But if f_{1} = x^{2} + \sqrt{5} x + 1, f_{2} and g the same as before, then p = 2 is a fail prime.

1.4 Good primes

All but finitely many primes must be good. This is because if one would run the Euclidean algorithm in characteristic 0, it would be a finite computation, and so there can only be finitely many conditions on the primes and each condition only excludes finitely many primes (see lemma 5).

Of course we will not run the Euclidean algorithm in characteristic 0, so this does not tell us which primes to use. But this is not a problem because to guarantee correctness of the algorithm, just as in Brown’s algorithm, all we need to do is to avoid the lc-bad primes. Experiments show that random primes are good with high probability. Hence, even if there was an oracle that quickly provided good primes, it would not noticeably improve the running time.

1.5 Motivation for the algorithm

The goal of this paper is to present an efficient modular GCD algorithm over a field L that consists of multiple extensions over ℚ. Suppose the largest numerator or denominator in g is c. To reconstruct g by computing g mod P = p_{1} ⋯ p_{m} using primes p_{1}, …, p_{m}, if we want log(P) = O(log(c)), that is, if we want the number of primes used to be proportional to the size of the coefficients in g, then we are forced to

1. Not use a primitive element to convert to a single extension, which is expensive and can cause a blowup in the size of the coefficients. This problem is well known, e.g. see [1].

2. Not invert lc(f_{2}), which can also cause a blowup, and can also be more expensive than computing g.

3. Use rational reconstruction. Otherwise a denominator bound would be necessary, but such bounds are generally too large. The defect bound (usually the (reduced [6]) discriminant), which is part of the denominator bound, is usually also too large.

Encarnacion’s paper confirms and deals with these items. As a result, Encarnacion’s algorithm is the fastest algorithm for a single extension. As for item 1, his paper deals only with a single extension, but he does illustrate that modifying that extension (making α_{1} an algebraic integer) is not efficient. But if modifying one extension α_{1} is not efficient, then modifying n extensions (replacing it by a primitive element) is certainly not efficient.

Our goal is to generalize Encarnacion’s algorithm to multiple extensions, without using a primitive element. A technical difficulty is how to generalize the discriminant test, which we did a preprint [6], where we also gave a formula for the reduced discriminant. It turned out, however, that a discriminant test is not necessary, it can be omitted.

Omitting the discriminant test simplifies the algorithm but has no noticeable impact on the running time. So for the single extension case, our algorithm is essentially the same as Encarnacion’s algorithm. For the multiple extension case our algorithm provides an obvious efficiency improvement because the primitive element conversion alone can easily cost more time than our entire algorithm, especially if the gcd is small.

In this paper we only treat univariate polynomials f_{1}, f_{2} ∈ L[x], but our implementation handles the multivariate case as well.

2. THE EUCLIDEAN ALGORITHM OVER A RING

Let α_{1}, …, α_{n} be algebraic numbers. Let L_{i} = ℚ(α_{1}, …, α_{i}) and L = L_{n}. Let d_{i} be the degree of α_{i} over L_{i−1}. The dimension of L as a ℚ-vector space is d_{∗} := d_{1} ⋯ d_{n}. A basis of L is:

M := {∏_{i=1}^{n} α_{i}^{e_{i}} | 0 ≤ e_{i} < d_{i}}.

Let R̃ be the set of all ℤ-linear combinations of M and let R̃_{i} = R̃ ⋂ L_{i}. Let m_{i} be the minimal polynomial of α_{i} over L_{i−1}. The degree of m_{i} is d_{i}, m_{i} is monic (the leading coefficient is lc(m_{i}) = 1) and m_{i}(α_{i}) = 0. The coefficients of m_{i} are in L_{i−1}. Let l_{i} be the smallest positive integer such that the coefficients of l_{i}m_{i} are in R̃_{i−1}. Denote 𝔽_{p} = ℤ/pℤ and l_{∗} = l_{1} ⋯ l_{n}.

In general R̃ is not a ring. For example, α_{1} ∈ R̃, but α_{1}^{d_{1}} is not in R̃ unless l_{1} = 1. When a, b ∈ R̃, to compute the product ab ∈ L we replace α_{1}, …, α_{n} by variables z_{1}, …, z_{n}, then multiply a,b as polynomials, and after that take the remainder modulo the polynomials m_{1}(z_{1}), …, m_{n}(z_{n}). During this computation we only divide a bounded number of times by l_{1}, …, l_{n}. Hence, if k is a sufficiently large integer, then l_{∗}^{k}ab ∈ R̃ for all a, b ∈ R̃.

If a ∈ L then define the denominator of a as the smallest positive integer den(a) such that den(a)a ∈ R̃. Note that R̃, and hence den(a), depends on the choice of α_{1}, …, α_{n}. For example, if α_{1} = \sqrt{8} and a = \frac{1}{2} α_{1} then den(a) = 2. For a ∈ L one has a ∈ R̃ ⟺ den(a) = 1, in particular den(0) = 1. Define

R_{p} = {a ∈ L | den(a) ≢ 0 mod p}  (1)
= \{ \frac{a}{m} | a ∈ R̃, m ∈ ℤ, m ≢ 0 mod p\}. (2)

If a, b ∈ L then den(ab) divides den(a)den(b)l_{∗}^{k} for some k. Hence, if p ∤ l_{∗} then R_{p} is a ring. We will always assume that p does not divide l_{∗} so that R_{p} is a ring (if p | l_{∗} then p is an lc-bad prime). Denote

ℤ_{(p)} = R_{p} ⋂ ℚ = \{ \frac{a}{m} | a, m ∈ ℤ, m ≢ 0 mod p\}.

Then R_{p} is a ℤ_{(p)}-module with basis M. Define

\overline{R} = R_{p}/pR_{p}.

If a ∈ R_{p} then we use the notation \overline{a}, or also a mod p, for the image of a in \overline{R}. If a ∈ L, then (primes that divide l_{∗} are always excluded)

\overline{a} is defined ⟺ a ∈ R_{p} ⟺ p ∤ l_{∗}den(a).

If \overline{a} is defined we will say that a can be reduced mod p.

Now \overline{R} is a ring and also an 𝔽_{p}-vector space with basis M mod p. We can do the following identifications:

R_{p} = R̃ ⊗_{ℤ} ℤ_{(p)}, L = R̃ ⊗_{ℤ} ℚ, and \overline{R} = R̃ ⊗_{ℤ} 𝔽_{p} (3)

If a ∈ L then a is a unit in R_{p} if and only if both a and 1/a are in R_{p} (whenever we write 1/a it is implicitly assumed that a ≠ 0). This is equivalent to p ∤ l_{∗}den(a)den(1/a). If a ∈ L we will call a a unit mod p if a ∈ R_{p} and \overline{a} is a unit in \overline{R}. The following lemma shows that these two notions are equivalent.

Lemma 1. Let a ∈ R_{p}. Then a is a unit in R_{p} if and only if \overline{a} is a unit in \overline{R}.

Proof: If a is a unit in R_{p} then a and 1/a are in R_{p}, hence \overline{a} and \overline{1/a} are defined, and since a ↦ \overline{a} is a ring homomorphism R_{p} → \overline{R} one sees that \overline{1/a} is the inverse of \overline{a}. Hence \overline{a} is a unit in \overline{R}. Conversely, assume \overline{a} is a unit. Then a ≠ 0 so we can take b := 1/a ∈ L. To finish the proof we need to show that b ∈ R_{p}. Take the smallest integer k for which c := bp^{k} ∈ R_{p}. Since k is minimal, we have \overline{c} ≠ 0 but then \overline{a}\overline{c} is the product of a unit and a nonzero element in \overline{R} and hence nonzero. But \overline{a}\overline{c} equals \overline{abp^{k}} = \overline{p^{k}} so \overline{p^{k}} ≠ 0, hence k = 0, so b ∈ R_{p} and a is invertible in R_{p}.

If f ∈ L[x] then the denominator den(f) is defined as the smallest positive integer such that den(f)f ∈ R̃[x]. Now f ∈ R_{p}[x] if and only if p ∤ den(f)l_{∗}. The polynomial \overline{f} is the image of f in \overline{R}[x], and is defined if and only if f ∈ R_{p}[x], in which case we will say that f can be reduced mod p. Furthermore, if f and \overline{f} have the same degree (when lc(f) is nonzero mod p) then we will say that f reduces properly mod p. If p is not an lc-bad prime it means that f_{1}, f_{2} can be reduced mod p, and that f_{2} reduces properly mod p.

Let 0 ≤ i ≤ j ≤ n and a ∈ L_{j}. Multiplication by a is an L_{i}-linear map ψ : L_{j} → L_{j}. The characteristic polynomial cp_{i}^{j}(a) ∈ L_{i}[x] of a over the extension L_{j} : L_{i} is defined as the characteristic polynomial of this linear map. The trace Tr_{i}^{j}(a) of a over L_{j} : L_{i} is the trace of ψ and the norm N_{i}^{j}(a) of a over L_{j} : L_{i} is the determinant of ψ. Whenever we do not mention the extension L_{j} : L_{i} it is assumed to be L : ℚ (so i = 0 and j = n) in which case we write Tr(a), N(a), cp(a). Now the integral closure of ℤ in L is

O = {a ∈ L | cp(a) ∈ ℤ[x]}.

This is a ring (see [5]), and the elements of O are called the algebraic integers in L. We will use the following notation for the integral closure of ℤ_{(p)} in L

O_{p} = {a ∈ L | cp(a) ∈ ℤ_{(p)}[x]}.

Suppose a ∈ L and m = den(cp(a)). Then by definition a ∈ O_{p} if and only if m ≢ 0 mod p. The characteristic polynomial of ma is in ℤ[x], hence ma ∈ O and hence

O_{p} = \{ \frac{a}{m} | a ∈ O, m ∈ ℤ, m ≢ 0 mod p\}.  (4)

Lemma 2. If 0 ≤ i ≤ j ≤ n and a ∈ O_{p} ⋂ L_{j} then a is a unit in O_{p} if and only if N_{i}^{j}(a) is a unit in O_{p}. In particular, a ∈ O_{p} is a unit if and only if N(a) ∈ ℚ is a unit in ℤ_{(p)}, in other words, both numerator and denominator of N(a) are not divisible by p. The same is also true for R_{p}.

Remark: If p ∤ l_{∗} then R_{p} ⊆ O_{p} and the lemma implies that if a ∈ R_{p} and 1/a ∈ O_{p} then 1/a ∈ R_{p}.

Proof: The L_{i}-linear map ψ : L_{j} → L_{j} that corresponds to multiplication by a is defined over O_{p}, i.e. the entries of the matrix of ψ are in O_{p}. If N_{i}^{j}(a), the determinant of ψ, is a unit in O_{p} then the matrix is invertible over O_{p}. So then ψ^{−1}(1) ∈ O_{p}, so 1/a ∈ O_{p}. Conversely, if a is invertible in O_{p} then ψ is an invertible linear map, so its determinant must be a unit. Now N(a) = N_{0}^{n}(a) ∈ L_{0} = ℚ and ℚ ⋂ O_{p} = ℤ_{(p)} so the second statement follows. The proof for R_{p} is the same, although as always p must not divide l_{∗} so R_{p} is a ring.

Note that one can check if a ∈ R_{p} is invertible, and if so, compute its inverse, with linear algebra over ℤ_{(p)} or over its field of fractions ℚ. The matrix of the system to be solved is the matrix of ψ. The same also holds for \overline{a} ∈ \overline{R}, whenever it is invertible, its inverse can be computed with linear algebra over 𝔽_{p}. But instead of solving linear equations, we will use the extended Euclidean algorithm to calculate inverses in \overline{R}. However, this can increase the number of fail primes because the calculation can fail even if \overline{a} is invertible. This is not a serious problem because the number of fail primes will still be finite (see section 1.4).

In the following, let R be a commutative ring with identity 1 ≠ 0. For a univariate polynomial f ∈ R[x] define monic(f) as follows: If f = 0 then monic(f) = 0. If f ≠ 0 and if the leading coefficient lc(f) ∈ R of f is a unit, then define monic(f) = lc(f)^{−1}f. If f ≠ 0 and lc(f) is not a unit then define monic(f)=“failed”.

If f_{1}, f_{2} ∈ R[x] then the monic gcd is defined as a polynomial g ∈ R[x] such that g = monic(g) and for every polynomial h one has: h | f_{1} and h | f_{2} if and only if h | g. It is easy to show that if a monic gcd of f_{1}, f_{2} exists, then it is unique. The well-known Euclidean algorithm over R works as follows.

Euclidean algorithm.

Input: a list (f_{1}, f_{2}) of two univariate polynomials with coefficients in R. Output: Either a message “failed” or the monic gcd.

1. Set r_{1} = f_{1}, r_{2} = f_{2}, i = 2.

2. If r_{2} = 0 then set r_{1} = monic(r_{1}). If r_{1} = “failed” then return “failed”.

3. If r_{i} = 0 then return r_{i−1}.

4. Set r_{i} = monic(r_{i}). If r_{i} = “failed” then return “failed”.

5. Set r_{i+1} to be the remainder of r_{i−1} divided by r_{i}.

6. Set i = i + 1 and go back to Step 3.

Remark on a shortcut: Suppose that r_{i} in step 3 is a nonzero constant. Some implementations of the Euclidean algorithm over a field will then take a shortcut: stop the computation, the output is 1. Over a ring we should not use this shortcut because that would invalidate lemma 3 below. This plays a role because our algorithm will not test if p divides the discriminant. We may only use the shortcut if r_{i} is a unit. For r_{i} ∈ \overline{R} we can test that efficiently by computing N(r_{i}) mod p (see lemmas 1,2).

Denote GCD_{R}(f_{1}, f_{2}) as the output of this algorithm. If GCD_{R}(f_{1}, f_{2}) ≠ “failed” then the sequence of polynomials r_{1}, …, r_{m} with r_{m−1} ≠ 0, r_{m} = 0, is called the monic polynomial remainder sequence of f_{1}, f_{2}.

Lemma 3. If g = GCD_{R}(f_{1}, f_{2}) and g ≠ “failed” then the ideal (r_{i−1}, r_{i}) = R[x]r_{i−1} + R[x]r_{i} remains the same during each step. In particular (f_{1}, f_{2}) = (g) which implies:

1. There exist s, t ∈ R[x] such that g = sf_{1} + tf_{2}.

2. f_{1} and f_{2} are divisible by g.

3. g is the monic gcd of f_{1} and f_{2}.

Proof: When we make r_{i} monic, we divide by a unit, which does not change the ideal. In step 6 we increase i so we must show that (r_{i−1}, r_{i}) = (r_{i}, r_{i+1}) which is clear because r_{i+1} is the remainder of r_{i−1} modulo r_{i}. Hence (f_{1}, f_{2}) = (r_{1}, r_{2}) = (r_{m−1}, r_{m}) = (g, 0) = (g). So g ∈ (f_{1}, f_{2}) which is part 1, f_{1}, f_{2} ∈ (g) which is part 2. Finally, every h that divides both f_{1} and f_{2} divides any element of (f_{1}, f_{2}) in particular it divides g. Since g is monic it satisfies precisely the definition of the monic gcd.

Remark: If GCD_{R}(f_{1}, f_{2}) ≠ “failed” then the extended Euclidean algorithm, which calculates s and t as well as g will not fail either.

Let d = GCD_{R}(f_{1}, f_{2}) be the output of the Euclidean algorithm. If all leading coefficients during the computation are units then the algorithm succeeds, the monic gcd exists and equals d = r_{m−1}. If there is no monic gcd in R[x] then d = “failed”. If a monic gcd g does exist then it is not necessarily true that the algorithm will find it; the output d is then either g or “failed”. A situation where the output is “failed” even when a monic gcd exists is given in the following lemma.

Lemma 4. Suppose p ∤ l_{∗} and f_{1}, f_{2} ∈ R_{p}[x]. Then f_{1}, f_{2} ∈ O_{p}[x]. Suppose a monic gcd g ∈ O_{p}[x] exists and that g ∉ R_{p}[x]. Then GCD_{O_{p}}(f_{1}, f_{2}) = “failed”.

Proof: If p ∤ l_{∗} then α_{1}, …, α_{n} ∈ O_{p}, hence R_{p} ⊆ O_{p} so f_{1}, f_{2} ∈ O_{p}[x]. Since GCD_{R_{p}}(f_{1}, f_{2}) = “failed”, when we run the Euclidean algorithm over R_{p} we will encounter a leading coefficient in R_{p} that is not a unit in R_{p}. But according to the remark after lemma 2, if a ∈ R_{p} is not a unit in R_{p} then it is also not a unit in O_{p} and hence the algorithm fails over O_{p} as well.

If the ring R in the Euclidean algorithm is a field L, then the output is never “failed”, so GCD_{L}(f_{1}, f_{2}) is always the monic gcd of f_{1}, f_{2} ∈ L[x].

Lemma 5. Suppose f_{1}, f_{2} ∈ L[x] and r_{1}, …, r_{m} ∈ L[x] is the monic polynomial remainder sequence. Let lc_{1}, …, lc_{m−1} in L be the leading coefficients that we divided by in steps 2 and 4. For all but finitely many primes the following holds:

1. f_{1}, f_{2} ∈ R_{p}[x], and lc_{1}, …, lc_{m−1} are units in R_{p}.

2. r_{1}, …, r_{m} ∈ R_{p}[x] and \overline{r_{1}}, …, \overline{r_{m}} is the monic polynomial remainder sequence of \overline{f_{1}}, \overline{f_{2}}.

3. p is a good prime which means: The monic gcd of \overline{f_{1}}, \overline{f_{2}} exists, will be found by the Euclidean algorithm, and equals \overline{g} where g ∈ L[x] is the monic gcd of f_{1}, f_{2}.

Proof: Part 1 holds for all primes that do not divide any of the following: l_{∗}, den(f_{1}), den(f_{2}), den(lc_{i}), den(1/lc_{i}) for i < m. Since these are finitely many integers, all nonzero, we see that part 1 holds for all but finitely many primes. The only divisions in the Euclidean algorithm are divisions by lc_{i}, so if the input is in R_{p}[x] and all lc_{i} are units in R_{p}, then all polynomials in the GCD_{L}(f_{1}, f_{2}) computation are in R_{p}[x]. Induction shows that \overline{r_{1}}, …, \overline{r_{m}} is precisely the monic polynomial remainder sequence of \overline{f_{1}}, \overline{f_{2}}, so part 2 follows from part 1. Part 3 follows from part 2.

Since we will only run the Euclidean algorithm in \overline{R}[x] for various primes p, and not in L[x], we do not know the values of lc_{i}. So the lemma does not tell us which primes are good, it only says that all but finitely many primes are good. We now investigate the relation between GCD\overline{R}(\overline{f_{1}}, \overline{f_{2}}) and GCD_{L}(f_{1}, f_{2}) when p is not an lc-bad prime.

Theorem 1. Let f_{1}, f_{2} ∈ L[x] and let g ∈ L[x] be the monic gcd. Assume p ∤ l_{∗}den(f_{1})den(f_{2}), f_{2} ≠ 0 and lc(f_{2}) ≢ 0 mod p, so p is not an lc-bad prime. Let d = GCD\overline{R}(\overline{f_{1}}, \overline{f_{2}}). If d ≠ “failed” then

deg(d) ≥ deg(g).

Furthermore, if deg(d) = deg(g) then g reduces properly mod p and d = \overline{g}.

Remark: The theorem says that if p is not lc-bad then p is either fail, unlucky, or good. This implies that if lcbad primes are avoided then the modular GCD algorithm is correct.

Proof: lc(f_{2}) ≢ 0 mod p, so if we assume d ≠ “failed” then lc(f_{2}) must be a unit mod p, see step 4 in the Euclidean algorithm. There exist (see lemma 3) s_{0}, t_{0} ∈ R_{p}[x] such that

\overline{s_{0}}\overline{f_{1}} + \overline{t_{0}}\overline{f_{2}} = d.

Now take a monic polynomial d_{0} ∈ R_{p}[x] such that d = \overline{d_{0}}. Then we have

s_{0}f_{1} + t_{0}f_{2} ≡ d_{0} mod p.

We will apply Hensel lifting to increase the modulus p to a higher power of p. Define (starting with i = 1)

h_{i} = (s_{i−1}f_{1} + t_{i−1}f_{2} − d_{i−1})/p^{i} ∈ R_{p}[x]

and let q_{i}, r_{i} ∈ R_{p}[x] be the quotient and remainder of h_{i} divided by d_{0} (this division works because d_{0} is monic). Then define

s̃_{i} = s_{i−1} − p^{i}q_{i}s_{0}, t̃_{i} = t_{i−1} − p^{i}q_{i}t_{0}, d_{i} = d_{i−1} + p^{i}r_{i}.

Then

s̃_{i}f_{1} + t̃_{i}f_{2} ≡ d_{i} mod p^{i+1}.

Now s̃_{i}, t̃_{i} can have higher degrees than s_{i−1}, t_{i−1}. To remedy this, do the following. For j ∈ {1, 2} denote f_{j,d} ∈ R_{p}[x] as a polynomial whose modular image equals \overline{f_{j}}/d. Take q_{i}s_{0} mod p, and divide it by \overline{f_{2,d}} ∈ \overline{R}[x]. This division works because the leading coefficient of \overline{f_{2,d}} is lc(f_{2}) mod p, which is invertible. Take q, r ∈ R_{p}[x] such that \overline{q}, \overline{r} are the quotient and remainder of this division. Take q, r in such a way that they have the same degree as \overline{q}, \overline{r}. Then define

s_{i} = s_{i−1} − p^{i}r, and t_{i} = t_{i−1} − p^{i}(q_{i}t_{0} + qf_{1,d}),

and we still have

s_{i}f_{1} + t_{i}f_{2} ≡ d_{i} mod p^{i+1}.

We can now increase i and do the next Hensel step, and continue in this way. Because deg(r) < deg(\overline{f_{2,d}}) and deg(r_{i}) < deg(d_{0}), the degrees of s_{i} and d_{i} will be bounded as i increases, and hence the degree of t_{i} mod p^{i+1} is bounded as well. So when i → ∞, the limit ŝ, t̂, d̂ of s_{i}, t_{i}, d_{i} exists in the ring R̂_{p}[x] defined below. Denote ℤ̂_{p} as the ring of p-adic integers. ℤ̂_{p} is the completion of ℤ_{(p)} with respect to the p-adic valuation norm. Let ℚ̂_{p} be the field of p-adic numbers, the field of fractions of ℤ̂_{p}. Denote L̂_{p} = R_{p} ⊗_{ℤ_{(p)}} ℚ̂_{p} = L ⊗_{ℚ} ℚ̂_{p}. This is in general not an integral domain because minimal polynomials can become reducible when one replaces ℚ by a larger field ℚ̂_{p}. Denote R̂_{p} = R_{p} ⊗_{ℤ_{(p)}} ℤ̂_{p}. Now R̂_{p} and L can be viewed as subrings of L̂_{p} and

R_{p} = R̂_{p} ⋂ L  (5)

After doing infinitely many Hensel steps we find ŝ, t̂, d̂ ∈ R̂_{p}[x] such that

ŝf_{1} + t̂f_{2} = d̂.

Now d̂ is monic and deg(d̂) = deg(d_{0}) = deg(d) because the p^{i}r_{i}, i = 1, 2, …, that we added to d_{0} have smaller degree than d_{0}. The polynomials f_{1}, f_{2} are elements of L[x]g ⊆ L̂_{p}[x]g. Hence ŝf_{1} +t̂f_{2}, which equals d̂, is a also an element of L̂_{p}[x]g. But d̂ ≠ 0 so

deg(d) = deg(d̂) ≥ deg(g).

If the degrees are the same then d̂ = g because g is the only monic element of L̂_{p}[x]g of that degree. Equation (5) then implies g ∈ R_{p}[x] (recall that d̂ ∈ R̂_{p}[x] and g ∈ L[x]). So g can be reduced mod p. Hence g reduces properly mod p because it is monic. The theorem now follows because d equals d̂ mod p, which equals g mod p.

3. IMPLEMENTATION

In this section we describe our implementation of the modular GCD algorithm for multivariate polynomials over L. There are several multivariate “modular” GCD algorithms over ℚ that one may consider extending to work over L. We have completed an implementation of Brown’s algorithm (see [3]) which uses rational reconstruction and trial division (see [8]), and have begun work on an implementation of Zippel’s algorithm (see [11]). We have encountered three bottlenecks on real problems, namely, (i) rational reconstruction, (ii) the trial divisions, and (iii) extensions of low degree. We will address (i) and (ii) in this paper. Problem (iii) is addressed in [8].

We first give details of the data structure that we use for multivariate polynomials over L. The data structure is designed to make the modular GCD algorithm fast. It supports n ≥ 0 extensions over ℚ and 𝔽_{p}. To fix notation, recall that L = ℚ(α_{1}, …, α_{n}) where α_{i} is algebraic over L_{i−1} = ℚ(α_{1}, …, α_{i−1}), and m_{i}(z_{i}) ∈ L_{i−1}[z_{i}] is the minimal polynomial for α_{i} over L_{i−1}. Let R = L[x_{1}, …, x_{k}]. Let f_{1}, f_{2} be non-zero polynomials in R and let g be their monic GCD.

3.1 A Data Structure for R

In [9], Stoutemyer asked the question “Which polynomial representation is best?” (for a general purpose computer algebra system). Based on his data, he concluded that the recursive dense representation was best overall, a conclusion that ran contrary to the general belief that one must use a sparse representation. In our context, we have additional reasons to choose this representation. Our input polynomials to the modular GCD algorithm are multivariate polynomials in x_{1}, …, x_{k} and z_{1}, …, z_{n}. Because the modular GCD algorithm is recursive and because the extensions must also be defined recursively if they are dependent, a recursive data structure will minimize data structure overhead. Because we compute the GCD modulo machine primes p_{1}, p_{2}, …, most of the work takes place in the last variable, i.e. in the ring 𝔽_{p}[z_{1}]. Since the bottom level of the recursive dense representation is a dense vector of machine integers, this yields the best representation for arithmetic in 𝔽_{p}[z_{1}]. We now describe the data structure <poly> using a BNF notation with some examples.

                                                          <poly> ::= POLYNOMIAL( <ring>, <data> )
                                                          <ring> ::= [ <char>, <vars>, <exts> ]
                                                          <char> ::= <nonnegative integer>
                                                          <data> ::= <rational number> | <immediate integer>
                                                                                       | vector(<data>)
                                                          <vars> ::= vector(<variables>)
                                                          <exts> ::= vector(<data>)

The characteristic of the ring is encoded by <char> and <exts> is a vector of the minimal polynomials. Thus the ring for the polynomial is encoded in the data structure. Since this information is identical for polynomials in the same ring it should be stored once so that the cost of storing the ring information is one word.

We impose the following restriction which is a key property of the recursive dense representation; a zero coefficient at any level in the data structure is represented by the immediate integer 0 (or nil pointer). This means that every algorithm must treat 0 as a special case. This exceptional case does not bother us greatly because in the implementation of most operations 0 is a special case anyway.

The bottom of the data structure is a word of storage which is either a pointer to a rational number or an immediate integer. In our Maple implementation, immediate integers are signed integers of 30 bits in length, hence, one bit is used to distinguish them from pointers. In the examples below, vectors are indicated by square brackets.

Example 1: The representation of the polynomial z^{4} − 10z^{2} + 1 in characteristic 0 and characteristic 3 is

                                                         POLYNOMIAL( [0,[z],[]], [1,0,-10,0,1] )
                                                         POLYNOMIAL( [3,[z],[]], [1,0,2,0,1] )

The empty vector [] indicates that there are no extensions and the data in both these examples is a vector of machine integers. Allowing one word as a header word for the POLY-NOMIAL structure and for each vector, the storage requirement for both polynomials is 16 words (count one word for POLYNOMIAL and each [ in the above). Since the ring information can be shared between polynomials over the same ring, a more accurate count is that 9 words are required. From now on we will not count the storage for the ring.

Example 2: The representation of the polynomial x^{3}−zx+ z^{2} in ℚ[z][x] and ℚ[z]/(z^{2} − 2) is

 POLYNOMIAL([0,[x,z],[]], [[0,0,1],[0,-1],0,[1]])
 POLYNOMIAL([0,[x,z],[[-2,0,1]]],[[2],[0,-1],0,[1]])

In the data structure, polynomials are reduced modulo the m_{i} on input. The storage requirement is 17 and 15 words respectively.

Example 3: The recursive dense data structure is not sparse, but neither is it truly dense. On sparse polynomials, the storage requirement is still very good. Consider the sparse polynomial 1 + 2x^{n} + 3y^{n} + 4z^{n} where n = 3. Our data structure for this polynomial is

 POLYNOMIAL(R,[[[1,0,0,4],0,0,[3]],0,0,[[2]]]);

This is 24 words (not counting the storage for the ring). In general it is 15 + 3n words. One of the main sparse representations for polynomials that is used in AXIOM is a linked list of pairs where each pair is a pointer to a coefficient and a pointer to a monomial where the monomial x^{i}y^{j}z^{k} would be stored as an exponent vector [i, j, k]. Thus each non-zero term of the polynomial requires 2 + 2 + 4 = 8 words of storage. On our example this would be 35 words, allowing 3 words for the top level of the data structure. Of course, this is not truly a sparse data structure because the monomial representation is not sparse. Nevertheless, on this example, the recursive dense representation uses less storage for n ≤ 6.

Example 4: Multiple extensions are handled in the ob-

vious way. The polynomial x + \sqrt{1/5} y + \sqrt{1 + \sqrt{1/5}} is represented by

 POLYNOMIAL([0, [x,y,z2,z1], [[[-1,-1],0,[1]],
  [-1, 0, 5]]], [[[0,[1]], [[0,1]]], [[[1]]]])

where the two minimal polynomials m_{1}(z_{1}) and m_{2}(z_{2}) have been replaced by m̌_{1} and m̌_{2}. We remark that although this makes the test for whether a prime p in the modular GCD algorithm divides den(m_{i}) easy, and makes reduction of the minimal polynomials modulo p easy, after having reduced the minimal polynomials modulo p, one should make them monic over 𝔽_{p} so that we do not repeatedly invert their leading coefficients.

3.2 Rational Reconstruction

If one naively applies rational reconstruction after computing the GCD modulo each prime p, the cost of the rational reconstruction may become the bottleneck asymptotically as well as in practice. Suppose g = x + a/b where a and b are integers of m digits in length. Then we need O(m) primes to reconstruct g. The total cost of Chinese remaindering will be O(m^{2}) but the total cost of rational reconstruction will be O(m^{3}). The reason is that Chinese remaindering can be done “incrementally” after each prime in O(m) time but, as far as we know, rational reconstruction cannot. We resolve this problem by attempting rational reconstruction after 1, 2, 3, 5, 8, 13, … primes. This ensures that the cost of rational reconstruction is also O(m^{2}).

We now state the modular GCD algorithm. As a preprocessing step, we compute f̌_{1} and f̌_{2}, that is, we cancel any rational scalar from the input polynomials before proceeding. We do not compute f̃_{1} or f̃_{2} which can cause a blowup. Let ic(f) denote the integer content of a polynomial. This is the rational constant c such that f/c = f̌, the semi-associate of f. For example, if f = 10/3x^{2} − 15 then ic(f) = 5/3 and f̌ = 2x^{2} − 9.

Modular GCD algorithm. Input: Non-zero f_{1}, f_{2} ∈ L[x]. Output: g, the monic GCD.

1. Set n = 0, f_{1} = f_{1}/ic(f_{1}) and f_{2} = f_{2}/ic(f_{2}).

2. Take a new prime p that is not lc-bad.

3. Let d be the result of the Euclidean algorithm mod p. If d = “failed” then go back to Step 2.

4. If d = 1 then return 1.

5. If n = 0 or deg(d) < deg(G) then set G = d, m = p, n = 1 and go to Step 9.

6. If deg(d) > deg(G) then go back to Step 2.

7. Let G be the result of Chinese remaindering on G mod m and d mod p. Set m = mp, n = n + 1.

8. If n is not a Fibonacci number go back to Step 2.

9. Apply rational reconstruction to obtain h ∈ L[x] from G mod m. If this fails, go back to Step 2.

10. Apply the Euclidean algorithm mod one more prime p to obtain d and test if h mod p = d. If not, go back to Step 4.

11. Trial division. If h|f_{1} and h|f_{2} then return h, otherwise go back to Step 4.

The algorithm will terminate as soon as we have enough good primes to reconstruct the coefficients of g. We will (i) state the expected running time in terms of f̌_{1} and f̌_{2}, (ii) not include the cost of the trial divisions, and (iii) assume classical (quadratic) algorithms for integer and polynomial arithmetic.

Let m be the number of good primes needed to reconstruct g. The algorithm as stated has been designed so that m is proportional to the length of largest integer in the coefficients of g. We will assume that the probability that a prime is good is high, so that m is close to the actual number of primes that were used. This assumption is true in practice. However, for theoretical completeness of the complexity estimate, we would need to determine some B = B(f_{1}, f_{2}) such that if p > B then the probability that p is good is greater than some constant. We did not determine such B because this issue would not have consequences for the algorithm in practice (one hardly ever encounters primes that are not good).

The other quantities appearing in the running time are (i) D, the degree of the number field L, (ii) N, the degree of the largest input polynomial, (iii) n the degree of the GCD, and (iv) M the size of the largest integer coefficient in f̌_{1} and f̌_{2}. The average running time of the modular GCD algorithm is O(mMND + mN^{2}D^{2} + m^{2}nD)) where the three contributions are for reducing f̌_{1} and f̌_{2} modulo m primes, applying the Euclidean algorithm m times, and the reconstruction time, respectively. If we consider the case where M = 2m and N = 2n corresponding to a GCD problem where the cofactors are the same size as the GCD, then this simplifies to O(MND(M + ND)).

3.3 Trial Division

Another bottleneck of the modular GCD algorithm is the trial divisions. If h is the result of rational reconstruction then we must check that h|f_{1} and h|f_{2} to show that h = g. Because these trial divisions can be expensive, we have considered abandoning trial divisions altogether in favor of a probabilistic result, that is, check that result of rational reconstruction agrees, say, with the GCD modulo five additional primes instead of one. However, in many applications where one computes GCDs, for example, normalizing a rational function, one wants to compute also the cofactors f_{1}/g and f_{2}/g, hence, the divisions cannot be avoided.

We can use either classical division or a modular division algorithm. If g is small in size compared with f_{1} and f_{2} then classical division is asymptotically faster. On the other hand, if g is of similar size to its cofactors f_{1}/g and f_{2}/g then a modular division algorithm will be asymptotically faster. When dividing f_{1} and f_{2} by h over L using the classical division algorithm, a significant improvement (we saw a speedup of a factor of 10 on one large example) can be obtained if one avoids fractions as much as possible. Notice that the leading coefficient of ȟ in the modular GCD algorithm is an integer. If also l_{i} = den(m_{i}) = 1, which is often the case, then the entire division algorithm can be completed using only integer arithmetic. If l_{i} ≠ 1 for some i then the division algorithm can still be modified to avoid fractions. We will describe how to do this for univariate polynomials with one field extension with minimal polynomial M.

Algorithm Fraction Free Trial Division. Input: A, B ∈ ℚ[x, z], M ∈ ℤ[z] : B ≠ 0, lc_{x} B ∈ ℚ, and deg M ≥ 1. Output: Q = A/B mod M if B|A mod M; “failed” otherwise.

Set m = deg_{x} A, n = deg_{x} B and d = deg_{z} M. Set i_{a} = ic(A) and a = A/i_{a}. Set i_{b} = ic(B) and b = B/i_{b}. Set l_{b} = lc_{x} b and l_{m} = lc_{z} M. Set s = 1, r = a, and q = 0. While r ≠ 0 and m ≥ n do

Set l_{r} = lc_{x} r. Note that l_{r} ∈ ℤ[z]. Set g = GCD(ic(l_{r}), l_{b}) and l_{r} = l_{r}/g.

Set s = (l_{b}/g) × s.

Set t = l_{r} × x^{m−n} and q = q + t/s.

Set r = (l_{b}/g) × r − t × b.
Set p = 1.

While r ≠ 0 and deg_{z} r ≥ d do Set l_{r} = lc_{z} r. Note that l_{r} ∈ ℤ[x]. Set g = GCD(ic(l_{r}), l_{m}) and l_{r} = l_{r}/g. Set t = l_{r}z^{deg_{z} r−m} and p = p × (l_{m}/g). Set r = (l_{m}/g) × r − t × M.

Set s = s × p.
Set m = deg_{x} r.

If r ≠ 0 then output “failed”. Set Q = (i_{a}/i_{b}) × q and output Q.

The algorithm first makes the inputs A and B primitive over ℤ. We claim that each time round the outer loop r and q satisfy a = bq + cr for some scalar c ∈ ℚ and r has integer coefficients. The outer loop reduces the degree of the remainder r in x. In the outer loop we multiply r by the smallest possible integer so that lc_{x} r, a polynomial in ℤ[z], will be divisible by lc_{x} b. The inner loop then reduces the remainder r modulo M. In the inner loop we multiply r by the smallest integer so that lc_{z} r, a polynomial in ℤ[x], will be divisible by lc_{z} M. The integers s and p are multipliers. They keep track of the integer factors of lc_{x} b and lc_{z} M, respectively, that r was multiplied by so that the quotient Q may be correctly computed from q.

Acknowledgement

We acknowledge the help of Craig Pastro with parts of the implementation, and the referees for valuable comments.

4. REFERENCES

[1] J. A. Abbott, R. J. Bradford, J. H. Davenport, The Bath Algebraic Number Package, Proceedings of SYMSAC ’86, ACM press (1986), pp. 250–253.

[2] R. J. Bradford, Some Results on the Defect, Proceedings of ISSAC ’89, ACM press (1989), pp. 129–135.

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

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

[5] E. Hecke, Lectures on the Theory of Algebraic Numbers, Springer Graduate Texts in Mathematics 77, (1981).

[6] M. van Hoeij, M. Monagan, A Modular GCD algorithm over Number Fields presented with Multiple Extensions, FSU preprint 02-03, (2002).

[7] L. Langemyr, S. McCallum, The Computation of Polynomial GCD’s over an Algebraic Number Field, J. Symbolic Computation 8 (1989), pp. 429–448.

[8] M. B. Monagan, A. D. Wittkopf, On the Design and Implementation of Brown’s Algorithm over the Integers and Number Fields, Proceedings of ISSAC ’2000 (2000), ACM Press, pp. 225–233.

[9] D. Stoutemyer, Which Polynomial Representation is Best?, Proceedings of the 1984 Macsyma User’s Conference, 1984.

[10] P. Wang, M. J. T. Guy, J. H. Davenport, p-adic Reconstruction of Rational Numbers, in SIGSAM Bulletin, 16, No 2 (1982).

[11] R. Zippel, Probabilistic algorithms for sparse polynomials, Proceedings of EUROSAM ’79, Springer-Verlag LNCS, 2 (1979), pp. 216–226.