A NEW HENSELIAN CONSTRUCTION AND ITS APPLICATION TO POLYNOMIAL GCDS OVER DIRECT

PRODUCTS OF FIELDS

François Boulier  Marc Moreno Maza  Cosmin Oancea

Abstract

Let K be a field of multivariate rational functions over an infinite field. Let L be a direct product of fields extending K by a tower of algebraic extensions. We propose a modular algorithm for computing polynomial gcds over L based on a Hensel lifting strategy. The rational reconstruction is avoided by guessing the denominator of the gcd before the lifting step. The algorithm is probabilistic but succeeds with probability one. Our implementation shows a significant improvement with respect to other methods.

Introduction

Let k be an infinite field and let R be an integral domain which is either the ring ℤ of integers or a ring k[t_{1}, …, t_{n}] of multivariate polynomials over k. Let L_{0} = K be the field of fractions of R. Let x_{1} < ⋯ < x_{p} be ordered variables and let m_{1}, …, m_{p} be nonconstant polynomials in R[x_{1}, …, x_{p}]. Let i be any integer in the range 1 ⋯ p. We assume that the greatest variable occurring in m_{i} is x_{i} and that the leading coefficient lc(m_{i}, x_{i}) of m_{i} w.r.t. x_{i} lies in R. Let L_{i} be the simple algebraic extension defined by m_{i} over L_{i−1}. In other words, the ring L_{i} is the residue class ring L_{i} = L_{i−1}[x_{i}]/⟨m_{i}⟩. We assume that the polynomial m_{i} generates a radical ideal in L_{i−1}[x_{i}]. We define L = L_{p}. Let y be an additional variable and let f_{1}, f_{2} be polynomials in L[y] with regular (i.e. not zero-divisors) leading coefficients. It is known that every ring L_{i} is a direct product of fields. Therefore, it makes sense to consider the gcds of f_{1}, f_{2} in L[y] (although a monic gcd of f_{1}, f_{2} may not always exist even if f_{1}, f_{2} are both monic). See [MR95] for more detail.

Langemyr and McCallum provides in [LM89] an algorithm for computing gcds in L[y] for the case where L is a field, R = ℤ and p = 1. This algorithm is improved by Encarnación in [Enc94]. In [MR95] the authors proposed an algorithm for the general case where L is not necessarily a field, p may be greater than 1 and R is any integral domain. This non-modular algorithm combines subresultant theory and the following remark.

In order to compute a gcd of f_{1}, f_{2} in L[y] one can apply the D5 principle [DDD85] to the Euclidean algorithm in L[y]. That is, running the Euclidean algorithm in L[y] as if L was a field but checking that the leading coefficient of every intermediate remainder is regular. When a zero-divisor is met, then a decomposition L_{1} ⊕ L_{2} of L has been found and the computations should split into two branches: one in L_{1}[y] and one L_{2}[y]. We denote by g_{monic} the polynomial of L[y] obtained by this process. Again, the polynomial g_{monic} is not necessarily monic but for every field F factor of L the projection of g_{monic} in F [y] is monic.

Let us regard now the set M = {m_{1}, …, m_{p}} as a regular chain in R[x_{1}, …, x_{p}]. We define U = R[x_{1}, …, x_{p}]/Sat(M) where Sat(M) is the saturated ideal of M. See [ALM99] for the notions related to regular chains. One can run the subresultant algorithm in U[y] as if U was an integral domain but checking that the leading coefficient of every intermediate pseudo-remainder is regular. If no zero-divisor is met during this process, then we obtain a polynomial g_{subres} ∈ U[y] such that g_{monic} and g_{subres} are equal in L[y] up to a multiplicative factor which is a regular element of L. This alternative approach is introduced in [Mor99].

Monagan and van Hoeij have proposed two modular algorithms for the case where L is a field. The first one in [HM02] assumes that R = ℤ holds. The second one in [HM04] applies to the case where R = k[t_{1}, …, t_{n}] but assumes p = 1. Both algorithms are based on the Chinese remaindering algorithm and the rational reconstruction.

In this paper, we propose an algorithm which applies to case where L is not necessarily a field, p may be greater than 1 and R = k[t_{1}, …, t_{n}] (as those of [MR95] and [Mor99]) but which is a modular algorithm (as those of [HM02] and [HM04]). A preliminary implementation shows that our new algorithm can process problems that were out of the scope of the algorithms of [MR95] and [Mor99]. We rely on a new Henselian construction presented in Section 1. Our approach avoids the rational reconstruction by guessing the denominator of the gcd, as discussed in Section 2. In addition, our algorithm is probabilistic (in the sense that it may run forever) but succeeds (in a finite number of steps) with probability one. The proof (not reported here) of this latter fact relies on subresultant theory, in particular on the results of [El 03], and a key notion that we introduce now.

Checking the regularity of an element in U is made by running the subresultant algorithm in (R[x_{1}, …, x_{ℓ−1}]/Sat(m_{1}, …, m_{ℓ−1}))[x_{ℓ}] for some ℓ ∈ {1, …, p}. This leads to the following concept. We call the history of g_{subres} the sequence denoted by history(g_{subres}) of the degrees of all intermediate pseudo-remainders appearing when computing g_{subres} (including those pseudo-remainders involved when testing the regularity of an element). Similarly we define history(g_{monic}) and observe that history(g_{subres}) and history(g_{monic}) coincide.

1 A new Henselian construction

Let M_{0} be a maximal ideal of R of the form ⟨t_{1} −t_{1,0}, …, t_{n} −t_{n,0}⟩. We denote by R_{M_{0}} the residue class field R/M_{0} and by M_{0} = {m_{1,0}, …, m_{p,0}} the image of the regular chain M in R_{M_{0}}[x_{1}, …, x_{p}]. We assume that no lc(m_{i}, x_{i}) is null modulo M_{0} such that M_{0} is a regular chain in (R/M_{0})[x_{1}, …, x_{p}]. Observe that M_{0} is also a lexicographical Gröbner basis for the monomial ordering induced by x_{1} < ⋯ < x_{p} and that Sat(M_{0}) is just the ideal generated by M_{0}. Finally, we denote by U_{M_{0}} the residue class ring R_{M_{0}}[x_{1}, …, x_{p}]/Sat(M_{0}).

Theorem 1.1. Assume that Sat(M_{0}) is radical (as it is the case for Sat(M)) such that U_{M_{0}} is a direct product of fields. Let f ∈ R[x_{1}, …, x_{p}][y] and f_{0} be its image in R_{M_{0}}[x_{1}, …, x_{p}][y]. Let c_{g}, c_{h} ∈ R[x_{1}, …, x_{p}] be units of L such that lc(f, y) = c_{g} c_{h}. Let d_{1}, …, d_{n} be positive integers. Let g_{0} and h_{0} be polynomials in R_{M_{0}}[x_{1}, …, x_{p}][y] with respective degrees d_{g} and d_{h} and such that g_{0} and h_{0} are relatively prime in U_{M_{0}}[y], their product g_{0}h_{0} is equal to f_{0} in U_{M_{0}}[y], and their leading coefficients lc(f_{0}, y) and lc(g_{0}, y) are units in U_{M_{0}} equal respectively to c_{g} and c_{h} in U_{M_{0}}. Then, there exists at most one couple (g, h) of polynomials in L[y] such that

(i) f equals g h in L[y],

(ii) g_{0} and h_{0} are the respective images of g and h in U_{M_{0}}[y],

(iii) lc(g) = c_{g} and lc(h) = c_{h},

(iv) For every i = 1, …, n we have max(deg(g, t_{i}), deg(h, t_{i})) ≤ d_{i}.

In addition, if it exists, such couple (g, h) can be constructed effectively in R[x_{1}, …, x_{p}][y].

For two polynomials g, h ∈ R[x_{1}, …, x_{p}][y] the equation f = g h in L[y] means f equals g h modulo Sat(M). That is, there exist polynomials r_{1}, …, r_{p} ∈ R[x_{1}, …, x_{p}][y] and a product λ of (powers of) the lc(m_{i}, x_{i}) such that

λ (f − g h) = r_{1}m_{1} + ⋯ + r_{p}m_{p}.  (1)

We assume that we have an initial guess for λ. If the construction to be described below does not succeed in computing g and h this implies either that no couple (g, h) of polynomials in L[y] satisfying (i), (ii), (iii), (iv) exists or that λ was not adequate. However, the construction will show how λ must be corrected for a next try. This feature validates the whole process. Now, let λ_{0} be the image of λ modulo M_{0}. Because of the properties of M_{0} described earlier, one can compute r_{1,0}, …, r_{p,0} ∈ R_{M_{0}}[x_{1}, …, x_{p}][y] such that

λ_{0} (f_{0} − g_{0} h_{0}) = r_{1,0}m_{1,0} + ⋯ + r_{p,0}m_{p,0}.  (2)

The fact that g_{0} and h_{0} are relatively prime in U_{M_{0}}[y] means that there exist polynomials σ_{0}, π_{0} ∈ R_{M_{0}}[x_{1}, …, x_{p}][y] and polynomials ρ_{0}, τ_{0}, θ_{1,0}, …, θ_{p,0} ∈ R_{M_{0}}[x_{1}, …, x_{p}] such that

σ_{0}g_{0} + π_{0}h_{0} = ρ_{0} and ρ_{0}τ_{0} = 1 + θ_{1,0}m_{1,0} + ⋯ + θ_{p,0}m_{p,0}.  (3)

Note that all these polynomials can be effectively computed, for instance by the algorithms of [MR95]. For simplicity, we sketch the remaining of the construction of g and h for the case n = 1. The general case can be obtained by the techniques discussed in Section 6.7 of [GCL92]. We define q = t_{1} − t_{1,0}. Then we set g^{(1)} = (c_{g} − lc(g_{0}, y))y^{d_{g}} + g_{0}, h^{(1)} = (c_{h} − lc(h_{0}, y))y^{d_{h}} + h_{0}, r_{1}^{(1)} = r_{1,0}, …, r_{p}^{(1)} = r_{p,0}. For 1 ≤ k ≤ d_{1} + 1, we shall construct g^{(k)}, h^{(k)}, r_{1}^{(k)}, …, r_{p}^{(k)} ∈ R[x_{1}, …, x_{p}][y] such that g^{(k)}, h^{(k)}, r_{1}^{(k)}, …, r_{p}^{(k)} are equal modulo q^{k} to the polynomials g, h, r_{1}, …, r_{p} to be computed, and g^{(k)}, h^{(k)}, r_{1}^{(k)}, …, r_{p}^{(k)} are equal modulo q to g^{(1)}, h^{(1)}, r_{1}^{(1)}, …, r_{p}^{(1)}. To this end, we define

F (u, v, w_{1}, …, w_{p}) = λ(f − u v) − w_{1}m_{1} − ⋯ − w_{p}m_{p}

where u, v, w_{1}, …, w_{p} are any polynomials in R[x_{1}, …, x_{p}][y]. Proceeding by induction on k ≥ 1 and using the Relations (1) and (2) one can easily show that

F (g^{(k)}, h^{(k)}, r_{1}^{(k)}, …, r_{p}^{(k)}) ≡ F (g, h, r_{1}, …, r_{p}) ≡ 0 mod q^{k}.  (4)

Expressing F (g^{(k+1)}, h^{(k+1)}, r_{1}^{(k+1)}, …, r_{p}^{(k+1)}) with the Taylor expansion of F around g^{(k)}, h^{(k)}, r_{1}^{(k)}, …, r_{p}^{(k)} and using Relation (4) leads to

ϵ_{k} = λ_{0}h_{0}g_{k} + λ_{0}g_{0}h_{k} + r_{1,k}m_{1,0} + … + r_{p,k}m_{p,0}  (5)

where g_{k}, h_{k}, r_{1,k}, …, r_{p,k} are the coefficients of q^{k} in g^{(k)}, h^{(k)}, r_{1}^{(k)}, …, r_{p}^{(k)} and where

ϵ_{k} = eval(\frac{E_{k}}{q^{k}} , t_{1} = t_{1,0}) and E_{k} = F (g^{(k)}, h^{(k)}, r_{1}^{(k)}, …, r_{p}^{(k)}).  (6)

This leads to the following Diophantine equation over the direct product of fields U_{M_{0}}

ϵ_{k} = λ_{0}h_{0}g_{k} + λ_{0}g_{0}h_{k}.  (7)

Since the leading coefficients of g_{k} and h_{k} are set to be c_{g} and c_{h} we can assume that we are looking for (g_{k}, h_{k}) such that deg(g_{k}) < d_{g} = deg(g_{0}) and deg(h_{k}) < d_{h} = deg(h_{0}). Therefore, it is easy to check that Equation (7) admits a unique solution which can be computed by means of Relation (3). (Observe that we have the equality ϵ_{k} = ρ_{0}τ_{0}ϵ_{k} in U_{M_{0}}[y].) Finally, the coefficients r_{1,k}, …, r_{p,k} are obtained from Equation (5).

If some E_{k} becomes null, then the construction succeeds by returning (g, h) = (g^{(k)}, h^{(k)}). If no E_{k} is null, this may be because λ is not big enough. To deal with this difficulty, the error E_{k} should be computed by pseudo-dividing f − g^{(k)}h^{(k)} w.r.t. M and comparing the actual λ used in this pseudo-division with the guessed λ. If the actual λ does not divide the guessed λ, then the construction should be started again with the correct λ.

2 Polynomial gcds over direct products of fields

Assume that the polynomials f_{1} and f_{2} admit a monic g_{monic} gcd in L[y]. The case where this assumption does not hold can be handled by splitting L. We describe now an algorithm which computes α ∈ R[x_{1}, …, x_{p}] and g, h ∈ R[x_{1}, …, x_{p}][y] such that

(1) α is a regular element in U,

(2) g is equal to g_{monic} in L[y] up to a multiplicative constant which is a unit of L,

(3) αf = gh holds in U[y] and therefore in L[y].

The key idea is to compute the denominator α before computing g. The computation of α is based on the following remark. Let r_{1} = f_{1}, r_{2} = f_{2}, …, r_{i} = q_{i+1} \overline{r_{i+1}} + r_{i+2}, …, r_{ℓ} be the Euclidean remainder sequence of f_{1} and f_{2} in L[y] where r_{i} = u_{i}\overline{r_{i}} and u_{i} = lc(r_{i}). With e_{i} = d_{i} − d_{i−1} + 1 the following equalities are fraction-free: u_{i+1}^{e_{i+1}}r_{i} = u_{i+1}^{e_{i+1}−1}q_{i+1}r_{i+1} + u_{i+1}^{e_{i+1}}r_{i+2}. Therefore, we can choose α = u_{2}^{e_{2}}u_{3}^{e_{3}} ⋯ u_{ℓ}^{e_{ℓ}} mod Sat(T). In practice, the degrees of α w.r.t. t_{1}, …, t_{n} are smaller than those of g. This is why reconstructing α before g is an efficient approach. We sketch now the main steps of our algorithm.

Step 1 Get an initial guess for the bounds d_{1}, …, d_{n} of Theorem 1.1. Get an initial guess N for the number of maximal ideals M_{0} that are needed in order to compute α by interpolation. Lastly, initialize a table T that we shall use for storing the history of a polynomial gcd of f_{1} and f_{2} in U_{M_{0}}[y], for every maximal ideal M_{0} that we use.

Step 2 Choose a maximal M_{0} such that the hypothesis of Section 1 hold.

Step 3 Compute α_{0} ∈ R_{M_{0}}[x_{1}, …, x_{p}] regular element in U_{M_{0}} and g_{0}, h_{0} ∈ R_{M_{0}}[x_{1}, …, x_{p}][y] such that g_{0} is a gcd of f_{1}, f_{2} in U_{M_{0}}[y] and α_{0}f = g_{0}h_{0} holds in U_{M_{0}}[y]. The computation of α_{0} is based on the above remark for α.

Step 4 If the history of g_{0} matches at least N − 1 other gcd computation histories in T then interpolate their α_{0}’s into an α and go to Step 5 otherwise go to Step 2.

Step 5 Apply the construction of Section 1. If it succeeds then return (α, g, h) otherwise revise (i.e. increase) the estimation of d_{1}, …, d_{n}, N and go to Step 2.

Definition 2.1. We say that M_{0} is a

lc-bad prime if one of the leading coefficients lc(m_{1}, x_{1}), …, lc(m_{p}, x_{p}) is zero modulo M_{0} or U_{M_{0}} is not a product of fields, or one of the leading coefficients lc(f_{1}, y), lc(f_{1}, y) is not regular modulo M_{0},

fail prime if M_{0} is not lc-bad but a subresultant gcd could not be computed modulo M_{0} because a zero-divisor was found,

unlucky prime if M_{0} is not a lc-bad or fail prime but the history of the subresultant gcd of f_{1}, f_{2} in U_{M_{0}}[y] is different from the history of the subresultant gcd of f_{1}, f_{2} in U[y] (the one we called g_{subres} above),

good prime if M_{0} is not a lc-bad prime, nor a fail prime, nor an unlucky prime.

Theorem 2.2. Assume M_{0} is not a lc-bad prime. Then, M_{0} is a good prime with a probability one and the subresultant gcd of f_{1}, f_{2} computed for this specialization is the image of g_{subres}. In particular, for the case R = k[t] (i.e. n = 1) there is only a finite number of specializations that are not good primes.

Theorem 2.3. Assume M_{0} is a good prime. Then, the polynomial g_{0} (computed in U_{M_{0}}[y]) divided by lc(g_{0}, y) is the image modulo M_{0} of the polynomial g_{monic} (computed in L[y]).

References

[ALM99] P. Aubry, D. Lazard, and M. Moreno Maza. On the theories of triangular sets. J. Symb. Comp., 28(1-2):105–124, 1999.

[DDD85] J. Della Dora, C. Discrescenzo, and D. Duval. About a new method method for computing in algebraic number fields. In Proc. EUROCAL 85 Vol. 2, volume 204 of Lect. Notes in Comp. Sci., pages 289–290. Springer-Verlag, 1985.

[El 03] M. El Kahoui. An elementary approach to subresultants theory. J. Symb. Comp., 35:281–292, 2003.

[Enc94] M. J. Encarnación. On a modular algorithm for computing gcds of polynomials over algebraic number fields. In Proc. ISSAC’94, pages 58–65. ACM Press, 1994.

[GCL92] K.O. Geddes, S.R. Czapor, and G. Labahn. Algorithms for Computer Algebra. Kluver Acadenic Publishers, 1992.

[HM02] M. van Hoeij and M. Monagan. A modular gcd algorithm over number fields presented with multiple extensions. In Teo Mora, editor, Proc. ISSAC 2002, pages 109–116. ACM Press, July 2002.

[HM04] M. van Hoeij and M. Monagan. Algorithm for polynomial gcd computation over algebraic function fields. (accepted for ISSAC 2004), 2004.

[LM89] L. Langemyr and S. McCallum. The computation of polynomial greatest common divisors over an algebraic number field. J. Symbolic Comp, 8:429–448, 1989.

[Mor99] M. Moreno Maza. On triangular decompositions of algebraic varieties. Technical Report TR 4/99, NAG Ltd, Oxford, UK, 1999. Presented at the MEGA-2000 Conference, Bath, England.

[MR95] M. Moreno Maza and R. Rioboo. Polynomial gcd computations over towers of algebraic extensions. In Proc. AAECC-11, pages 365–382. Springer, 1995.

François Boulier  Marc Moreno Maza  Cosmin Oancea
Université Lille 1  Université Lille 1  University of Western Ontario
59655 Villeneuve d’Ascq  &  London N6A 5B7
France  University of Western Ontario  Canada
boulier@lifl.fr  moreno@orcca.on.ca  oancea@orcca.on.ca