Computing Parametric Geometric Resolutions

Éric Schost Laboratoire GAGE, École polytechnique F-91128 Palaiseau, France schost@gage.polytechnique.fr

Abstract

Consider a polynomial system f of n equations in n unknowns that depends on some parameters. Suppose that for a generic instantiation of the parameters this system has a finite number of solutions, and that these solutions are not singular. We propose a probabilistic algorithm to compute a geometric resolution of the generic fiber with coefficients that are rational functions of the parameters, give its complexity and show various examples of applications.

1 Introduction

This paper is concerned with the resolution of some polynomial systems with parameters. Given an effective field k, of characteristic either zero or big enough, we will consider polynomials f_{1}, …, f_{n} in k[P_{1}, …, P_{m}, X_{1}, …, X_{n}]. The variables P are thought as parameters, and the variables X as unknowns.

We suppose that for a generic specialization of the variables P, the zeros of f_{1}, …, f_{n} are isolated and nonsingular. We call generic fiber the fiber above the prime (0) ⊂ k[P_{1}, …, P_{m}] ([9]). Our aim is to give a description of this fiber, i.e. a parametric resolution of the system, with coefficients in k(P_{1}, …, P_{m}).

We propose a probabilistic algorithm that does so; the idea is to lift the generic fiber starting from a lucky specialization of the parameters. This lies in the continuity of the work of the TERA group, particularly of the preprint [15]. In fact, both papers make use of techniques that were developed in the papers [11], [12].

We first show the existence of a resolution with good degree properties, then exhibit an algorithm to obtain this resolution and give its complexity. We have implemented the algorithm in Magma [1], on the basis of the Kronecker package developed by G. Lecerf [13],[2]. Thus, a fair size of this paper is devoted to the results of our computations.

2 Hypotheses, notation and results

Let k be an effective field and P_{1}, …, P_{m}, X_{1}, …, X_{n} be indeterminates. Call K the field k(P_{1}, …, P_{m}). We consider some polynomials f_{1}, …, f_{n} in k[P_{1}, …, P_{m}, X_{1}, …, X_{n}] that generate an ideal I in k[P_{1}, …, P_{m}, X_{1}, …, X_{n}]; the associated variety in 𝔸m+n ; \overline{k} will be called V. We do not suppose that the variety V is in Noether position with respect to the variables P. Let π be the projection morphism:

\begin{aligned}π : & V & → & 𝔸m ; \overline{k} \\ & (p, x) & ↦ & p. \end{aligned}

We suppose that π is dominant and generically finite, let Δ be the Zariski-closed subset of 𝔸m ; \overline{k} above which π is not finite, and U the closure of V −π^{−1}(Δ); then U is our object of interest. We suppose that the ideal [I : I(Δ)^{∞}] is radical, which means that the generic fiber is reduced.

Using the definition of degree given in [14], we call

• d_{P,X} (resp. d_{X}) the maximum of the degrees (resp. in the variables X) of the input polynomials;

• D_{P,X} the degree of U, with D_{P,X} ≤ d_{P,X}^{n};

• D_{X} the degree of the generic fiber, with D_{X} ≤ d_{X}^{n};

• δ ≤ d_{X}^{n−1}, the maximum of the degrees of the intermediate varieties V (f_{1}(p, x), …, f_{i}(p, x)) ⊂ \overline{k}^{n}, where p is a generic parameter, for 1 ≤ i ≤ n − 1.

The solutions of a zero-dimensional system can be represented by means of an univariate representation. This representation bears the name Rational Univariate Representation in [19], or Geometric Resolution in [11],[12],[15],[13]. Our positive-dimensional problem is actually zero-dimensional over K, and we want a resolution (t, q, w) of the generic fiber with coefficients in K. This consists in:

• a separating element t of the extension K → K[U];

• its monic minimal polynomial q ∈ K[T] of degree D_{X};

• the parametrization w of the points of U, of the form q′(t)x_{i} = w_{i}(t), 1 ≤ i ≤ n. The w_{i} belong to K[T] and have degree strictly less than D_{X}.

The separating element can be chosen with coefficients in k. Indeed, if t is a separating element for some fiber that fulfills all conditions of genericity, then t separates the points of all fibers above an open subset of 𝔸m ; \overline{k}, as the generic fiber is reduced.

The polynomials f_{1}, …, f_{n} are given by a Straight-Line Program (SLP in the sequel) S. We will call

• L the size of S itself, and

• ℓ the size of the specialization of S on a generic parameter, which is less than L.

M(d) will be the cost of multiplying two univariate polynomials of degree d on the field k, in terms of operations in the field k. M(d) can be taken O(d log d), when the characteristic is zero or big enough ([8]). M(d, m) will be the cost of m-variate polynomial multiplication, truncating at precision d; this can be taken less than M((2d)^{m}) ([8]).

Theorem 1 We can compute the resolution of the generic fiber by a probabilistic algorithm in time

O(n(nℓ + n^{4})M(d_{X}δ)^{2} + (nL + n^{4})M(D_{X})M(4D_{P,X}, m)).

The probabilistic aspect is twofold. First we choose a random parameter, which must satisfy genericity conditions. We distinguish three cases, which have different consequences on the behavior of the algorithm (see 4.5).

Theorem 2 Let Ω be a subset of k of cardinality N. Let p = (P_{1}, …, P_{m}) ∈ Ω^{m}. Then

• the probability that the fiber above p is not finite is less than md_{P,X}^{n}/N;

• the probability that the fiber above p is not reduced is less than m(2(D_{P,X} − D_{X})(D_{X} − 1) + D_{X}(D_{X} − 1))/N;

• the probability that the fiber above p has not generic degree is less than m(D_{P,X} − D_{X})/N.

The other point is the resolution of a zero-dimensional system, where we use the probabilistic algorithm given in [13]. The probability of failure of a similar algorithm is given in [16], and is not developed here.

Finally, note that there are other solutions to obtain the parametric resolution. We now assume that the last variable X_{n} is a separating element for the generic fiber. This is by no means a restriction, for we can always get to this case through a change of variables with coefficients in k.

• A Gröbner basis in the polynomial ring k[P_{1}, …, P_{m}, X_{1}, …, X_{n}] for some elimination order gives more than the resolution of the generic fiber: it computes a resolution of all components of the variety V, so this output is often way too big.

• A Gröbner basis with lexicographic order X_{1} > ⋯ > X_{n} in the ring K[X_{1}, …, X_{n}] yields an univariate representation which closely resembles our previous definition. The difference lays in that the Gröbner basis’ parametrizations have much bigger coefficients, as they encode the inverse of q′ modulo q.

• A last possibility is to compute a Rational Univariate Representation from a Gröbner basis in K[X_{1}, …, X_{n}] with any order, using for instance the algorithms presented in [19].

We refer the reader to section 6.4, in which most of these possibilities are experimented on a particular example.

3 Existence of a resolution

We first show the existence of an eliminating polynomial for any variable, with coefficients in k[P_{1}, …, P_{m}] and suitable degree, then deduce the existence of a resolution. The next section will give an algorithm to compute this resolution.

Proposition 1 For each i in {1, …, n} there is a polynomial Q_{pol} ∈ k[P_{1}, …, P_{m}][X_{i}] that belongs to I(U). The total degree of Q_{pol} is less than D_{P,X}, its partial degree in X_{i} is less than D_{X}.

Proof. The proof closely follows that of [21, Proposition 1]. Consider the morphism

\begin{aligned}ϕ : & U & → & 𝔸m+1 ; \overline{k} \\ & (p, x) & ↦ & (p, x_{i}). \end{aligned}

We call W the closure of Im(ϕ), which is an hypersurface

of 𝔸m+1 ; \overline{k}.

We first give a bound on the degree of W. Let H be a line in 𝔸m+1 ; \overline{k} such that #(H ∩ W) equals deg W. Build

H′ ⊂ 𝔸m+n ; \overline{k} the cylinder H × \overline{k}^{n−1}. Then H′ ∩ U is the set {(p, x) ∈ U|(p, x_{i}) ∈ H ∩ W}. The construction of U implies that H′ ∩ U is finite, and so of cardinal less than D_{P,X}. It follows that #(H ∩ W) is also less than D_{P,X}, and finally, deg W ≤ D_{P,X}.

Let Q_{pol} be a square-free polynomial in \overline{k}[P_{1}, …, P_{m}, X_{i}] such that W = Q_{pol}^{−1}(0), so that deg Q_{pol} = deg W ≤ D_{P,X}. Notice that for (p, x) ∈ U, Q_{pol}(p, x_{i}) = 0. One can see Q_{pol} as a monovariate polynomial in X_{i} over \overline{k}[P_{1}, …, P_{m}]; we then have the following properties:

• The degree of Q_{pol} in X_{i} is the degree of the projection of a generic fiber of U on the X_{i}-axis, so it is less than the degree D_{X} of the generic fiber.

• The content of Q_{pol} is 1.

We now have to show that Q_{pol} actually belongs to k[P_{1}, …, P_{m}, X_{i}]. Following the notations in [21], we call B the algebra k[U], then B′ the product B ⊗_{k[P_{1},...,P_{m}]} K and finally B′′ the tensor product B′ ⊗_{k} \overline{k}, which is a finite dimensional vector-space of dimension D_{X} over K′ = \overline{k}(P_{1}, … P_{m}). We denote by M_{X_{i}} the endomorphism of multiplication by X_{i} in B′′, and by m_{X_{i}} its minimal polynomial. As the endomorphism M_{X_{i}} is actually defined over K, its minimal polynomial has coefficients in K. Consider Q_{rat} ∈ K′[X_{i}], the polynomial Q_{pol} divided by its leading coefficient in X_{i}. Then m_{X_{i}} divides Q_{rat}. As they have the same degree, these two polynomials are equal, so that Q_{rat} has coefficients in K. As the content of Q_{pol} in \overline{k}[P_{1}, …, P_{m}] is 1, Q_{pol} can be recovered as the numerator of Q_{rat}, and also has its coefficients in k. This concludes the proof. □

Proof of Theorem 2. All assumptions are consequences of Zippel-Schwartz’ lemma [23], which states that the number of zeros in a set of the form Ω^{m} of a m-variate polynomial of partial degree in every variable bounded by d is less than md|Ω|^{m−1}. Consequently, in each case, the proof amounts to enclosing bad parameters in an hypersurface of suitable degree.

• The non-finite fibers project on an hypersurface of degree less than d_{P,X}^{n}.

• The non-reduced fibers project on a hypersurface which is included in the discriminant of the eliminating polynomial of a generic separating element. It can be checked that the degree of this discriminant is less than 2(D_{P,X} − D_{X})(D_{X} − 1) + D_{X}(D_{X} − 1), using the properties of homogeneity and weighted homogeneity of the discriminant.

• The fiber has non-generic degree when the leading coefficient of the eliminating polynomial of a generic separating element, as defined in Proposition 1, cancels; this coefficient has degree less than D_{P,X} − D_{X}. □

Using well-known techniques ([17], [18]), one can recover a whole resolution of the generic fiber through the computation of the characteristic polynomial of a generic separating element. This is summarized in the following proposition:

Proposition 2 There exists a resolution of the generic fiber of the form

\begin{aligned}Q(T) & = & 0 \\ \frac{∂Q}{∂T} (T)X_{i} & = & W_{i}(T), & i = 1, …, n, \end{aligned}

where Q, W_{1}, …, W_{n} belong to K[T]. The polynomial Q is monic and has degree D_{X}; the W_{i} have degree strictly less than D_{X}. All the coefficients of these polynomials have numerator and denominator of degree less than D_{P,X}.

Proof. Consider the characteristic polynomial Q_{rat} of a generic primitive element T = ∑ λ_{i}X_{i}. Its coefficients have suitable degree in the parameters, according to the previous Proposition. Its derivative with respect to the variable λ_{i} yields the parametrization W_{i} we seek. Q_{rat} is polynomial in λ_{i}, so that differentiation with respect to λ_{i} does not increase the degree of the rational functions in the parameters P that appear. □

Note that this construction implies that the parametrizations and the eliminating polynomial are of similar degrees and heights.

4 Description of the algorithm

Our algorithm relies on a technique of specialization and “unspecialization” of the parameters, which is similar to the one in [15]. It is based on the computation of successive approximations of the resolution in k[[P_{1} − p_{1}, …, P_{m} − p_{m}]], for a parameter p = (p_{1}, …, p_{m}). This approximation process appears in the papers [11], [12], during the intersection process of the geometric resolution.

We use the following algorithmic tools:

• Newton Lifting([13]) Consider a polynomial system F_{1}, …, F_{n} in a ring A[X_{1}, …, X_{n}], and I an ideal of A. Suppose that the jacobian of this system is invertible on the points defined by F = 0, and that we have a geometric resolution of the system F = 0 modulo I. Then the Newton algorithm allows to compute a resolution modulo any power of the ideal I. Here, I will be of the form (P_{1} − p_{1}, …, P_{m} − p_{m}) in the ring k[P_{1}, …, P_{m}].

• Rational Reconstruction Two similar questions of rational reconstruction will be asked. First, we let A be the ring k[P_{1}, …, P_{m}], K its quotient field and I ⊂ A the ideal (P_{1} − p_{1}, …, P_{m} − p_{m}). Let now q be an element of K defined at (p_{1}, …, p_{m}), with numerator and denominator of degrees bounded by η. We want to recover q from its I-adic expansion (i.e. its Taylor expansion at (p_{1}, …, p_{m})).

Later on, we will ask the same question in the ring ℤ: we want to recover a rational number q = a/b from its reduction modulo some number p. This time, we call η the maximum of the heights of a and b, and the precision of the expansion will be the height of p. In both cases, one can recover q from its expansion at precision 2η using Padé approximants. This process will be called the rational reconstruction of q. When the expansion is not known with enough precision, the reconstruction can be impossible.

Here is the algorithm, which divides in two parts:

Compute_Resolution:
Input: the polynomials f=(f1,...,fn)
Output: a resolution (t,Q,W) of the generic fiber
#--------------------------------------#
p=(p1...pm)<-random parameter
(t,q,w)<-resolution of the fiber above p
Return Lift_Generic_Fiber(f,p,(t,q,w))
#--------------------------------------#

The routine Lift Generic Fiber is the following:

Lift_Generic_Fiber:
Input: the polynomials f=(f1,...,fn)
       a parameter p=(p1,...,pm)
       a resolution (t,q,w) of the fiber above p
Output: a resolution (t,Q,W) of the generic fiber
#--------------------------------------#
k<-1
While Not Finished do
    (t,q,w)<-Lift(f,p,(t,q,w))
         modulo (P1-p1,...,Pm-pm)^(2^k)
    Finished,(t,Q,W)<-StopCriterion(f,(t,q,w))
    k<-k+1
end While
Return (t,Q,W)
#--------------------------------------#

We enter the κ-th step of the main loop with a resolution of the fiber at precision 2^{κ−1}, and lift it to precision 2^{κ} via a Newton lifting step. The subroutine StopCriterion first tries to build a rational reconstruction of all the coefficients of the resolution, then tests the reconstruction on the input system. To this effect, we use a subroutine

Rational Reconstruction which outputs true and a recon-

struction of all the coefficients if it exists, false otherwise.

StopCriterion:
Input: the polynomials f=(f1,...,fn)
       (t,q,w): resolution with coefficients
          in k[[P1-p1,...,Pm-pm]]
Output: b: boolean
       (t,Q,W): resolution with coefficients
          in k(P1,...,Pm)
#--------------------------------------#
b,(Q,W)<-Rational_Reconstruction(q,w)
If b and (f(P,W/Q')=0 mod Q)
    Return true,(t,Q,W)
else
    Return false
#--------------------------------------#

This subroutine can be replaced by a probabilistic one, where one tests the resolution specialized on a new random parameter. The probability of error could once again be bounded using Zippel-Schwartz’ lemma.

We now show that the algorithm actually finishes.

Proposition 3 From a generic choice of (p_{1}, …, p_{m}), the algorithm computes a resolution of the generic fiber of the variety U in at most ⌈log_{2}(D_{P,X})⌉ + 1 steps.

Proof. Let p = (p_{1}, …, p_{m}) be a random parameter and t be a separating element for the fiber above p. Then t is a separating element in an open neighborhood of p, so there is a unique resolution of the generic fiber with t as separating element, and a unique resolution of the generic fiber with coefficients in k[[P_{1} − p_{1}, …, P_{m} − p_{m}]] for this separating element; these coefficients coincide with the power series expansions of those of the aforementioned resolution around p.

The fiber above p is reduced, so we can apply the Newton algorithm to lift the specialized resolution to any precision. As this fiber has generic degree D_{X}, κ steps of the main loop of the algorithm yield the resolution of the generic fiber with coefficients in k[[P_{1} − p_{1}, …, P_{m} − p_{m}]] at precision 2^{κ}. We know that these series are expansions of rational functions, with degrees of numerators and denominators bounded by D_{P,X}. Thus it is enough to compute the series at precision 2D_{P,X} to recover the rational functions, so that the resolution is computed in at most ⌈log_{2}(D_{P,X})⌉ + 1 steps. □

We now give the complexity of each step of the algorithm. All complexity results are stated in terms of operations in the field k.

4.1 The specialized resolution

The first part of the algorithm consists in the resolution of the fiber above a random choice p = (p_{1}, …, p_{m}) of the parameters. We have to solve a zero-dimensional system, producing as output a geometric resolution of the chosen fiber. As usual, the output (t, q, w) consists in a separating element t of the points of the fiber above p, its minimal polynomial q ∈ k[T] of degree D_{X} and the parametrization of the coordinates of the form q′(t)x_{i} = w_{i}(t), 1 ≤ i ≤ n.

This resolution can be computed in time polynomial in the previously defined quantities, following the algorithms in [11]. We state here the theorem given in [13]:

Proposition 4 The resolution of the fiber above p can be computed in time

O(n(nℓ + n^{4})M(d_{X}δ)^{2})

by a probabilistic algorithm.

In some applications, the degree of a generic fiber is much less than the degree D_{P,X} of the variety U. This first part of the algorithm is then scheduled to be the less timeconsuming.

Finally, any other procedure that outputs a geometric resolution is as well suited to this first step of the algorithm.

4.2 Lifting the resolution

The second step is an iterative unspecialization process that relies on the following symbolic Newton lifting procedure:

                                                        Lift:
                                                        Input: the polynomials f=(f1,...,fn)
                                                               the parameter p=(p1,...,pm)
                                                               a resolution (t,q,w) of the fiber
                                                                   f(P,.)=0 modulo (P1-p1,...,Pm-pm)^(2^k)
                                                        Hypothesis: the jacobian of f(p,X)
                                                                    is invertible on the fiber
                                                        Output: a resolution (t,q,w) of the fiber
                                                                   f(P,.)=0 modulo (P1-p1,...,Pm-pm)^(2^(k+1))

Proposition 5 Given a parameter p = (p_{1}, …, p_{m}), and a resolution of the fiber of U above p, the κ-th step of Newton’s lifting algorithm takes time

O((nL + n^{4})M(D_{X})M(2^{κ}, m)).

Proof. The paper [13] gives the complexity of the Newton procedure in the more general context of lifting modulo the powers of any ideal I of the coefficients ring. The κ-th Newton step costs O((nL + n^{4})M(D_{X})a(2^{κ})), where a(h) is the cost of arithmetic modulo I^{h}. In our setting, a(h) is the cost of m-variate series multiplication at precision h, i.e. M(h, m), so that our κ-th step costs O((nL + n^{4})M(D_{X})M(2^{κ}, m)). □

4.3 Rational Reconstruction

We want to obtain a resolution with coefficients in k(P_{1}, …, P_{m}) with numerators and denominators of degrees less than 2^{κ−1}, from a resolution modulo (P_{1} − p_{1}, …, P_{m} − p_{m})^{2^{κ}}. To this effect, we compute Padé approximants of all the coefficients, if they exist. We first give the cost of the reconstruction of a single coefficient in the mono- and multivariate cases, then take all these coefficients into account.

• In the monovariate case, a rational approximant of a series known at precision 2^{κ} can be computed in O(κM(2^{κ})) operations in the base ring, see for instance [8].

• In the multivariate case, numerous versions of Padé approximants coexist. We propose here a straightforward continuation of the monovariate case. Consider a truncated series s ∈ k[[P_{1}, …, P_{m}]] - we have translated (p_{1}, …, p_{m}) to zero to simplify things a bit. We get back to a monovariate situation by introducing s̃ = s(P_{1}t, …, P_{m}t), for a new variable t. Now s̃ can be seen as element of k(P_{1}, …, P_{m})[[t]] of precision 2^{κ}, and we can compute an approximant ã of s̃ with numerator and denominator of degree less than 2^{κ−1} in t. Substituting t = 1 in ã gives a rational function whose power series expansion coincide with s. We perform computations on homogeneous coefficients, truncated at precision 2^{κ−1}, so if there is a solution for the multivariate approximation problem with correct degree, we can find it in O(M(2^{κ})κM(2^{κ−1}, m − 1)).

Proposition 6 The total cost of rational reconstruction at step κ is

O(nD_{X}κM(2^{κ})M(2^{κ−1}, m − 1)).

4.4 Putting everything together

We now show the bound on the whole complexity of our algorithm. It is the sum of the costs of the zero-dimensional resolution plus the sum of all steps of the lifting. The total number of lifting steps is bounded by ⌈log_{2}(D_{P,X})⌉ + 1 and the sum of the costs of all steps can be dominated by twice the cost of the last step, if we use the inequalities M(x, m) ≤ M(2x, m)/2. Furthermore, the cost of rational reconstruction is dominated by the cost of the Newton iteration. All this yields the following upper bound:

Proposition 7 The total complexity of the algorithm is

O(n(nℓ + n^{4})M(d_{X}δ)^{2} + (nL + n^{4})M(D_{X})M(4D_{P,X}, m)).

4.5 Troubleshooting

Finally, we detail what happens when the random parameter is not generic:

• Infinite and non-reduced fibers are detected during the zero-dimensional solving process.

• Fibers with non-generic degree are likely to make the reconstruction impossible, so we detect that at the end of the algorithm. Still, it is possible that a whole component disappears, so that we compute a resolution of all remaining ones.

5 Practical strategies

This section presents possible improvements of the main algorithm. These features do not drastically modify the complexity class of the problem, but have some practical impact on the examples.

5.1 Modular arithmetic

When a resolution is requested on the rational field, we perform the computations modulo some powers of prime numbers. Note that all computations in the previous algorithm can be done modulo p^{h} for all primes p but a finite number. So, we first compute the generic resolution modulo powers of different primes, then apply the Chinese Remainder Theorem followed by a rational reconstruction on ℚ on each coefficient of all the rational functions that appear.

We denote by A cat B the concatenation of the lists A and B. Then the whole algorithm becomes:

Compute_Resolution_Rational :
Input: the polynomials f=(f1,...,fn)
Output: a resolution (t,Q,W) with rational
coefficients of the generic fiber
#--------------------------------------#
p=(p1,...,pm)<-random parameter
(t,q,w)<-resolution with rational coefficients
of the fiber above p
h<-height of (q,w)
pn<-random prime number
While Not Finished do
Lr<-Lr cat
[Lift_Generic_Fiber(f,p,(t,q,w)) mod pn^h]
Lp<-Lp cat [pn]

     pn<-NextPrime(pn)
     Finished,(t,Q,W)<-StopCriterion(Lp,Lr)
end While
Return (t,Q,W)
#--------------------------------------#

This time, StopCriterion tests if it is possible to build a rational reconstruction on ℚ of all the coefficients of the rational functions that appear in the resolutions, and if the reconstruction satisfies the equations f_{1}, …, f_{n}. If it is the case, it returns true and the reconstruction; otherwise it returns false.

We compute all resolutions with height of order the height of the specialized fiber. Indeed, it is likely that the process of specialization has increased the height of the resolution, so that the generic fiber should have smaller coefficients than the specialized fiber. Then all chances are that a single prime number is enough to build the resolution with rational coefficients. An analysis of the probability of choosing bad reduction primes can be done, using arithmetic Bézout inequalities.

It is noteworthy that we could also use a lifting process to obtain a resolution on the rational field, but our implementation uses the former interpolation technique for practical reasons.

5.2 Factorization

Suppose that the minimal polynomial of the specialized fiber splits into i irreducible factors of degrees D_{X1}, …, D_{Xi}, so that there exist as many resolutions that represent the irreducible components of this fiber. Then, even when the generic fiber is irreducible, it is possible to take profit of this factorization: the lifting is done on all smaller resolutions, which are combined before each call to StopCriterion. Thus the term M(D_{X}) = M(∑ D_{Xi}) in the complexity of the lifting step is replaced by ∑ M(D_{Xi}) - recall that for a naive multiplication algorithm, M(d) = O(d^{2}).

Furthermore, the generic fiber may actually factor in several components; suppose this is the case, and that we are only interested in some of those components. The specialized fiber also factors, into possibly more factors. If it is possible to identify the “trace” of the generic factorization on the specialized fiber, then the lifting can be done on only the components of interest. See example 6.3 for an application.

5.3 Computing outside a given hypersurface

Suppose we want the points of U lying outside a given hypersurface H. It suffices to remove the points of H that intersect the specialized fiber, and perform the lifting on what remains. The algorithm for geometric resolution given in [13] allows to remove these points quite easily; other zerodimensional solving processes also give this possibility, using for instance Rabinowicz’ trick.

For instance, if the generic fiber is not reduced, it is possible to obtain all its simple points. It suffices to remove the multiple points of the specialized fiber. This is the case in example 6.1.

6 Examples

This last section gathers examples of various sources where we have applied our techniques.

The algorithm is implemented in Magma, on the basis of the Kronecker package developed by G. Lecerf [2], which gives the resolution of the specialized system and also most functionalities necessary for our algorithm. Consequently, these computations will be referred to as “Kronecker” in the time-tables.

We also tried Gröbner Bases computations, as in each example we could find a separating variable. We used Magma, which allows for such computations on many fields, including rational function fields.

Our computations were done on the machines of the UMS Medicis [3]. The notation EV6 stands for a Compaq Alpha EV6 XP/1000 500 Mhz processor, and EV56 for Digital PWS/500au - Alpha EV56 500 Mhz processor.

All systems are not displayed for lack of space. These equations are available at request.

6.1 Number of points of a Jacobian

This example relates computations that were performed with P. Gaudry and R. Harley for their paper [10]. The framework is the determination of the number of points of the jacobian of a given hyperelliptic curve over a (big) finite field k. Following Schoof’s algorithm for elliptic curves, the algorithm is based on explicit computation of points of p^{i}- torsion of this group, for various primes p. This part relates the computations for the case p = 2.

There is a polynomial system that gives a point P_{i+1} of 2^{i+1}-torsion from the knowledge of a point P_{i} of 2^{i}-torsion. This polynomial system has 4 equations in 6 variables. The variables split in the 4 “coordinates” of the point P_{i+1} plus 2 parameters that are functions of P_{i}.

Given a point P_{1}, the resolution of this system gives the point P_{2}, which is in turn fed to a similar system... As i increases, these points live in extensions of k of increasing degrees, so that the resolution of the instantiated systems gets longer and longer. It is preferable to compute for once the parametrized resolution, then to specialize it when needed.

Here are the first lines of the SLP that encodes the system, which was obtained using the optimize function in Maple [4]. The coefficients of this program belong to the finite field k = 𝔽_{p}, where p is the first prime greater than 10^{19}. The first two variables are the parameters.

Vars:=[c1,c2,u1,u2,v0,v1];
Syst:=function(Vars)
 f1:=3141592653589793238;
 f2:=4626433832795028841;
 f3:=9716939937510582097;
 f4:=4944592307816406286;
 f5:=2089986280348253421;
 t1 := u1^2;  t2 := t1^2;  t3 := t1*u1;
 t4 := f1*t3; t5 := 3*u2; t8 := f1*u2;
 t10 := v0^2; t13 := u2^2; t14 := f2*u2;
 t20 := 2*t13; t23 := f1*t13; t26 := v1^2;
 .........
end function;
Cost : l=251 multiplications         209 additions
       L=253 multiplications         209 additions

The following table gives the time necessary to compute the generic fiber. In this process, we somehow changed the algorithm we presented. Instead of lifting both parameters at once in k[[C_{1} − c_{1}, C_{2} − c_{2}]] we lift them in directions of the form C_{1} = c_{1} + C, C_{2} = c_{2} + iC, for various i, thus computing in k[[C]]; then we interpolate the results. This brings back the problem to monovariate computations, at the theoretical loss of a factor two. Still, multivariate computations required unexplainable amounts of memory in Magma, so the monovariate computations was preferred.

Kronecker EV6 Gröbner Magma EV6
160 min.  ∞

• The column “Kronecker” indicates the total runningtime of our algorithm.

• The column “Magma” stands for the computation of a lexicographic Gröbner basis with coefficients in k(c_{1}, c_{2}). Computations were stopped after 3 days of user-time, and used more than 900 Mo of memory.

6.2 The system “Hawes4”

This example was borrowed from J.C. Faugère’s homepage [5] (see also the database Polydata [6]). This system has 6 equations plus two inequations in the 8 variables a, b, c, x, y_{1}, z_{1}, y_{2}, z_{2}. Thanks to J.C. Faugère, we know that the original question was “We wish to eliminate x, y_{1}, z_{1}, y_{2}, z_{2} from the system ignoring the trivial solutions of y_{1} = y_{2} and z_{1} = z_{2}”. This amounts to computing the minimal polynomial of say c (which is not separating) with coefficients in ℚ(a, b). Once we have the generic fiber, this is easy. A straightforward way to do so is the computation of a resultant; a more efficient solution is presented in [22]. Here is this system, under a not expanded and “easy-to-

evaluate” form, followed by the cost of the SLP we used.

                                                        Vars:=[a,b,c,x,y1,z1,y2,z2];
                                                        Syst:=[
                                                         3*z1^2+y1^2+b,
                                                         y2^2+3*z2^2+b,
                                                         2*c*y1+x+2*z1*y1+5*y1^4+3*a*y1^2,
                                                         2*c*y2+x+5*y2^4+3*a*y2^2+2*y2*z2,
                                                         -c*y2^2+c*y1^2-x*y2+x*y1+z1^3+z1*y1^2+b*z1-y2^5
                                                         -a*y2^3-y2^2*z2-z2^3-b*z2+y1^5+a*y1^3,
                                                          -6*c*y1^2*z1*(3*x*z2-45*y2^4*z2+y2^3-3*y2*z2^2-b*y2)
                                                         +6*c*y2^2*z2*(3*x*z1-3*z1^2*y1-45*z1*y1^4+y1^3-b*y1)
                                                         -2*(3*z1^2+30*z1*y1^3+9*a*z1*y1-y1^2)*
                                                         (-9*a*y2^2*z2+3*x*z2-45*y2^4*z2+y2^3-3*y2*z2^2-b*y2)+
                                                         2*(30*y2^3*z2-y2^2+9*a*y2*z2+3*z2^2)*
                                                         (-9*a*z1*y1^2+3*x*z1-3*z1^2*y1-45*z1*y1^4+y1^3-b*y1)
                                                        ];
                                                        Ineq:=[y1-z1,y2-z2];
                                                        Cost : l=66 multiplications         50 additions
                                                               L=76 multiplications         52 additions

This example is relevant of the behavior of our algorithms with respect to the representation of the input polynomial. The last polynomial of the system was given under the expanded form:

                                                        P6:=-2*y2^2*y1^3+6*z2^2*y1^3+6*z1^2*b*y2-2*y1^2*b*y2
                                                        +18*z2^2*x*z1+18*z1^2*y2*z2^2+540*z1*y1^3*a*y2^2*z2+
                                                        +...

Cost : l=161 multiplications         52 additions
L=191 multiplications         57 additions

The compact representation was found manually. It could be argued that this part is not automatizable, but it is clear that the system was originally written in a compact way. This has an important consequence on the computation times, see the comments hereafter. As in the previous example, we did not lift both parameters at once.

Kronecker EV6 Gröbner Magma EV6 FGb PPro
15 + 4.5 h  ∞  12 h

• The column “Kronecker” indicates the total runningtime of our algorithm. The timing splits in the computation of first the generic fiber, then of the minimal polynomial of c. This computation using the expanded form of the system takes app. 50 h.

• The column “Magma” stands for the computation of a lexicographic Gröbner basis with coefficients in ℚ(a, b). As in the previous example, they were stopped after several days of user-time, and used more than 650 Mo.

• The column “FGb” gives the computation time as shown on J.C. Faugère’s homepage [5], on a PC Pentium Pro 200 Mhz. This directly gives the minimal polynomial of c via the computation of an elimination ideal.

6.3 Automatic geometric theorem proving

Most questions of automatic theorem proving fit naturally in our framework, as the indeterminates split in two groups of parameters and variables. The question can then often be brought to test whether some polynomial functions cancel on the generic fiber of a variety of positive dimension, disregarding degenerate cases. We propose a straightforward algorithm: compute the generic fiber, then test whether the conclusions hold. This will be illustrated by a simple theorem of plane geometry.

This geometric proposition was pointed out to us by G. Lecerf. Given a vector AB, we define a function f_{AB}(M) as the abscissa of M on the line (MO) such that O ∈ (AB) and ⟨MO, AB⟩ is orthogonal, direct.

Proposition 8 Let ABCD be a plane quadrilateral. If the sum f_{AB}(M) + f_{BC}(M) + f_{CD}(M) + f_{DA}(M) is a constant function of M, then ABCD is a parallelogram.

Our description of this problem uses 8 variables describing the quadrilateral ABCD, that are linked by 4 equations plus 2 inequations. These equations define a 4-dimensional variety, taking the first four variables as parameters, and the conclusion amounts to the cancellation of 2 polynomial functions on the generic fiber of this variety. Here are the times necessary to compute the generic fiber:

Kronecker EV6 Gröbner Magma EV6
37 s.  3 h.

One then verifies that the quadrilateral ABCD in the generic fiber actually is a parallelogram. This example shows the good behavior of our method:

• The specialized fiber factors in two components: one is a well-behaved parallelogram ABCD, the other one is of the form ABCB. One verifies that the generic fiber comprises this uninteresting component, so that the lifting can be done only on the interesting, wellbehaved component.

• Some degenerate cases naturally appear as denominators of the coefficients in the resolution. One needs not bother about them in advance, as they are obtained through the process of rational reconstruction during the resolution.

These nice features will appear on other examples. We have tried our method on other problems, such as Pappus’ theorem, with the same success, still keeping in mind that our algorithm is probabilistic.

6.4 Deformations of singular hypersurfaces

The theoretical background of this section can be found in the paper by F. Rouillier, M.F. Roy and M. Safey [20] and references therein. These authors address the problem of finding one point in each connected component of a hypersurface H = P^{−1}(0) ⊂ ℝ^{n}. This is achieved by considering the critical points on H of the function d_{A}(M) = ∥AM∥^{2}, for some point A, which is done by computing the (complex) zero-set of the system

P = 0, gradP//AM.

We suppose we have a generic enough point A, so that this system is zero-dimensional if and only if the hypersurface H has a finite number of singularities. When H has an infinite number of singularities, one introduces an infinitesimal parameter ϵ and studies the critical points on the level sets P = ϵ, i.e. the system

P = ϵ, gradP//AM.

Seeing this system as parametrized by ϵ, Sard’s theorem implies that the generic fiber is finite, so that our algorithm gives a parametric resolution of the generic fiber above ϵ.

We experimented on the following example, which is borrowed from M. Safey’s homepage [7]. The point of this example is not so much in the computation times as in the sizes of the various outputs. Here is the system, in a conventional, dense representation. The following cost is that of the SLP that we used to evaluate it.

Vars:=[e,t2,t3,t4];
Syst:=[t2^6+t2^4*t3^2+t2^4*t4^2-t2^2*t3^4-2*t2^2*
t3^2*t4^2-t2^2*t4^4+t3^6+3*t3^4*t4^2+3*t3^2*t4^4
+t4^6-e,
4*t2^5*t3-12*t2^5+2*t2^4*t3+8*t2^3*t3^3-8*t2^3*t3^2+
8*t2^3*t3*t4^2-8*t2^3*t4^2-4*t2^2*t3^3-4*t2^2*t3*t4^2
-8*t2*t3^5+4*t2*t3^4-16*t2*t3^3*t4^2+8*t2*t3^2*t4^2
-8*t2*t3*t4^4+4*t2*t4^4+6*t3^5+12*t3^3*t4^2+6*t3*t4^4,
4*t2^5*t4-18*t2^5+2*t2^4*t4+8*t2^3*t3^2*t4-12*t2^3
*t3^2+8*t2^3*t4^3-12*t2^3*t4^2-4*t2^2*t3^2*t4-4*t2^2
*t4^3-8*t2*t3^4*t4+6*t2*t3^4-16*t2*t3^2*t4^3+12*t2
*t3^2*t4^2-8*t2*t4^5+6*t2*t4^4+6*t3^4*t4+12*t3^2
*t4^3+6*t4^5];
Cost : l=25 multiplications          12 additions
       L=25 multiplications          13 additions

We experimented various approaches on this example. The following table gives the running times and the sizes of the outputs.

Kronecker Kronecker+Gröbner Gröbner 1 Gröbner 2
Magma  Magma
EV56  EV56  EV56  EV56
170 s.  52 s.  119 s.  364 s.
5 Ko  5 Ko  351 Ko  27579 Ko

• The first two columns give the running time and size of the output of our package. In the first case, the zerodimensional instantiation was solved using Kronecker, whereas in the second case it was solved using Magma’s built-in Gröbner facilities.

• The column Gröbner 1 stands for the computation of a lexicographic Gröbner basis in ℚ(ϵ)[t_{3}, t_{4}, t_{2}] in Magma (t_{2} is a separating element). As expected, the size of the output is much more important, due to the size of the coefficients in the parametrization. One can check that these parametrizations multiplied by q′ modulo q yield the much smaller representation we propose.

• Finally, the column Gröbner 2 is the computation of a Gröbner basis of the system in ℚ[t_{4}, t_{3}, t_{2}, ϵ]. It could be argued that things could be better if we added some condition ϵ ≠ 0, and this objection is sound in this particular example. We only show this as an example of what could happen if the degeneracy conditions were not known a priori.

References

[1] http://www.maths.usyd.edu.au:8000/u/magma/.

[2] http://kronecker.medicis.polytechnique.fr/.

[3] http://www.medicis.polytechnique.fr/.

[4] http://www.maplesoft.com.

[5] http://posso.lip6.fr/˜ jcf/.

[6] http://www.polydata.medicis.polytechnique.fr/.

[7] http://www-calfor.lip6.fr/˜ safey/.

[8] Bini, D., and Pan, V. Polynomial and matrix computations. Progress in theoretical computer science. Birkhäuser, 1994.

[9] Eisenbud, D. Commutative Algebra with a view toward Algebraic Geometry. Springer, 1995.

[10] Gaudry, P., and Harley, R. Counting points on hyperelliptic curves over finite fields. Submitted to ANTS IV, 2000.

[11] Giusti, M., Haegele, K., Heintz, J., Montana, J., Morais, J., and Pardo, L. Lower bounds for diophantine approximations. In Journal of Pure and App. Algebra (Proceedings of MEGA’96) (1997), vol. 117.

[12] Giusti, M., Heintz, J., Morais, J., Morgenstern, J., and Pardo, L. Straight-line programs in geometric elimination theory. Journal of Pure and App. Algebra 124 (1998).

[13] Giusti, M., Lecerf, G., and Salvy, B. A Gröbner free alternative for polynomial system solving. To appear in Proc. FOCM’99.

[14] Heintz, J. Definability and fast quantifier elimination in algebraically closed fields. Theor. Comput. Sci. 24, 3 (1983).

[15] Heintz, J., Krick, T., Puddu, S., Sabia, J., and Waissbein, A. Deformation techniques for efficient polynomial equation solving. Submitted.

[16] Heintz, J., Matera, G., and Waissbein, A. On the time-space complexity of geometric elimination procedur es. Manuscript of Universidad Favaloro, Buenos Aires, Argentina, 1999.

[17] Macaulay, F. S. The Algebraic Theory of Modular Systems. Cambridge University Press, 1916.

[18] Renegar, J. On the computational complexity and geometry of the first-order theory of the reals. Part I. Journal of Symbolic Computation 13, 3 (1992).

[19] Rouillier, F. Solving zero-dimensional systems through the rational univariate representation. AAECC, 9 (1999).

[20] Rouillier, F., Roy, M.-F., and Safey, M. Finding at least one point in each connected component of a real algebraic set defined by a single equation. To appear in Journal of Complexity.

[21] Sabia, J., and Solernó, P. Bounds for traces in complete intersections and degrees in the Nullstellensatz. Applicable Algebra in Engineering Communications and Computing (AAECC Journal) 6 (1996).

[22] Shoup, V. Efficient computation of minimal polynomials in algebraic extensions of finite fields. In Proc. ISSAC’99 (1999).

[23] Zippel, R. Probabilistic algorithms for sparse polynomials. In Proceedings EUROSAM’ 79 (1979), no. 72 in LNCS, Springer.