A Modular Resultant Algorithm for Number Fields ^{∗}

M.B. Monagan, S.C. Tse and A.D. Wittkopf

Department of Mathematics, Simon Fraser University, Burnaby, British Columbia, V5A 1S6 Canada.

{monagan,stephent,wittkopf}@cecm.sfu.ca

ABSTRACT

We present a modular algorithm for efficiently computing the Sylvester resultant of two multivariate polynomials over an algebraic number field. The algorithm is implemented in Maple using a new data structure which is designed to support modular algorithms. The algorithm and the implementation support multiple field extensions.

1. INTRODUCTION

Let R be a commutative ring and let f(x) = a_{n}x^{n} + ... + a_{1}x + a_{0} and g(x) = b_{m}x^{m} + ... + b_{1}x + b_{0} be polynomials in R[x]. The resultant of f(x) and g(x), denoted res_{x}(f, g), is the determinant of Sylvester’s matrix. Assuming division in R, the most direct approach to compute res_{x}(f, g) is to compute the determinant of Sylvester’s matrix using Gaussian elimination. If n ≥ m, this requires O(mn^{2}) arithmetic operations in R. Faster, the Euclidean algorithm can be modified to compute the resultant in O(mn) operations.

If R is a polynomial ring in one or more variables, as is the case with most applications of resultants in computer algebra, these algorithms lead to computations in the fraction field of R where the cost of polynomial gcds is relatively expensive. The gcds can be avoided in both cases by using fraction-free algorithms which make only exact polynomial divisions in R. The best fraction-free modification of the Euclidean algorithm for computing res_{x}(f, g) is the subresultant algorithm [5]. However, in general, in all of these algorithms, an expression swell occurs in the intermediate steps and it is more efficient to apply a modular approach. The first person to do this was Collins in [1] for R = ℤ[y_{1}, y_{2}, ..., y_{k}]. To compute the resultant r(y_{1}) = res_{x}(f, g) for R = ℤ[y_{1}] using Collin’s algorithm, or the modular resultant algorithm as it is often called, one computes the resultant modulo enough primes so that the answer over ℤ can be obtained by applying the Chinese remainder theorem.

∗This work is supported by NSERC of Canada and the MI-TACS NCE of Canada.

For each prime p_{i}, one computes the resultant in ℤ_{p_{i}}[y_{1}] at enough points y_{1} = β_{j} so that the resultant can be obtained by polynomial interpolation. This reduces the computation to many resultants in one variable over a finite field ℤ_{p} where no expression swell occurs. Collins proves that provided the primes and evaluation points are chosen so that the leading coefficients a_{n} and b_{m} do not vanish, the resultant can be recovered when sufficiently many primes and evaluation points are used. A bound on the largest integer coefficient and on the degree of any y_{i} appearing in the resultant can be obtained from inspecting Sylvester’s matrix. We refer to [3] for a better coefficient bound.

In this paper we consider the case R = K[y_{1}, y_{2}, ..., y_{k}] where K = ℚ(α_{1}, α_{2}, ..., α_{t}) and the α_{i} are algebraic over ℚ, that is, R is a polynomial ring over an algebraic number field. One application is Trager’s algorithm [9] for factoring polynomials in K[x] where resultants are used to compute norms of polynomials. Also, resultants over number fields will arise naturally when solving systems of polynomial equations where the coefficients involve algebraic numbers. The main contributions we make in the paper are:

1. We generalize the modular algorithm of Collins [1] to compute resultants of multivariate polynomials over number fields with one or more field extensions (Sections 2 and 3).

2. We improve the efficiency of Collin’s algorithm in practice by using a probabilistic algorithm to compute the degree and also the low degree of the resultant in all variables. We give an error analysis for the case R = ℚ(α_{1})[y_{1}, y_{2}] (Section 4).

3. To avoid rational reconstruction in the recovery of the coefficients of the resultant, we pre-compute a multiple of the denominator of the resultant (Section 5).

4. We use primitive elements to convert multiple field extensions into one large extension after reducing modulo the prime for improved efficiency. The algorithm for the conversion uses gcds and resultants. We give a running time analysis for it (Section 6).

5. We have implemented the algorithms in a new recursive dense data structure in Maple [4] that is designed to speed up modular algorithms. We give some details of the data structure and implementation (Section 7).

We remark that for the case R = 𝔽_{q}[y_{1}, ..., y_{k}], where 𝔽_{q} is a finite field with q elements, the only difficulty in generalizing Collins’ algorithm is when q is too small, that is, there are not enough evaluation points in 𝔽_{q} to interpolate some y_{i}. This can be resolved by taking evaluation points from 𝔽_{q^{l}}, where q^{l} is chosen to be sufficiently large.

Before we look at the number field case we make some remarks about the efficiency of Collin’s algorithm in comparison with the subresultant algorithm. Specifically, we want to justify our use of a probabilistic algorithm. Consider the following problem for R = ℤ[y]. Suppose deg_{x}(f), deg_{x}(g), deg_{y}(f) and deg_{y}(g) are all n, and the length of the integer coefficients in f and g is O(n). Thus the size of f and g is O(n^{3}) and the size of res_{x}(f, g) is O(n^{4}) (neglecting logarithmic factors). One can show that the time complexity of Collin’s algorithm is O(n^{6}) and that of the subresultant algorithm is O(n^{10}). To give the reader a feel for the kind of improvement that one obtains on a specific example, for f, g ∈ ℤ[y][x] of total degree 10, with random 10 digit coefficients, the Maple implementation of Collin’s algorithm (`resultant/modular2`) is 63 times faster than the Maple implementation of the subresultant algorithm (`resultant/subres`). As impressive as these results appear, we make the following critical observations.

Suppose f = c_{n}x^{n} + c_{n−1}x^{n−1} + ... + c_{1}x + c_{0} and g = f + 1. If c_{n} = 1 then res_{x}(f, g) = 1. If the subresultant algorithm is used it will terminate in one step since rem_{x}(f, g) = 1. If Collin’s algorithm is used the bounds for the number of primes and evaluation points will depend on the size of the coefficients and n. We can see that Collin’s algorithm will take arbitrarily longer than the subresultant algorithm on this class of problems. Do we care about such degenerate cases? We say yes, for it has been our experience that the degree bounds and coefficient bounds on many real problems are often very bad. This increases the cost of the modular resultant algorithm but has no effect on the subresultant algorithm. Space prevents us from properly substantiating this claim, however, we give two examples.

Example 1.1 (Cyclotomic Polynomials^{1}): Let Φ_{n}(x) denote the n’th cyclotomic polynomial. Let f = Φ_{n}(x) and g = Φ_{m}(x). For 100 < m < n ≤ 200 we find that res_{x}(f, g) = 1! But for m = 197 and n = 199, Hadamards bound on the determinant of Sylvester’s matrix is a 453 digit integer! Even for the exceptional case of m = 100, n = 200 where the resultant res_{x}(f, g) = 2^{40}, Hadamards bound is over three times longer than the resultant.

Example 1.2 (The Logistic Map^{2}): Another way to construct resultants with degree and coefficient bounds which are poor estimates is to compose polynomials. Let f(x) = ax(1 − x) for parameter a. Let g(x) = f^{(n)}(x). The bifurcation points of f are the roots of the polynomial r(a) = res_{x}(g(x)−x, g′(x)+1). These resultants have much smaller coefficients than expected and are divisible by a large power of a.

1This example is given to us by Chris Smith, Department of Mathematics, University of Edinburgh. 2This example is given to us by Robert Corless, Department of Applied Mathematics, University of Western Ontario.

To improve the performance of the modular resultant algorithm we will require that the number of primes and the number of evaluation points used is proportional to the size of res_{x}(f, g) and not dependent on the inputs f and g. To achieve this we will use a probabilistic approach for both the coefficients and the degree of the resultant.

We make one further observation. In our experience, if R is a polynomial ring, res_{x}(f, g) is almost always denser than f and g. This contrasts with the gcd problem where the gcd(f, g) is usually sparser than the f and g. The reason for this is the following: when R is a polynomial ring, the resultant is an integer linear combination of products of polynomials and the product of two polynomials is usually denser than the inputs. Therefore, whereas sparse algorithms have practical advantages for the gcd problem, they are of much less value for the resultant problem. Consequently, we do not pursue a sparse approach here.

2. RESULTANTS OVER NUMBER FIELDS

Let deg_{x}(f) denote the degree of f in x and lc_{x}(f) denote the leading coefficient of f in x. Let p be a prime and let ϕ_{p} : R[x_{1}, …, x_{k}] → R_{p}[x_{1}, …, x_{k}] denote the modular homomorphism ϕ_{p}(f) = f mod p.

Let ℚ[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩ denote a number field over ℚ with t extensions. Each z_{i} is an extension variable and each m_{i}(z_{i}) ∈ ℚ[z_{1}, …, z_{i−1}]/⟨m_{1}, …, m_{i−1}⟩[z_{i}] is the minimal polynomial for z_{i}. Likewise let ℤ_{p}[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩ denote a finite ring over ℤ_{p} with t extensions.

Figure 1 is a homomorphism diagram for computing multivariate resultants over number fields. It provides an overview of our modular algorithm. We fix some of the details:

NRES computes multivariate resultants over number fields. N is the number of primes whose product must be large enough to reconstruct the rational coefficients of the resultant. In the algorithm, rather than computing a bound on the size of the coefficients, we simply run until the current value for the resultant does not change for five additional primes.

CRES computes multivariate resultants over finite rings by collapsing extensions. σ_{i}(z̃) is the substitution for z_{i} in term of z̃. It is discussed in Section 6.

PRES computes multivariate resultants over finite rings with zero or more extensions. The algorithm over ℤ_{p} is from [1] but extended to support extensions over ℤ_{p}. It uses dense evaluation and interpolation where D_{i} ≥ deg_{y_{i}}(r) is a degree bound on the degree of y_{i} of the resultant.

URES computes univariate resultants over rings with division using the Euclidean algorithm. It is discussed in Section 3.

REMARK In the diagram the inputs of PRES and URES are over ℤ_{p}[z̃]/⟨m̃⟩ (a finite ring with only one extension). This is because CRES collapses many extensions into one extension. However, our implementation of PRES and URES support multiple extensions, including finite fields.

ℚ[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩[y_{1}, …, y_{k}][x] ⟶^{NRES} ℚ[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩[y_{1}, …, y_{k}]
ϕ_{p}↓^{p∈{p_{1},...,p_{N} }}  ^{Chinese Remainder}↑^{Rational Reconstruction}
ℤ_{p}[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩[y_{1}, …, y_{k}][x] ⟶^{CRES} ℤ_{p}[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩[y_{1}, …, y_{k}]
|_{z_{i}=σ_{i}(z̃)}↓^{i=1,... ,t}  ^{|_{z̃=z_{1}+···+z_{t}}} ↑
ℤ_{p}[z̃]/⟨m̃⟩[y_{1}, …, y_{k}][x]  ⟶^{PRES}  ℤ_{p}[z̃]/⟨m̃⟩[y_{1}, …, y_{k}]
|_{y_{i}=β_{ij} ∈ℤp[z̃]/⟨m̃⟩}↓^{i=1,... ,k; j=0,... ,D_{i}}  ^{Polynomial}↑^{Interpolation}
ℤ_{p}[z̃]/⟨m̃⟩[x]  ⟶^{URES}  ℤ_{p}[z̃]/⟨m̃⟩

Figure 1: Homomorphism Diagram for Multivariate Resultants

3. THE EUCLIDEAN ALGORITHM OVER A RING

Algorithm URES implements the Euclidean algorithm over any field F. It computes the Euclidean polynomial remainder sequence of f, g ∈ F [x], or simply the remainder sequence, r_{0}, …, r_{n+1} ∈ F [x] satisfying

r_{0} = f,
r_{1} = g,
r_{i+1} = rem_{x}(r_{i−1}, r_{i}), and
r_{n+1} = rem_{x}(r_{n−1}, r_{n}) = 0.

Let d_{0}, …, d_{n} be the degrees of r_{0}, …, r_{n} in x. Then we can compute the resultant of f, g as follows [2]:

res_{x}(f, g) = lc_{x}(r_{n})^{d_{n−1}} ∏_{i=0}^{n−2}(−1)^{d_{i}d_{i+1}}lc_{x}(r_{i+1})^{d_{i}−d_{i+2}}

To avoid the coefficient swell, it is desirable to apply the Euclidean algorithm over a finite ring by first reducing the inputs f and g modulo some prime p. If F is a number field, since a minimal polynomial may become reducible modulo a prime, algorithm NRES may call URES over a finite ring with zero divisors which can cause URES to fail. For example, consider the finite ring ℤ_{5}[u]/⟨m⟩ where m = 3u^{3} − 1 is the minimal polynomial^{3}. Note that over ℤ_{5}, m = 3(u + 2)(u^{2} + 3u + 4) and thus ℤ_{5}[u]/⟨m⟩ is not a field. Let f = x^{3} + u^{2}x^{2} + 1 and g = x^{2} + x + u + 4 be polynomials in ℤ_{5}[u]/⟨m⟩[x]. We compute res_{x}(f, g) by URES as follows,

r_{0} = f = x^{3} + u^{2}x^{2} + 1
r_{1} = g = x^{2} + x + u + 4
r_{2} = rem_{x}(r_{0}, r_{1}) = 4(u + 2)(u + 4)x + u^{2} + u + 3
r_{3} = rem_{x}(r_{1}, r_{2}) = ⋯

3By the common definition, minimal polynomials are monic and irreducible. For efficiency, instead of monic polynomials with fractional coefficients, we always work with their associates over ℤ. For brevity, instead of characteristic polynomials, we call them minimal polynomials even when they are reducible, in case of finite rings.

The computation of r_{3} fails because lc_{x}(r_{2}) = 4(u + 2)(u + 4) is a zero divisor. If this occurs, a failure is generated; we simply trap the failure in NRES and try another prime. If the primes are large (typically, our primes are greater than 10^{9}), the proportion of zero divisors of the ring will be small, hence, failures will occur infrequently. The following theorems guarantee the correctness and termination of the modular resultant algorithm NRES for one field extension and one variable x. They can be generalized to a number field with k > 1 extensions. Proofs are based on [4].

Theorem 1 (Correctness). Let f = a_{s}x^{s}+⋯+a_{0}, g = b_{t}x^{t}+⋯+b_{0} ∈ ℚ(α)[x] = ℚ[z]/⟨m⟩[x] where m ∈ ℤ[z] is the minimal polynomial for α and where s, t > 0 and a_{s}, b_{t} ≠ 0. Let p be a prime such that

1. ϕ_{p}(lc_{z}(m)) ≠ 0,

2. f̄ = ϕ_{p}(f) and ḡ = ϕ_{p}(g) exist,

3. ϕ_{p}(a_{s}) ≠ 0 and ϕ_{p}(b_{t}) ≠ 0, and

4. URES(f̄, ḡ) terminates without failure.

Let r = res_{x}(f, g), r̄ = res_{x}(f̄, ḡ), and ω = URES(f̄, ḡ). Then r̄ = ϕ_{p}(r) and ω = ϕ_{p}(r).

Proof. Let M be the Sylvester matrix of f, g and M̄ be the Sylvester matrix of f̄, ḡ. Since ϕ_{p}(a_{s}) ≠ 0 and ϕ_{p}(b_{t}) ≠ 0, deg_{x}(f̄) = deg_{x}(f) = s and deg_{x}(ḡ) = deg_{x}(g) = t. Hence M and M̄ have the same dimension of s+t. By definition, each entry M̄_{ij} of M̄ is the image of the entry M_{ij} of M, that is, each M̄_{ij} = ϕ_{p}(M_{ij}). Hence M̄ = ϕ_{p}(M). Therefore, r̄ = det(M̄) = det(ϕ_{p}(M)) = ϕ_{p}(det(M)) = ϕ_{p}(r) because det is homomorphic. In addition, if URES(f̄, ḡ) terminates, ω = r̄ = ϕ_{p}(r) because all arithmetic operations in URES are homomorphic.

Lemma 2 (Existence of Polynomials). Let f ∈ ℚ[x] or f ∈ ℚ[z]/⟨m⟩[x] where m is the minimal polynomial. Let p be a prime. Let D : ℚ(α)[x] → ℤ denote the integer denominator, which is the smallest integer such that D(f)×f ∈ ℤ[α][x]. Then ϕ_{p}(f) exists if and only if D(f) ≢ 0 mod p.

Lemma 3 (Non-trivial Elements). Let f = a_{n−1}α^{n−1} + ⋯ + a_{1}α + a_{0} ∈ ℚ(α) where n is the degree of α. Let p be a prime such that ϕ_{p}(f) exists. Let each u_{i} denote the numerator of a_{i}. Let C : ℚ(α) → ℤ denote the greatest common divisor of u_{0}, …, u_{n−1}. Then ϕ_{p}(f) ≠ 0 if and only if C(f) ≢ 0 mod p.

Lemma 4 (Existence of Inverses). Let f ∈ ℚ(α) and p be a prime. Let N : ℚ(α) → ℚ denote the norm function. Then ϕ_{p}(f)^{−1} exists if and only if N(f) ≢ 0 mod p.

Theorem 5 (Termination). If a prime satisfies the conditions in Theorem 1 we say it is a good prime, otherwise, we say it is a bad prime. All but finitely many primes are good.

Proof. Intuitively, if one would apply the Euclidean algorithm directly over ℚ(α), the computation would be finite and hence the number of divisions is also finite. For each of the divisions there are only finitely many primes modulo which the inverse of an element of ℚ(α) does not exist. Therefore, all but finitely many primes are good.

By Lemma 2, f̄, ḡ exist if and only if D(f), D(g) ≢ 0 mod p. By Lemma 3, ϕ_{p}(a_{s}), ϕ_{p}(b_{t}) ≠ 0 if and only if C(a_{s}), C(b_{t}) ≢ 0 mod p. Since all lc_{z}(m), D(m), D(f), D(g), C(a_{s}) and C(b_{t}) are finite integers, there are only finitely many primes that divide any of them. Hence the first three conditions hold for all but finitely many primes.

Let r_{0}, …, r_{n+1} ∈ ℚ(α)[x] be the remainder sequence of f, g, as described at the beginning of this section. Let r̄_{0}, …, r̄_{k+1} ∈ ℤ_{p}(α)[x] be the remainder sequence of f̄, ḡ if the sequence exists. We prove by induction on k that, for all but finitely many primes, r̄_{i} exists and r̄_{i} = ϕ_{p}(r_{i}) for i = 0, …, n + 1.

For i = 0 and i = 1, r̄_{0} = f̄_{s} = ϕ_{p}(a_{s}) = ϕ_{p}(r_{0}) and r̄_{1} = ḡ_{t} = ϕ_{p}(b_{t}) = ϕ_{p}(r_{1}). Assume the statement is true for all i ≤ n. Given r̄_{i−1} and r̄_{i} we obtain r̄_{i+1} by computing the remainder of r̄_{i−1} by r̄_{i}. Let c_{i} be the leading coefficient of r_{i} in x. r̄_{i+1} exists and r̄_{i+1} = ϕ_{p}(r_{i+1}) if and only if the inverse of ϕ_{p}(c_{i}) exists. By Lemma 4, ϕ_{p}(c_{i})^{−1} exists if and only if N(c_{i}) ≢ 0 mod p. Since N(c_{i}) is a finite rational, there are only finitely many primes that divide it. Hence, for all but finitely many primes, r̄_{i} exists and r̄_{i} = ϕ_{p}(r_{i}) for i = 0, …, n + 1.

Theorem 1 tells us that if algorithm URES terminates without failure for a prime p, it outputs ϕ_{p}(r), the resultant modulo p. Theorem 5 tells us that URES terminates without failure for all but finitely many primes. Therefore, algorithm NRES will eventually have the resultant modulo enough primes to reconstruct the resultant over ℚ.

Below is the pseudo-code for NRES assuming that the input polynomials are over ℤ, the minimal polynomials are monic over ℤ and the coefficients are multivariate polynomials in Y = y_{1}, ..., y_{k}. In Section 5 we extend it to ℚ. Algorithm PDEG applies a probabilistic algorithm to obtain sharp degree bounds, which is discussed in Section 4.

Algorithm NRES: Number Field Resultant.

Input: f, g ∈ K[Y][x]\{0}, K = ℤ[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩ Output: r = res_{x}(f, g) ∈ K[Y]

1. {initialization}

(a) {probabilistic degree bounds} Π_{θ} ← PDEG(f, g)

(b) {big modulus prime} p ← 10^{8}

(c) {image agreement count} η ← 0

(d) {resultant so far} r ← null

2. {main computation} while η < 5 do

(a) {find next good prime} do p ← the next prime after p until ϕ_{p}(lc_{x}(f)), ϕ_{p}(lc_{x}(g)) ≠ 0 and ϕ_{p}(lc_{z_{1}}(m_{1})), …, ϕ_{p}(lc_{z_{t}}(m_{t})) ≠ 0

(b) {solve image problem} r̄ ← CRES(ϕ_{p}(f), ϕ_{p}(g), Π_{θ}) if r̄ = FAIL then goto step 2a.

(c) {lift modulus images} r_{0} ← CHINESE-REMAINDER([r, r̄])

(d) {stop if five agreements in a row} if r_{0} = r then η ← η + 1 else η ← 1 and r ← r_{0}

3. output r

4. DEGREE BOUNDS

Algorithm PRES requires a degree bound on each variable of the resultant in order to interpolate that variable. Because degree bounds derived from Sylvester’s matrix can be off, we will develop a probabilistic algorithm to obtain sharp degree bounds before we compute the resultant modulo a prime. We make a further observation. Very often in practice the resultant is divisible by large powers of some variables. Let r = res_{x}(f, g), d = deg_{y}(r) and ldeg_{y}(r) be the largest integer l such that y^{l}|r. If we know l in advance, it is easy to modify the interpolation process so that we need only d − l + 1 points to interpolate y.

Let us illustrate our algorithm with an example over ℤ. Let f(x, y, z) = y^{2}z^{2}x − y^{3}z^{2}x^{2} + 2 y^{3}zx^{3} − y^{2}zx^{2} − y^{3}x^{4} − x, g(x, y, z) = y^{2}z^{2} − 2 y^{3}z^{2}x + 6 y^{3}zx^{2} − 2 y^{2}zx − 4 y^{3}x^{3} + 1 and let r(y, z) = res_{x}(f, g). We have

r = (y^{2}z^{2} + 1) (y^{2}z^{2} + 2 yz − 5)^{2} (y^{2}z^{2} + 4 yz + 5) y^{9}.

We compute r̄ = res_{x}(ϕ_{p} f(x, y, α), ϕ_{p} g(x, y, α)) for a suitable prime p and α chosen at random from ℤ_{p} using degree bounds for deg_{y}(r(y, z)) based on Sylvester’s matrix.^{4} Then deg_{y}(r̄) ≤ deg_{y}(r) and ldeg_{y}(r̄) ≥ ldeg_{y}(r). If p is large then the probability of equality here is high. We repeat this to obtain degree estimates for deg_{z}(r) and ldeg_{z}(r).

4We could also use a probabilistic algorithm here, similar to that for NRES, such that we interpolate resultant images at different evaluation points until there are several agreements in a row.

In this example we have deg_{y}(r) = 17, deg_{z}(r) = 8, ldeg_{y}(r) = 11, and ldeg_{z}(r) = 0. Since the degree bounds we obtain from Sylvester’s matrix are deg_{y}(r) ≤ 20 and deg_{z}(r) ≤ 11, the number of univariate resultants computed over ℤ_{p} is reduced from (20 + 1) × (11 + 1) = 252 to (17 − 11 + 1) × (8 − 0 + 1) = 63 for a savings of a factor of 4. If we do not take into account ldeg_{y}r we would need (17 + 1) × (8 + 1) = 162 univariate resultants and our savings is a factor of only 1.56.

Let us determine the probability of obtaining a wrong degree bound for one field extension, i.e. for K = ℚ(α). Let S = ℤ_{p}[u]/⟨m⟩ be a finite ring where m ∈ ℤ_{p}[u] is the minimal polynomial of degree q. Let f, g ∈ S[y, z][x]. We want to apply the probabilistic algorithm to obtain the degree bound on y of r = res_{x}(f, g) ∈ S[y, z]. Let deg_{y}(r) = d. We can write the resultant as r = c_{d}(z)y^{d}+⋯+c_{1}(z)y+c_{0}(z) where each c_{i}(z) ∈ S[z]. Let deg_{z}(c_{d}) = n. Let f̄, ḡ ∈ S[y][x] be f, g evaluated at z = α for a random evaluation point α ∈ S. Let r̄ = res_{x}(f̄, ḡ) ∈ S[y] with deg_{y}(r̄) = d̄. Let E denote the probability that the degree bound is wrong, that is, E = Prob(d̄ < d). This happens when α is a root of c_{d}(z). The following lemma tells us what E is and suggests we should do.

Lemma 6 (Number of Roots in Finite Rings). Let p be a prime, m ∈ ℤ_{p}[u] of degree q > 0, and S = ℤ_{p}[u]/⟨m⟩. Let f ∈ S[z] have degree n. Let N be the number of roots of f in S.

1. If m is irreducible, then N ≤ n.

2. If m is square-free, then N ≤ n^{q}.

3. Otherwise N ≤ np^{q−1}.

REMARK The bound in case 3 cannot be improved significantly in general for the polynomial z^{q} has p^{q−1} roots in ℤ_{p}[u]/⟨u^{q}⟩.

Proof. For Case 1, S is a finite field. The statement is equivalent to the Fundamental Theorem of Algebra.

For Case 2, S is a simple finite ring. The worst case for E occurs when m has q distinct linear factors, that is, m = (u − β_{1}) ⋯ (u − β_{q}). The statement then follows from the Chinese remainder theorem.

For Case 3, S is an arbitrary finite ring. The worst case for E occurs when m has only one linear factor of multiplicity q, that is, m = (u − β)^{q}. We present a counting argument on the number of roots using p-adic lifting as follows. Let S^{i} = ℤ_{p}[u]/⟨(u − β)^{i}⟩ for i = 1, …, q. Consider the base case when q = 1. Since u − β is irreducible as in Case 1, f has at most n roots in S^{1}. Let α_{1}, …, α_{n} be the roots of f in S^{1}. Using p-adic lifting, we can lift these roots from S^{1} to S^{2} = ℤ_{p}[u]/⟨(u − β)^{2}⟩ so that f has at most np^{1} roots in S^{2} constructed by p-adic lifting. We can continue to lift to S^{q} = ℤ_{p}[u]/⟨(u − β)^{q}⟩ so that f has at most np^{q−1} roots in S^{q} = S constructed by p-adic lifting.

We cannot always choose a prime p such that m is irreducible modulo p (Case 1). But we can with high probability choose p such that m is square-free modulo p (Case 2). Case 2 does not hold when r = res_{u}(m,′ m) ≡ 0 mod p. Since r is an integer of a finite length, Case 2 holds for all but finitely many primes.

Theorem 7 (Error Probability in Simple Rings). Let S be as in Lemma 6. If m is square-free then E ≤ (n/p)^{q}.

Proof. As above, let c_{d} ∈ S[z] of degree n be the leading coefficient of r in y. Note that there are p^{q} possible evaluation points in S. Since S is a simple finite ring, by Lemma 6, c_{d} has at most n^{q} roots in S. Therefore, E ≤ n^{q}/p^{q} = (n/p)^{q}.

This is a good result. For example, for n = 100, p > 10^{8}, and q = 4, the probability of an error is less than 1/10^{24}. Repeating this with a second prime q > 10^{8} we would reduce the error to less than 1/10^{48}.

REMARK When the inputs are over ℚ, it is possible that c_{d}(z) becomes 0 modulo the prime we choose and thus the degree bound is wrong for all random points α ∈ S. This happens when the prime divides all the coefficients of c_{d}(z). Hence, instead of using the same prime but different β in algorithm PDEG, we use new primes to reduce the probability that this happens.

REMARK When the inputs f, g ∈ S[y_{1}, ..., y_{k}][x] where k > 2, Lemma 6 and Theorem 7 can be extended by replacing S[z] in Lemma 6 and Theorem 7 by S[z_{1}, ..., z_{k−1}] and generalizing the result of Schwartz [8] on the number of roots of f ∈ ℤ_{p}[z_{1}, ..., z_{k−1}] to f ∈ S[z_{1}, ...., z_{k−1}].

5. DENOMINATOR MULTIPLES

In Section 2, NRES computes resultants when the coefficients of the input and the minimal polynomials are over ℤ. In order to extend NRES to work over ℚ, we first make the inputs primitive over ℤ, that is, we clear the fractions from the input by multiplying through by a suitable integer and divide out by any common factor of the numerators of the rational coefficients. However, the denominator problem with the minimal polynomials is not as easy to resolve. Consider the number field ℚ[w, z]/⟨m_{1}, m_{2}⟩ where m_{1} = 3w^{3}−1 ∈ ℚ[w] and m_{2} = 2z^{3}−3zw^{2}+w+3 ∈ ℚ[w]/⟨m_{1}⟩[z]. Let f = x^{3} + x + 2wz and g = x^{2} + (1 + w)x + w + z be polynomials in ℚ[w, z]/⟨m_{1}, m_{2}⟩[x]. Note that the coefficients of f, g are over ℤ but those of res_{x}(f, g) = \frac{7}{6} + 2 z + (2 − \frac{20}{3} z + \frac{7}{2} z^{2}) w +(−2 + 9 z + 10 z^{2}) w^{2} are not. The denominator 6 of the resultant is introduced from divisions by the minimal polynomials which have non-unit leading coefficients. To recover a resultant with fractions we can apply either rational reconstruction or, if we know a multiple of the denominator, integer reconstruction. Both reconstructions have a quadratic time complexity^{5}, but the former has

5Rational reconstruction has a cubic time complexity if it is applied after each image; [4] point out that, if it is applied after each F_{i} image (the Fibonacci numbers F_{i}: 1, 2, 3, 5, 8, 13, …), the reconstruction has a quadratic complexity. a big constant. Hence, we will use integer reconstruction in the following.

To compute a multiple of the denominator 6 in the example, we start with the minimal polynomial of the last extension, which is m_{2}(w, z) in the example. First, we compute the degree bound on z of the resultant, d_{z} = 3 using the probabilistic degree bound^{6} in Section 4. Then we know there is at most κ_{z} = d_{z} − deg_{z}(m_{2}) + 1 = 3 − 3 + 1 = 1 division by m_{2}. Hence the contribution of the extension m_{2}(w, z) to the denominator multiple is lc_{z}(m_{2})^{−κ_{z}} = 1/2.

But, since m_{2} is in ℚ[w]/⟨m_{1}⟩[z], the division by m_{2} also contribute to the degree in w. Let σ_{w} denote this accumulated degree from m_{2}(w, z) for w. Working out the details^{7}, σ_{w} = deg_{w}(λ_{w})⌈κ_{z}/(deg_{z}(m_{2}) − deg_{z}(λ_{w}))⌉ = 2⌈2/(3 − 1)⌉ = 2 where λ_{w} = −3zw^{2} is the leading term of m_{2} in w. In general, we update the accumulated degrees of all other extension variables recursively.

Now we can proceed to compute the contribution of the extension m_{1}(w) with the accumulated degree σ_{w} = 2. Again, we compute the degree bound on w of the resultant, d_{w} = 4. Then, we compute the bound on the number of divisions by m_{1}, κ_{w} = (σ_{w} + d_{w}) − deg_{w}(m_{1}) + 1 = (2 + 4) − 3 + 1 = 4. Hence the denominator multiple from all extensions is lc_{w}(m_{1})^{−κ_{w}}lc_{z}(m_{2})^{−κ_{z}} = 1/(3^{4} × 2^{1}) = 1/162.

Here are the recursive formulas for the general case. Let f, g ∈ ℚ[z_{1}, …, z_{t}]/⟨m_{1}, …, m_{t}⟩[Y][x]. Let θ_{i} denote the degree and c_{i} denote the leading coefficient of m_{i} in z_{i}. Let λ_{j}^{k} denote the term of degree k of m_{j} and λ_{ij}^{k} denote the degree of λ_{j}^{k} in z_{i}. Let d_{i} denote the degree bound on z_{i} of res_{x}(f, g) . Let σ_{i} denote the accumulated degree and κ_{i} denote the bound on the number of divisions by m_{i}. Let η_{ij}^{k} denote the bound on the number of divisions by m_{i} due to λ_{j}^{k}. Let D denote the denominator multiple of res_{x}(f, g). Then,

η_{ij}^{k} = λ_{ij}^{k}⌈κ_{j}/(θ_{j} − k)⌉
σ_{i} = ∑_{j=i+1}^{t} max(η_{ij}^{0}, …, η_{ij}^{θ_{j} −1})
κ_{i} = max(0, σ_{i} + d_{i} − θ_{i} + 1)
D = c_{1}^{−κ_{1}} ⋯ c_{t}^{−κ_{t}}

6. COLLAPSING EXTENSIONS

Suppose the number field is of degree n = 2^{m}. The modular resultant algorithm will be considerably faster if the number field is presented by a single field extension of degree 2^{m} than by m quadratic field extensions. This is because of the overhead of storage management of small vectors and program control when working with extensions of low degree

6by ignoring the minimal polynomials and treating extension variable as polynomial variables 7consider the Newton’s polygon of the monomials of w, z and the monomial of the steepest slope; for example, given m_{1} = 3w^{7}+2 and m_{2} = 5z^{5}+w^{4}z^{4}+w^{6}z^{3}+1, the term w^{4}z^{4} “dominates” the term w^{6}z^{3} in calculating the accumulated degree for w in comparison with the arithmetic operations in ℤ_{p} that are being done. This overhead can be measured. We have timed Maple’s overhead (see Appendix) for polynomial multiplication and gcd operations in ℤ_{p}[z]/⟨m⟩[x]. We determined that the overhead for deg_{z}(m) = 4 is a factor of 15 to 16 and the overhead for deg_{z}(m) = 8 is a factor of 4 to 5. Therefore, we can get an non-asymptotic, but still very significant speed up for resultant computations by collapsing multiple extensions into one large extension. For a typical problem with a quadratic extension, this can speed up the entire resultant computation by a factor of 10. Although this causes severe expression swell over ℚ, the conversion is rather efficient over ℤ_{p}. We give an O(n^{3}) conversion algorithm. Another reason to collapse the extensions into one large extension would be to take advantage of fast multiplication and division in ℤ_{p}[z].

Let us illustrate the conversion algorithm with an example. Consider a finite ring ℤ_{7}[w, z]/⟨m_{1}, m_{2}⟩ where m_{1} = 3w^{3} − 1 ∈ ℤ_{7}[w] and m_{2} = 2z^{3} − 3zw^{2} + w + 3 ∈ ℤ_{7}[w]/⟨m_{1}⟩[z]. Let f = x^{3} + x + 2wz and g = x^{2} + (1 + w)x + w + z be polynomials in ℤ_{7}[w, z]/⟨m_{1}, m_{2}⟩[x]. We first compute a representation of the primitive element u = w + z of the two extensions^{8}, that is, we compute the minimal polynomial for u and the substitutions of w, z in u. One way to do this is to use linear algebra. Another, potentially faster way, is to use resultants and gcds. The idea, attributed to Collins by Loos in [6], is as follows. First, by computing a sequence of resultants, we eliminate the extension variables w, z from the minimal polynomials for w, z and the defining equation u − w−z. The last resultant in the sequence will be the minimal polynomial for u. Then we can recover the substitutions of w, z by computing in the new field ℤ_{7}[u]/⟨m⟩ the gcds of the minimal polynomials and the intermediate resultants.

Let us apply the algorithm to the example. First we compute a sequence of resultants as follows:

r_{2} = u − w − z mod ⟨m_{1}, m_{2}⟩ ∈ ℤ_{7}[u, w][z]
r_{1} = res_{z}(r_{2}, m_{2}) mod ⟨m_{1}⟩ ∈ ℤ_{7}[u][w]
= 4uw^{2} + (6 + 6u^{2})w + 5u^{3} + 6
r_{0} = res_{w}(r_{1}, m_{1}) ∈ ℤ_{7}[u]
= 6u^{9} + 3u^{6} + 5u^{4} + u^{3} + 6u^{2} + 3u + 1

The minimal polynomial for u is m = r_{0}. Let gcd_{x}(f, g) denote the monic greatest common divisor of f, g in x. We compute gcd_{w}(g, r_{1}) mod ⟨m⟩ = w−(3u^{8}+5u^{7}+2u^{6}+2u^{5}+ 3u^{4} + 3u^{2} + 3u + 6) and solve for w in ℤ_{7}[u]/⟨m⟩. Hence the substitution of w in u is 3u^{8} + 5u^{7} + 2u^{6} + 2u^{5} + 3u^{4} + 3u^{2} + 3u + 6 and that of z is u − w = 4u^{8} + 2u^{7} + 5u^{6} + 5u^{5} + 4u^{4} + 4u^{2} + 5u + 1 in ℤ_{7}[u]/⟨m⟩.

We make the substitution into f, g to obtain f̃ = x^{3} + x + 2u^{8}+2u^{7}+3u^{5}+3u^{4}+u^{3}+5u^{2}+6u+3, g̃ = x^{2}+(4u^{8}+2u^{7}+ 5u^{6} +5u^{5} +4u^{4} +4u^{2} +5u+2)x+u ∈ ℤ_{7}[u]/⟨m⟩[x]. Then we proceed in this ring of one extension of degree 9 to obtain r̃ = res_{x}(f̃, g̃) = 2u^{8} + 4u^{7} + 6u^{4} + 3u^{2} + 6u + 4 ∈ ℤ_{7}[u]/⟨m⟩. At last we back-substitute u = w + z into r̃ to obtain r = res_{x}(f, g) = 3w^{2}z^{2} + 2wz^{2} + 5wz + 5z^{2} + 5w + 2z + 5 ∈ ℤ_{7}[w, z]/⟨m_{1}, m_{2}⟩.

8At the end of this section we say what to do when u = w+z is not a primitive element.

Theorem 8 below is taken from Loos [6]. It can be generalized to a number field with an arbitrary number of extensions by iterating the procedure, that is, eliminating the last field extension until there is only one left.

Theorem 8 (Correctness). Let F [w, z]/⟨m_{1}, m_{2}⟩ be a number field where m_{1} ∈ F [w] of degree s and m_{2} ∈ F [w]/⟨m_{1}⟩[z] of degree t are the minimal polynomials. Let

r_{2} = u − w − z mod ⟨m_{1}, m_{2}⟩ ∈ F [u, w][z]
r_{1} = res_{z}(r_{2}, m_{2}) mod ⟨m_{1}⟩ ∈ F [u][w]
r_{0} = res_{w}(r_{1}, m_{1}) ∈ F [u]

Let m = r_{0}. If u = w+z is a primitive element then m is the minimal polynomial for u and F [u]/⟨m⟩ ≅ F [w, z]/⟨m_{1}, m_{2}⟩. Let σ : F [w, z]/⟨m_{1}, m_{2}⟩ → F [u]/⟨m⟩ be the mapping of the isomorphism. Then σ is given by

gcd_{w}(m_{1}, r_{1}) mod ⟨m⟩ = w − σ(w) ∈ F [u]/⟨m⟩[w].

Proof. Let α_{1}, …, α_{s} be the distinct roots of m_{1} and β_{1}, …, β_{t} be the distinct roots of m_{2}. By an alternative definition of resultants [2],

r_{1} = ∏_{1≤j≤t}(u − w − β_{j}) mod ⟨m_{1}⟩
= (u − w − β_{1}) ⋯ (u − w − β_{t}) ∈ F [u][w]
r_{0} = ∏_{1≤i≤s} ∏_{1≤j≤t}(u − α_{i} − β_{j})
= (u − α_{1} − β_{1}) ⋯ (u − α_{s} − β_{t}) ∈ F [u]

Let each γ_{ij} = α_{i} + β_{j} be a root of r_{0} = m. Since α_{i}, β_{j} are distinct, there are s×t distinct roots of m. Since both fields have the same dimension, F [u]/⟨m⟩ ≅ F [w, z]/⟨m_{1}, m_{2}⟩.

By construction, σ(w) = α_{i} = γ_{ij} − β_{j} for some i, j. Hence w − σ(w) = w − α_{i} is a factor of m_{1} and w − σ(w) = w − (γ_{ij} − β_{j}) = −(γ_{ij} − w − β_{j}) is also a factor of r_{1}. Suppose w − α̃ is another common factor of r_{1}, m_{1}. Since any factor of r_{1} is of the form w − (γ_{ij} − β̃) for some β̃, α̃ = γ_{ij} − β̃. This implies α̃ = α_{i} = σ(w) because α_{i}, β_{j} are distinct. Therefore, gcd_{w}(m_{1}, r_{1}) = w − σ(α).

The following theorem is stated for the equal degree case to obtain a meaningful measure for the time complexity of the conversion.

Theorem 9 (Time Complexity). Let deg_{w}(m_{1}) = \sqrt{n} and deg_{z}(m_{2}) = \sqrt{n}. Then deg_{u}(m) = [F [w, z]/⟨m_{1}, m_{2}⟩ : F] = n and the time complexity of the conversion is O(n^{3}) arithmetic operations in F.

REMARK This is the same complexity that one obtains if one uses linear algebra. However, it will be easier to obtain an asymptotic speedup using asymptotically fast polynomial operations than using asymptotically fast matrix operations. Moreover, we think an O(n^{2}) conversion algorithm may be possible using this approach. In particular, the minimal polynomial is computed in O(n^{2}).

Proof. There are three main computations in the conversion. The first two are computing r_{1} and gcd(m_{1}, r_{1}) mod ⟨m⟩ in Theorem 8. The third is computing the other substitutions σ(w^{i}z^{j}), for i, j = 0, …, \sqrt{n} − 1, given σ(w), σ(z). Note, σ(z) is easily computed from u = w + z.

Let us first determine the time complexity for computing r_{0} = res_{w}(r_{1}, m_{1}) ∈ F [u]. Note that deg_{u}(r_{0}) = n and deg_{w}(r_{1}) = deg_{w}(m_{1}) = \sqrt{n}. To interpolate the monic r_{0} we need n points for u, taking O(n^{2}) operations in F. There are n image problems in F [w]. The Euclidean algorithm in URES takes O(\sqrt{n} × \sqrt{n}) operations in F for each. Hence r_{0} can be computed in O(n^{2} + n × n) = O(n^{2}) operations in F.

To compute gcd_{w}(m_{1}, r_{1}) mod ⟨m⟩ ∈ F [u]/⟨m⟩[w], the Euclidean algorithm requires O(\sqrt{n}×\sqrt{n}) operations in F [u]/⟨m⟩ since deg_{w}(r_{1}) = deg_{w}(m_{1}) = \sqrt{n}. Using classical algorithms, arithmetic operations in F [u]/⟨m⟩ are O(n^{2}) since deg_{u}(m) = n. Hence computing σ(w) requires O(n × n^{2}) = O(n^{3}) operations in F.

Finally, given σ(w), σ(z), we compute other substitutions σ(w^{i}z^{j}) = σ(w)^{i}σ(z)^{j} mod ⟨m⟩ ∈ F [u]/⟨m⟩. This requires O(n) operations in F [u]/⟨m⟩, or O(n × n^{2}) = O(n^{3}) operations in F. Hence the overall time complexity of the conversion is O(n^{3}) operations in F.

Over Finite Rings: When the inputs are over finite rings, however, the conversion may fail because of division by a zero divisor. Let us illustrate the problem with an example. Consider a finite ring ℤ_{11}[w, z]/⟨m_{1}, m_{2}⟩ where m_{1} = w^{3} + 2 ∈ ℤ_{11}[w] and m_{2} = z^{2} + 1 ∈ ℤ_{11}[w]/⟨m_{1}⟩[z]. Note that m_{1} can be factored as (w^{2} + 4 w + 5) (w + 7) and thus ℤ_{11}[w, z]/⟨m_{1}, m_{2}⟩ is not a field. Let us try to convert it to ℤ_{11}[u]/⟨m⟩ anyway. To obtain the minimal polynomial for u, we compute the sequence of resultants as follows,

r_{2} = u − w − z mod ⟨m_{1}, m_{2}⟩ ∈ ℤ_{11}[u, w][z]
r_{1} = res_{z}(r_{2}, m_{2}) mod ⟨m_{1}⟩ ∈ ℤ_{11}[u][w]
= w^{2} + 9 wu + u^{2} + 1
r_{0} = res_{w}(r_{1}, m_{1}) ∈ ℤ_{11}[u]
= u^{6} + 3 u^{4} + 4 u^{3} + 3 u^{2} + 10 u + 5

Note that m = r_{0} = (u^{2} + 4 u + 8) (u^{2} + 3 u + 6) (u + 2)^{2} over ℤ_{11} and thus ℤ_{11}[u]/⟨m⟩ is not a field. To obtain the substitution of w in u, we compute gcd_{w}(m_{1}, r_{1}) by the Euclidean algorithm as follows,

s_{0} = m_{1} = w^{3} + 2
s_{1} = r_{1} = w^{2} + 9 wu + u^{2} + 1
s_{2} = rem_{x}(s_{0}, s_{1})
= 3(u + 2)(u + 9)w + 9 u^{3} + 9 u + 2
s_{3} = rem_{x}(s_{0}, s_{1}) = …

It fails in computing s_{3} because the leading coefficient of s_{2}, 3(u + 2)(u + 9), is a zero divisor of the ring.

Similar to computing resultants over finite rings in Section 3, the conversion will trap the failure and try another prime. Alternatively, we could use other primitive elements such as u = w + sz where s is an integer.

Primitive Elements: We have assumed that the sum of extension variables u = w + z is a primitive element. This is not true in general. For example, consider the number field ℚ(\sqrt{2}, \sqrt{3} − \sqrt{2}). Letting w = \sqrt{2} and z = \sqrt{3} − \sqrt{2}, if we define u = w + z = \sqrt{2} + (\sqrt{3} − \sqrt{2}) = \sqrt{3}, then we cannot represent \sqrt{2} in the new field ℚu. In the implementation, we check if substitutions of both extension variables exist and if not, try some other linear combination u = w + sz of the extension variables and repeat the conversion modulo another prime p. Since the primes is large, choosing s ∈ ℤ_{p} randomly from [100, …, p − 100] will give us a primitive element with high probability. Moreover, the cost is cheap since all arithmetic is being done modulo p.

REMARK If all the minimal polynomials of the input field are square-free over ℚ, we choose a prime and a linear combination of the extension variables such that the new minimal h(u) polynomial is square-free modulo the prime. This ensures the existence of the substitutions and ensures a low probability of obtaining wrong degree bounds in Section 4.

7. IMPLEMENTATION

The algorithms in this paper are implemented using the recden package in Maple [4]. The source codes and test examples can be downloaded from [10]. recden is a recursive dense representation of polynomials. Formally recden polynomials are encoded in Maple as follows:

  <poly>   ::= RPOLY( R::<ring>, C::<coeffs> )
  <ring>   ::= [ m::nonnegint, X::list(symbol),
                 E::list(<coeffs>) ]
  <coeffs> ::= 0 or list(<coeff>)
  <coeff>  ::= rational or <coeffs>

Thus a recden polynomial has two operands: the polynomial ring R and the coefficients C. The ring R consists of the modulus m (or the characteristic), the list of variables X, and the list of ring extensions over the base ring ℤ_{m} or ℚ. By convention, m = 0 means the ring is over ℚ. The data structure for C is recursive dense with the integer 0 indicating a zero coefficient. This implies a polynomial is a zero polynomial if and only if C = 0.

The other key requirement is that the bottom level of the data structure, which represents a polynomial in one variable, uses a dense array of coefficients. Thus the representation for ℤ_{p}[x], which is what most computations in the modular algorithms are reduced to, is an array of machine integers. Polynomial operations in ℤ_{p}[x] are carefully coded in Maple’s kernel in the C language for efficiency. See Monagan [7].

In the following we give examples of recden for polynomials over ℤ, a finite field, and a number field.

Example 7.1 (Integer Polynomial): Let f_{1} = x^{2} + y^{2} + 2 ∈ ℤ[y][x]. Its recden representation is:

      RPOLY( [0,[x,y],[]], [[2,0,1],0,[1]] )

Here m = 0 because f_{1} is over ℤ, X = [x,y] because x, y are the polynomial variables, and E = [] because there are no extensions. Also, C_{1} = [2,0,1] is the coefficient in ℤ[y] of x^{0}, that is, 2y^{0} + 0y^{1} + 1y^{2}. Similarly C_{2} = 0 (the integer 0) is the coefficient in ℤ[y] of x^{1} and C_{3} = [1] is that of x^{2}.

Example 7.2 (Finite Ring): Let F = ℤ_{3}[y]/⟨y^{4} + y + 1⟩ be a finite ring of cardinality 3^{4}. Let f_{2} = x^{2} +y^{2} +2 ∈ F [x]. Its recden representation is:

                                                          RPOLY([3,[x,y],[[1,1,0,0,1]]], [[2,0,1],0,[1]])

Here m = 3 because f_{2} is over ℤ_{3} and E_{1} = [1,1,0,0,1] is the encoding of the minimal polynomial 1y^{0} + 1y^{1} + 0y^{2} + 0y^{3} + 1y^{4}. Note that X and C are the same as those of f_{1} but they are interpreted differently. For f_{1} both x, y are polynomial variables, but for f_{2} only x is a polynomial variable while y is an extension variable. This is determined by the number of extensions in E. Likewise, for f_{1} we interpret C_{1} = [2,0,1] as a polynomial, but for f_{2} we interpret it as an element of a finite ring.

Besides the efficient representation, recden also provides a user-friendly interface. The input command for and the user display of f_{2} in Maple are:

                                                        > f2 := rpoly(x^2+y^2+2, [x,y], 3, y^4+y+1);

f_{2} := x^{2} + y^{2} + 2 mod ⟨y^{4} + y + 1, 3⟩

Example 7.3 (Number Field): Let F be the number field

ℚ(\sqrt[3]{5}, \sqrt{−2 − (\sqrt[3]{5})^{2}}). Note that F has two algebraic extensions. We can represent F by ℚ[w, z]/⟨w^{3} − 5, z^{2} + w^{2} + 2⟩. Let f_{3} = z + w + 1 ∈ F. Its input command, its user display and its recden representation are:

                                                        > f3 := rpoly(z+w+1, [z,w], [w^3-5,z^2+w^2+2]);

f_{3} := z + w + 1 mod ⟨z^{2} + w^{2} + 2, w^{3} − 5⟩

                                                        > lprint(f3);
                                                           RPOLY([0,[z,w],[[[2,0,1],0,[1]],[-5,0,0,1]]],
                                                                           [[1,1],[1]])

Here E_{2} = [-5,0,0,1] is the first extension with the minimal polynomial w^{3} − 5 and E_{1} = [[2,0,1],0,[1]] is the last extension with z^{2} + w^{2} + 2. Since the same encoding is used for any coefficients, the E_{1} of f_{3} (z^{2} + w^{2} + 2 ∈ ℚ[w]/⟨w^{3} − 5⟩[z]) is identical to the C of f_{1} (x^{2} + y^{2} + 2 ∈ ℤ[y][x]) and of f_{2} (x^{2} +y^{2} +2 ∈ ℤ_{3}[y]/⟨y^{4} + y + 1⟩[x]). That is, the encoding is independent of variables and base rings.

Code Snippet: To give the reader a taste of the implementation, we show the code for Algorithm URES in Figure 2. It computes the resultant of two polynomials using the Euclidean algorithm. The following example shows URES executing over the number field ℚ[y]/⟨y^{4} + y + 1⟩ then, after applying ϕ_{3} to the polynomials, over the finite ring ℤ_{3}[y]/⟨y^{4} + y + 1⟩.

URES := proc(ff,gg) local f,g,h,R,n,m,r,k;
## ff,gg are of positive degree in R[x]
(f,g) := (ff,gg);
(n,m) := (degrpoly(f),degrpoly(g));
R := projring(getring(f));
r := rpoly(1,R);
while m>0 do
h := remrpoly(f,g);
if iszerorpoly(h) then return h; fi;
k := degrpoly(h);
## r := r (-1)^(nm) lc(g)^(n-k)
r := mulrpoly(r,powrpoly(lcrpoly(g),n-k));
if irem(n*m,2)=1 then r := negrpoly(r) fi;
(f,g) := (g,h);
(n,m) := (m,k);
end do;
mulrpoly(r,powrpoly(lcrpoly(g),n));
end proc:

Figure 2: Maple Code for Univariate Resultants

> f := rpoly(x^2+y^2+2, [x,y], y^4+y+1);
> g := rpoly(x^3+x^2*y+x*y+2, [x,y], y^4+y+1);

f := x^{2} + y^{2} + 2 mod ⟨y^{4} + y + 1⟩ g := x^{3} + x^{2}y + xy + 2 mod ⟨y^{4} + y + 1⟩

> URES(f,g);

1 − 25y + 18y^{2} − 14y^{3} mod ⟨y^{4} + y + 1⟩

> URES(phirpoly(f,3),phirpoly(g,3));

1 + 2y + y^{3} mod ⟨y^{4} + y + 1, 3⟩

8. REFERENCES

[1] George Collins. The Calculation of Multivariate Polynomial Resultants. JACM 18(4), 1971.

[2] Keith Geddes, Stephen Czapor, and George Labahn. Algorithms for Computer Algebra, Kluwer Academic, 1996.

[3] A. Goldstein, G. Graham, A Hadamard type bound on the coefficients of a determinant of polynomials. SIAM Review, 16, pp. 394–395, 1974.

[4] Mark van Hoeij, Michael Monagan. A Modular GCD Algorithm over Number Fields Presented with Multiple Field Extensions. Proceedings of ISSAC ’02, pp. 109–116, 2002.

[5] Rüdiger Loos. Generalized Polynomial Remainder Sequences. Computer Algebra, Symbolic and Algebraic Computation, second edition, Springer Verlag, pp. 115–138, 1983.

[6] Rüdiger Loos. Computing in Algebraic Extensions. Computer Algebra, Symbolic and Algebraic Computation, second edition, Springer Verlag, pp. 173–188, 1983.

[7] Michael Monagan. In-place Arithmetic for Polynomials over ℤ_{n}. DISCO ’92, pp. 22–34, 1992 [8] Jacob Schwartz. Fast Probabilistic Algorithms for

Verification of Polynomial Identities. JACM27(4), 1980.

[9] Barry Trager, Algebraic Factoring and Rational Function Integration. Proceedings of ISSAC ’76, ACM Press, pp. 219–226, 1976.

[10] Stephen Chun-To Tse. Algorithms and Bounds for Resultants. Master Thesis, Computing Science Department, Simon Fraser University, 2002.

   http://www.sfu.ca/~stephent/resultant.

Appendix

The following timings show that, in Maple, when computing in R[x] where R = ℤ_{p}[z]/⟨m(z)⟩, there is a loss of efficiency of a factor of 15 to 16 more than expected in theory when the degree of m is 4 (compared with 256).

> p := 46327:
> k1 := 4;
                          k1 := 4
> m := z^4-10*z^2+1:
> c := () -> randpoly(z,degree=k1-1,dense):
> n1 := 100:
> f := randpoly(x,dense,degree=n1,coeffs=c):
> g := randpoly(x,dense,degree=n1,coeffs=c):
> m := modp1(ConvertIn(m,z),p):
> F := evalgf1(ConvertIn(f,z,x),m,p):
> G := evalgf1(ConvertIn(g,z,x),m,p):
> st := time():
> evalgf1(Multiply(F,G),m,p):
> mt1 := time()-st: st := time():
> evalgf1(Gcd(F,G),m,p):
> gt1 := time()-st:
> k2 := 256;
                         k2 := 256
> c := () -> randpoly(z,degree=k2-1,dense):
> m := z^128+c():
> n2 := 10:
> f := randpoly(x,dense,degree=n2,coeffs=c):
> g := randpoly(x,dense,degree=n2,coeffs=c):
> p := modp1(Prime(1)):
> m := modp1(ConvertIn(m,z),p):
> F := evalgf1(ConvertIn(f,z,x),m,p):
> G := evalgf1(ConvertIn(g,z,x),m,p):
> st := time():
> evalgf1(Multiply(F,G),m,p):
> mt2 := time()-st: st := time():
> evalgf1(Gcd(F,G),m,p):
> gt2 := time()-st:
> Loss1 := mt1/n1^2/k1^2  * n2^2*k2^2/mt2;
                    Loss1 := 14.86883354
> Loss2 := gt1/n1^2/k1^2  * n2^2*k2^2/gt2;
                    Loss2 := 16.09043949