Computing Parametric Geometric Resolutions

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

September 5, 2002

Abstract

Given a polynomial system of n equations in n unknowns that depends on some parameters, we define the notion of parametric geometric resolution as a means to represent some generic solutions in terms of the parameters.

The coefficients of this resolution are rational functions of the parameters; we first show that their degree is bounded by the Bézout number d^{n}, where d is a bound on the degrees of the input system. Then we present a probabilistic algorithm to compute a parametric resolution. Its complexity is polynomial in the size of the output and in the complexity of evaluation of the input system. The probability of success is controlled by a quantity polynomial in the Bézout number.

We present several applications of this process, notably to computations in the Jacobian of hyperelliptic curves and to questions of real geometry.

AMS Subject classification: 14Q20, 14Y05, 13P99, 68W30. Keywords: polynomial systems with parameters, complexity, theory of elimination, symbolic Newton operator.

1 Introduction

A variety of real-life problems can be modelized by polynomial systems involving free variables, or parameters. Even if a specialization of such a system may be easy to solve, a description of the generic solutions is typically hard to grasp. Still, the motivations for computing a generic description are numerous; we present various examples below. Our goal here is then to present an efficient, ready-to-implement, elimination procedure, adapted to such situations.

An introductory example. We can describe the main features of our approach on a simple example. Consider the following version of a serial robot, inspired by [30]: the robot has two segments of length 1, and is built as follows:

(x, y) θ_{1} θ_{2}

Figure 1: A simple robot arm

The coordinates (x, y) are thought as parameters, from which we want to recover the angles θ_{1}, θ_{2}. For i = 1, 2, take c_{i} = cos θ_{i} and s_{i} = sin θ_{i}. Then these values are related by the following polynomial system:

\{ \begin{matrix}c_{1} + c_{2} & = & x, \\ s_{1} + s_{2} & = & y, \\ c_{1}^{2} + s_{1}^{2} & = & 1, \\ c_{2}^{2} + s_{2}^{2} & = & 1. \end{matrix}

Our goal is to describe the solutions (c_{1}, s_{1}, c_{2}, s_{2}) in terms of (x, y) by the following kind of representation:

s_{2}^{2} − ys_{2} + \frac{1}{4} \frac{x^{4} + 2x^{2}y^{2} − 4x^{2} + y^{4}}{x^{2} + y^{2}} = 0 and \{ \begin{matrix}c_{1} & = & \frac{y}{x}s_{2} + \frac{1}{2} \frac{x^{2}−y^{2}}{x} , \\ c_{2} & = & −\frac{y}{x}s_{2} + \frac{1}{2} \frac{x^{2}+y^{2}}{x} , \\ s_{1} & = & −s_{2} + y. \end{matrix}

This representation gives a description of the solutions (c_{1}, s_{1}, c_{2}, s_{2}) for generic parameters values (x, y). Indeed, given any value of (x, y) in ℝ^{2} that does not cancel the denominators, s_{2} becomes the solution of a second-degree equation, with coefficients in ℝ. Then the values of (c_{1}, s_{1}, c_{2}) are given as functions of s_{2}.

The general case. We now proceed to describe the above process in greater generality. We suppose that we are given a polynomial system f = (f_{1}, …, f_{n}) in k[P_{1}, …, P_{m}, X_{1}, …, X_{n}], where k is any effective field. The variables P = (P_{1}, …, P_{m}) are thought as parameters, the variables X = (X_{1}, …, X_{n}) as unknowns, and we wish to compute a description of the solutions in terms of the parameters.

We restrict our study to the solutions that do not cancel the Jacobian determinant of the system f with respect to the variables X. These solutions will be called simple solutions.

We then define the notion of parametric geometric resolution as a description of these simple solutions by means of primitive element techniques. More precisely, the solutions will be described through the following encoding:

Q_{u}(u) = 0, { Q′_{u}(u)x_{1} = V_{1}(u), ; ⋮ ; Q′_{u}(u)x_{n} = V_{n}(u),

where

• x_{1}, …, x_{n} are the images of X_{1}, …, X_{n} in a suitable quotient algebra,

• u is a linear form in x_{1}, …, x_{n},

• Q_{u}, V_{1}, …, V_{n} are univariate polynomials whose coefficients are rational functions in P_{1}, …, P_{m}.

The polynomial Q_{u} is the minimal polynomial of the primitive element u. The polynomials V_{1}, …, V_{n} are parametrizations that give the values of the unknowns X in terms of the roots of the polynomial Q_{u}. This generalizes the representation of the solutions given in the introductory example.

We keep in mind that all coefficients of Q_{u}, V_{1}, …, V_{n} are rational functions in the parameters. In the general case, the introduction of the factor Q′_{u}(u) in the parametrizations will assure good degree properties for these coefficients.

Just as in the example, we also note that this representation gives a description of the generic solutions of the parametric system. Indeed, since the coefficients of Q_{u}, V_{1}, …, V_{n} are rational functions, we must avoid the parameter values that cancel one of their denominators.

Fields of application. Before presenting the content of this article, we describe some applications of these techniques. The last two examples cover some situations that are a priori non-parametric: yet, they can usefully be brought to our setting, through an appropriate shift of point of view.

• Specialization. A parametric resolution describes generic solutions. As such, it can be specialized on an open subset of the parameter space, where the denominators of the coefficients of the polynomials Q_{u}, V_{1}, …, V_{n} do not vanish. This enables to solve generic specializations of the input system f without further computations, for instance in software environments which do not provide elimination facilities.

We illustrate this possibility in Section 6 by the example of halving in the Jacobian of a genus 2 curve: dividing by 2 in a Jacobian amounts to solve a parametric polynomial system, whose parameters are the coordinates of the dividend. A parametric resolution of such a system was used in the genus 2 point-counting record of [23].

• Real geometry. As a particular case of parametric system, we may consider systems with infinitesimal coefficients, seeing the infinitesimals as parameters. Solving such systems is a cornerstone of many algorithms in real algebraic geometry [36, 37, 51]. These algorithms often require to study the limits of the solutions when the infinitesimals go to zero, for which a parametric resolution is well-suited.

As an example, we will show how to compute a family of critical points on deformations of a singular real hypersurface. The computation of such critical points is used as a subroutine in the algorithm of [51], which aims at computing one point on each connected component on a real hypersurface. Our implementation enabled a first comparison between several approaches for this question, applied to a practical interpolation problem in [52].

• Computing eliminating polynomials. Finally, many elimination questions can be reduced to our setting, and formulated as the computation of a suitable eliminating form for a polynomial function of the unknowns X_{1}, …, X_{n}. Such computations are straightforward once a parametric resolution is known.

This possibility is illustrated by an application coming from invariant theory [7]: the classification of a particular orbit space requires to compute the relation between 3 rational functions of 2 variables. We refer to Section 6 for the treatment of this question using our parametric formalism.

Overview of the article. The goal of this article is threefold. First, we estimate the complexity of a parametric resolution: we give bounds on the degrees of the polynomials that appear in such a representation. Next, we present a probabilistic algorithm for computing a parametric resolution, and estimate in a precise manner its complexity and probability of success. Finally, we present the applications mentioned above in greater detail.

• Complexity bounds. The first part of this article is devoted to prove that all coefficients that appear in a parametric resolution have degree bounded by an intrinsic geometric quantity, which itself is bounded by the Bézout number of the input system. This result is the continuity of [53, 27] and notably [34], and improves on the following important aspect.

The results obtained in the above references require the zero-set of the defining system f = (f_{1}, …, f_{n}) to be in Noether position with respect to the variables P_{1}, …, P_{m}. This condition implies that the number of solutions of all specializations of the system f is constant, if counted with multiplicities. In particular, this excludes all systems whose parametric resolution comprises denominators in the parameters.

Whereas this condition is not a limitation in the context of [53, 27, 34], it is a severe restriction for our situation. First, the validity of this condition is not easily tested for. Furthermore, many applications, for instance all those mentioned in the above paragraphs, do not fit into this setting. Our simple introductory example already gave a hint that the presence of denominators should be expected in many practical cases.

Treating the general case then requires a further geometric study, to control the complexity of the denominators in the coefficients. Correspondingly, if we want to follow the philosophy of elimination used in [27, 34], new algorithmic tools must be used.

• Algorithms. We propose an algorithm for a computing a parametric resolution, with complexity polynomial in the size of the output. This algorithm was implemented, and its practical behavior reflects its good theoretical complexity, enabling us to treat otherwise out-of-reach systems. The algorithm is probabilistic, and we have explicit estimates on the probability of success.

Here is a sketch of the method: the generic solutions will be obtained by successive approximations using a formal Newton operator.

The underlying paradigm is that solving a polynomial system over the base field k is a well-solved task. Thus, as input, we suppose that we are given a generic point (p_{1}, …, p_{m}) in the parameter space, and a description of the solutions of the system f specialized at (p_{1}, …, p_{m}).

This description has the form

q_{u}(u) = 0, { q′_{u}(u)x_{1} = v_{1}(u), ; ⋮ ; q′_{u}(u)x_{n} = v_{n}(u),

where q_{u}, v_{1}, …, v_{n} are univariate polynomials with coefficients in k, x_{1}, …, x_{n} are the algebraic variables and u is a linear form in x_{1}, …, x_{n}. Then there exists a parametric resolution composed of polynomials Q_{u}, V_{1}, …, V_{n}, such that q_{u}, v_{1}, …, v_{n} are the polynomials Q_{u}, V_{1}, …, V_{n}, with coefficients specialized at (p_{1}, …, p_{m}).

Our algorithm consists in lifting the dependency of the polynomials q_{u}, v_{1}, …, v_{n} in the parameters P_{1}, …, P_{m}; that is, we “unspecialize” the parameters. To this effect, we apply a formal Newton operator to [q_{u}, v_{1}, …, v_{n}]. This produces a sequence of polynomials whose coefficients are the successive power series expansions of the coefficients of [Q_{u}, V_{1}, …, V_{n}] at (p_{1}, …, p_{m}).

The idea of applying lifting techniques to solve polynomial systems can be traced back to the articles of Trinks [63] and Winkler [66]. It also underlies much of the recent work of the TERA group [27, 26, 25, 28, 34, 35], where the use of the Straight-Line Program encoding was the key to algorithms with good complexity.

As in the above references, we suppose that the input system is given by a Straight-Line Program. Our first contribution is then a generalization of the lifting operator of [28]. The computation sequence we propose becomes more lucid, the key point being to extend this operator to a wider class of representations, triangular representations.

Then, as in many algorithms relying on lifting techniques, the final step requires to go from a local description to a global one. In our context, this consists in recovering the coefficients of the parametric resolution, which are rational functions, from the knowledge of their power series expansions. This step was not necessary under the stronger hypotheses of [27, 26, 25, 28, 34, 35]; it is now made necessary by the presence of denominators in the coefficients.

To this effect, we propose an algorithm for the rational reconstruction of multivariate rational functions, which reduces to the usual Padé approximants computation in the univariate case. Its complexity is the best known to date for this question.

Summary.

• We define the notion of parametric geometric resolution, as a means to represent the generic solutions of a parametric system. We give bounds for the complexity of a parametric resolution, in a geometric context where previous results did not apply.

• We propose a probabilistic algorithm for computing such a parametric resolution, and work out a precise control of the probabilistic aspects. This algorithm is valid over any effective field, non necessarily perfect, as was often the case [53, 34, 28]. The complexity is polynomial in the size of the output.

• As intermediate results, we give an algorithm for the rational reconstruction of a multivariate rational function. We also extend the Newton operator given in [28] to a larger context, resulting in a simpler presentation of the computations.

• Finally, we demonstrate the use of these results by treating various real-life applications.

Related work. We have already mentioned that this article is in the continuity of the work of the TERA group [27, 26, 25, 28, 35] and notably of [34]. Let us mention other possible approaches.

• Zero-dimensional solving over a rational function field. The resolution of the system as a zero-dimensional problem over the rational function field k(P_{1}, …, P_{m}) leads to the same output as our algorithm.

For the resolution of zero-dimensional systems, we mention in particular the algorithm of geometric resolution [27, 26, 25, 28, 35]. Other approaches include the computation of Gröbner bases [9, 21], possibly followed by a Rational Univariate Representation [50]. We also mention the linear algebra methods, using the matrices introduced by Macaulay [44] or generalizations thereof [19, 47, 18].

The complexity of these zero-dimensional solving methods is not always known in terms of operations in the base field, which is here the rational function field k(P_{1}, …, P_{m}). Moreover, there is no obvious bound on the degrees of the rational functions that may appear through the computations. Thus it is quite difficult to estimate the complexity of this kind of approach in terms of operations in k, but practice reveals that they are quite costly.

• Exhaustive descriptions. A radically different approach to the question of parametric systems is to give an exhaustive description of the solutions, describing all possibilities of degeneracy. We mention in particular the techniques of dynamic evaluation [16, 29, 14], the comprehensive Gröbner bases [65] and the computation of parametric Gröbner bases proposed in [32] and [45].

Whereas the complexities of the dynamic evaluation method or of Montes’ algorithm are not known to us, the approaches of Grigoriev and Vorobjov and of Weispfenning are known to lead to algorithms of complexity of order d^{O(mn^{2})}, d being a bound on the degree of the input polynomials. A very crude estimation of our complexity result will turn out to be the better bound d^{O(mn)}. Still, the reader must keep in mind that the outputs of these algorithms differ from ours.

Acknowledgements. I wish to thank M. Giusti, who advised my PhD. Thesis, from which this work is taken. It is also a pleasure to thank J. Heintz, G. Lecerf, L.-M. Pardo and B. Salvy for many helpful comments.

2 Notations, main results

We now present our results in a more precise fashion. To this effect, we recall and introduce some notations used in the sequel.

The input system. The base field is denoted by k; we will denote by \overline{k} its algebraic closure. We consider a polynomial system f = (f_{1}, …, f_{n}) in k[P_{1}, …, P_{m}, X_{1}, …, X_{n}], where the indeterminates P = (P_{1}, …, P_{m}) are thought as parameters, or free variables, and X = (X_{1}, …, X_{n}) are thought as algebraic variables.

Given a value p = (p_{1}, …, p_{m}) in \overline{k}^{m}, we will denote by f(p, .) the system f where the indeterminates P are specialized at p.

For the sake of concision, the field \overline{k}(P_{1}, …, P_{m}) will be called K; nevertheless, we will remember that it is a field of rational functions, when geometric arguments or complexity statements are required. For similar concision imperatives, sums such as ∑ u_{i}x_{i} will always be taken for i in 1, …, n.

Geometric objects. Let I be the ideal generated by the polynomials (f_{1}, …, f_{n}) in \overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}]. We wish to exclude the locus where the Jacobian determinant jac(f, X) vanishes. To this effect, we consider J = (I : jac(f, X)^{∞}) the intersection of the primary components of I which do not contain a power of jac(f, X). The ideal J can be defined by polynomials with coefficients in k; we suppose that it is not the trivial ideal (1). The corresponding variety V = V(J) ⊂ 𝔸^{m+n}(\overline{k}) is our object of interest.

Let us denote by π : 𝔸^{m+n}(\overline{k}) → 𝔸^{m}(\overline{k}) the canonical projection on the parameter space 𝔸^{m}(\overline{k}). Then Lazard’s Lemma (see [8] and [46, Proposition 3.2]) shows the following fact:

Fact 1 The restriction of π to each irreducible component of V(J) is a dominant map, not necessarily finite but with generically finite fibers.

Thus, for a generic value p in 𝔸^{m}(\overline{k}), the system f(p, .) = 0 admits finitely many solutions in V, so this situation is actually zero-dimensional over the field K = \overline{k}(P_{1}, …, P_{m}).

We denote by J_{K} the ideal generated in K[X_{1}, …, X_{n}] by the polynomials in J, and by x_{1}, …, x_{n} the images of X_{1}, …, X_{n} modulo J_{K}. Lazard’s Lemma implies that J and J_{K} are radical ideals, and the above discussion states that J_{K} has dimension zero. We now present our basic way of representing such zero-dimensional objects.

Geometric resolutions. The notion of geometric resolution is defined for zero-dimensional systems in [27, 26, 25, 28, 35]. The definition in a general setting is as follows.

Let 𝔎 be any field, 𝔍 a zero-dimensional ideal of 𝔎[X_{1}, …, X_{n}] and x_{i} the image of X_{i} modulo 𝔍, for i in 1, …, n. Then a geometric resolution of the extension 𝔎 → 𝔎[X_{1}, …, X_{n}]/𝔍, if it exists, consists in:

• a primitive element u = ∑ u_{i}x_{i} of 𝔎 → 𝔎[X_{1}, …, X_{n}]/𝔍,

• its monic minimal polynomial Q_{u} ∈ 𝔎[U],

• a parametrization of the algebraic variables in terms of the primitive element. Following [4, 50, 28], we use in priority a parametrization of the form Q′_{u}(u)x_{i} = V_{i}(u) in 𝔎[X_{1}, …, X_{n}]/𝔍, where V_{i} is in 𝔎[U], for i = 1, …, n.

In our particular context, we call parametric geometric resolution, or parametric resolution, a geometric resolution of the extension K → K[X_{1}, …, X_{n}]/J_{K}. Thus, it has coefficients in the rational function field K.

The main results in this paper are the proof of the existence of a primitive element of K → K[X_{1}, …, X_{n}]/J_{K}, bounds on the degrees of the expressions that may appear in the corresponding parametric resolution, a study of the locus where its specialization fails, and a probabilistic algorithm to compute such a resolution.

Complexity notations. We will estimate the complexity of the input and the intrinsic complexity of the geometric objects V and π using the following notations.

• The polynomials (f_{1}, …, f_{n}) are of degree bounded by d, and given by a Straight-Line Program of size L (see [11] for a definition).

• We call deg_{π} the generic cardinality of the fibers of the restriction of π to V, that is the generic number of simple solutions of the specialized systems f(p, .) = 0.

We call deg_{V} the degree of the variety V, using the notion of affine degree given by Heintz in [33].

Using the Bézout inequality of [33], both deg_{π} and deg_{V} can be bounded by the Bézout number d^{n}.

The complexities of our algorithms will be measured using the following notations.

• The notation f ∈ O_{log}(g) means that there exists a constant a such that f is in O(g log(g)^{a}).

• M_{u}(D) denotes the cost of the multiplication of univariate polynomials of degree D, in terms of operations in the base ring. M_{u}(D) can be taken in O(D log D log log D), or O_{log}(D), using the algorithms of Schönhage and Strassen [56] and Schönhage [55].

• M_{s}(D, m) denotes the cost of m-variate series multiplication at precision D. This can be taken less than M_{u}((2D + 1)^{m}) using Kronecker’s substitution, see [39] and [64, ex. 16.16].

If the base field k has characteristic zero, this complexity is in O_{log}(M_{u}(\binom{D+m}{m})), i.e. linear in the size of the series, up to logarithmic factors; see [42].

We make the assumption that there exists a universal constant c < 1 such that M_{s}(D, m) ≤ cM_{s}(2D, m) holds for all D and m.

With these notations, our first result concerns the existence and the complexity of a parametric resolution:

Theorem 1 Let u = ∑ u_{i}x_{i} be a primitive element of K → K[X_{1}, …, X_{n}]/J_{K} with coefficients in \overline{k}, so that the following relations are satisfied in K[X_{1}, …, X_{n}]/J_{K}:

Q_{u}(u) = 0, { Q′_{u}(u)x_{1} = V_{1}(u), ; ⋮ ; Q′_{u}(u)x_{n} = V_{n}(u),

where

• the polynomials Q_{u} and V_{1}, …, V_{n} are in K[U],

• the polynomial Q_{u} is the monic minimal polynomial of u and the polynomials V_{i} have degree less than Q_{u}.

Then the polynomial Q_{u} has degree deg_{π}. We recall that K is \overline{k}(P_{1}, …, P_{m}). Then all numerators and the least common multiple of the denominators of all the coefficients of the polynomials Q_{u}, V_{1}, …, V_{n} have degree in P_{1}, …, P_{m} at most deg_{V} ≤ d^{n}. Furthermore, if the cardinality of k is greater than d^{n}(2d^{n} + nd + 1), then there exists a primitive element u = ∑ u_{i}x_{i} of K → K[X_{1}, …, X_{n}]/J_{K} with coefficients in k. In this case, the polynomials Q_{u}, V_{1}, …, V_{n} have their coefficients in the subfield k(P_{1}, …, P_{m}) of K.

The second theorem is of algorithmic nature. Its complexity statement is notably given in terms of a degree denoted by deg_{u}. This degree is defined as the maximum of the degrees of the coefficients that appear in the parametric resolution for a primitive element u. Using Theorem 1, deg_{u} is bounded by deg_{V}, and thus by d^{n}.

Theorem 2 Assume that the cardinality of k is greater than d^{n}(2d^{n} + nd + 1). There exists a probabilistic algorithm which computes a parametric resolution for a primitive element ∑ u_{i}x_{i} with coefficients in k, through the following steps:

• the first step consists in computing a geometric resolution of the simple zeros of the specialized system f(p, .) = 0, where p = (p_{1}, …, p_{m}) is a point in k^{m}.

• the second step is a formal Newton lifting process, which requires

O_{log}((nL + n^{4})M_{u}(deg_{π})M_{s}(4 deg_{u}, m) + nm^{2} deg_{π} M_{u}(deg_{u})M_{s}(4 deg_{u}, m − 1))

operations in k, where deg_{u} is the maximum of the degrees in P_{1}, …, P_{m} of the numerators and denominators of the coefficients of the polynomials Q_{u}, V_{1}, …, V_{n}.

The algorithm chooses 3m − 1 values in k, including p = (p_{1}, …, p_{m}); if Γ is a subset of k, and these values are chosen in Γ^{3m−1}, then the algorithm succeeds for all choices except at most 110nd^{4n}|Γ|^{3m−2}, for d ≥ 2 and n ≥ 2.

We do not give more details on the first step of this algorithm, the resolution of a zerodimensional system over k, as efficient solutions are already known (see Section 4.1). We concentrate on the lifting step, which is the most expensive. Since deg_{u} is bounded by d^{n}, the dependence of our complexity in d is of order d^{O(nm)} base field operations. Theorem 2 shows that the complexity of the lifting step is polynomial in the size of the output. More precisely, using the algorithms for fast univariate polynomial and power series arithmetic mentioned above, we obtain the following corollary. The proof is straightforward, using for instance [42, Lemma 3].

Corollary 1 If k has characteristic zero, then, in terms of operation in k, the complexity of the lifting step is in

O_{log} ((nL + n^{4} + nm^{2}) deg_{π} \binom{4 deg_{u} +m}{m}) .

The size of the output is within O(n deg_{π} \binom{deg_{u} +m}{m}) elements in k. Thus, the complexity of the lifting step is at most quartic in the size of the output.

Organization of the paper.

• In Section 3, we prove Theorem 1, and give estimates the degrees of a degeneracy locus, which will help quantify the probability of success of our algorithm.

• Section 4 presents the most important algorithmic tools we use, a new version of a formal Newton operator and a new algorithm for the reconstruction of a multivariate rational function.

• Section 5 presents the main algorithm, with the proof of its complexity and probability of success. This will prove Theorem 2.

• The algorithm is implemented in Magma [2], on the basis of the Kronecker package [28, 40]; the applications we treated and the practical behavior of the implementation are presented in Section 6.

3 Degree estimates for the parametric resolution

In this section, we establish the existence of a parametric resolution and prove the bounds on the degree in P_{1}, …, P_{m} of its coefficients. We also give a bound on the degree of a hypersurface in the parameter space which contains the points where a parametric resolution cannot be specialized; this result controls the probability of success of the algorithm given in Section 5.

The organization is as follows. Subsection 3.1 establishes some technical algebraic results, that are used in the sequel. In Subsection 3.2, we prove bounds on the complexity of the minimal polynomial of a function of X_{1}, …, X_{n}. We use this result in Subsection 3.3 to give the proof of Theorem 1. Finally, Subsection 3.4 is devoted to the study of the aforementioned degeneracy hypersurface.

Notations. Throughout this section, we use some additional notations. Let us also write again the definitions of the most important objects used up to now.

• The ideal I ⊂ \overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}] is generated by the polynomials f; J is its saturation with respect to the jacobian determinant of f. The variety V is the zero-set of J.

• We recall that K denotes the field \overline{k}(P_{1}, …, P_{m}), and J_{K} the extension of the ideal J in K[X_{1}, …, X_{n}]. For the sake of shortness, B denotes the finite-dimensional quotient algebra K[X_{1}, …, X_{n}]/J_{K}.

• We frequently use the notion of minimal and characteristic polynomial of elements of B, which we now recall. Let then g be in B, and U a new variable.

The minimal polynomial of g is the monic generator of the ideal {P ∈ K[U], P (g) = 0 in B}. This polynomial is denoted by Q_{g}; we will write Q_{g} = M_{g}/D_{g}, where M_{g} is primitive in \overline{k}[P_{1}, …, P_{m}][U], and D_{g} ∈ \overline{k}[P_{1}, …, P_{m}] is its leading coefficient.

The characteristic polynomial of g is the characteristic polynomial of the endomorphism of multiplication by g in B. It denoted by χ_{g}. We will take χ_{g} = Ξ_{g}/Θ_{g}, where Ξ_{g} is primitive in \overline{k}[P_{1}, …, P_{m}][U] and Θ_{g} ∈ \overline{k}[P_{1}, …, P_{m}] is its leading coefficient.

• If Q is a polynomial or rational function depending on indeterminates P_{1}, …, P_{m} and p = (p_{1}, …, p_{m}) is a point in \overline{k}^{m}, we call specialization of Q at p the polynomial Q with all coefficients specialized at p, if possible. It is denoted by Q(p), or Q(., p).

We extend naturally this denomination to the specialization of a parametric resolution at a point p in \overline{k}^{m}. We obtain a family of polynomials in \overline{k}[U].

• The notation u will denote a n-uple (u_{1}, …, u_{n}) in \overline{k}^{n}.

In Subsection 3.4, we introduce some new variables λ_{1}, …, λ_{n}. Then, just as above, a notation such as Ξ(., u) denotes a polynomial, here Ξ, where the variables (λ_{1}, …, λ_{n}) are specialized on the values (u_{1}, …, u_{n}).

For the sake of simplicity, all complexity results are stated in terms of the Bézout number d^{n}. More precise results could be obtained using the product of the degrees of the polynomials f_{1}, …, f_{n}.

3.1 Preliminaries

Since we make no assumption on the field k, some extra care is required concerning the separability of the extensions we consider; this is the object of the following lemma. As a consequence, we deduce the existence of a primitive element for the extension K → B.

Lemma 1 The extension K → B is a product of separable field extensions, and has degree deg_{π}.

Proof. Recall that the ideal J is the defining ideal of the variety V. Let us write the primary decomposition of J in \overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}], J = ∩_{ℓ}J_{ℓ}, where the sum is taken for indices ℓ in some set L. We now use this decomposition to reduce the proof to the prime case.

For each ℓ, let J_{K,ℓ} be the extension of J_{ℓ} in K[X_{1}, …, X_{n}] = \overline{k}(P_{1}, …, P_{m})[X_{1}, …, X_{n}]. Since J is radical, all ideals J_{ℓ} are prime. Fact 1 states that the restriction of the projection π to each V(J_{ℓ}) is dominant, so the ideals J_{ℓ} contain no element of \overline{k}[P_{1}, …, P_{m}]. The extended ideals J_{K,ℓ} then remain prime in \overline{k}(P_{1}, …, P_{m})[X_{1}, …, X_{n}]. From this, we easily deduce the isomorphism

\overline{k}(P_{1}, …, P_{m})[X_{1}, …, X_{n}]/J_{K,ℓ} ≃ 𝔣𝔯 (\overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}]/J_{ℓ}) ,

where 𝔣𝔯(.) denotes the fraction field of an integral ring.

Using Fact 1, we also deduce that the extended ideal J_{K} is the intersection of the extended ideals J_{K,ℓ} for ℓ in L, and that any two distinct extended ideals J_{K,ℓ} and J_{K,ℓ′} are coprime. These results yield the following sequence of isomorphisms:

B = \overline{k}(P_{1}, …, P_{m})[X_{1}, …, X_{n}]/J_{K}
≃ ∏_{ℓ∈L} \overline{k}(P_{1}, …, P_{m})[X_{1}, …, X_{n}]/J_{K,ℓ}
≃ ∏_{ℓ∈L} 𝔣𝔯 (\overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}]/J_{ℓ}) .

By construction, the Jacobian determinant jac(f, X) is invertible on a dense subset of each zero-set V(J_{ℓ}). Then the Jacobian criterion given in [17, Corollary 16.16] states that each of the field extensions

\overline{k}(P_{1}, …, P_{m}) → 𝔣𝔯 (\overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}]/J_{ℓ})

is separable (see [46] for a similar statement). Thus the first assertion of the lemma is proved; we now show that the degree of the extension K → B is indeed deg_{π}.

Using the separability condition obtained above, Proposition 1 in [33] shows that for each ℓ the degree of the extension

\overline{k}(P_{1}, …, P_{m}) → 𝔣𝔯 (\overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}]/J_{ℓ})

is the generic cardinality of the fibers of π restricted to V(J_{ℓ}). Thus the sum of their degrees is deg_{π}. This proves the second point of the lemma. □

Corollary 2

• An element in B is primitive for K → B if and only if its characteristic polynomial has no multiple root.

• There exists a primitive element of the extension K → B of the form ∑ u_{i}x_{i}, with coefficients in \overline{k}.

Proof. An element is primitive for K → B if and only if its minimal polynomial equals its characteristic polynomial. Separability implies that the minimal polynomial of an element in B has no multiple root, which proves the first assertion.

The second result is folklore, and follows from both facts that K → B is a product of separable field extensions, and that \overline{k} is an infinite subfield of K. See for instance the proof of [13, Theorem 2.1.5], which can be transcripted verbatim to the present situation. □

3.2 Degree of an eliminating polynomial

We now address our first complexity question: we consider the complexity of the minimal polynomial in B of a element g ∈ \overline{k}[X_{1}, …, X_{n}]. The following proposition is an extension of [53, Proposition 1]: in our situation, the presence of denominators in the minimal polynomial of the element g deserves special attention.

Proposition 1 Let g be a polynomial in \overline{k}[X_{1}, …, X_{n}]. Then the minimal polynomial Q_{g} ∈ K[U] of the image of g in B can be written M_{g}/D_{g}, where • M_{g} is primitive in \overline{k}[P_{1}, …, P_{m}][U], and D_{g} ∈ \overline{k}[P_{1}, …, P_{m}] is its leading coefficient;

• seen in \overline{k}[P_{1}, …, P_{m}, U], M_{g} has total degree at most deg_{V} deg g.

If g belongs to k[X_{1}, …, X_{n}], then M_{g} and D_{g} may be taken with coefficients in k.

Proof. Let φ be the morphism

\begin{matrix}V & → & 𝔸^{m+1}(\overline{k}) \\ (p, x) & ↦ & (p, g(x)). \end{matrix}

The closure W of its image is an hypersurface in 𝔸^{m+1}(\overline{k}) of degree at most deg_{V} deg g. We will prove that an equation defining this hypersurface yields the minimal polynomial of g in B.

Let thus M_{g} in \overline{k}[P_{1}, …, P_{m}, U] be a squarefree polynomial of degree at most deg_{V} deg g defining the hypersurface W. We will consider this polynomial in either \overline{k}[P_{1}, …, P_{m}, U] or \overline{k}[P_{1}, …, P_{m}][U]; we let D_{g} ∈ \overline{k}[P_{1}, …, P_{m}][U] be the leading coefficient of M_{g}, when M_{g} is considered univariate in U.

Let Q_{g} ∈ K[U] be the monic minimal polynomial of the multiplication by g in B. We first notice that the polynomial Q_{g} can be written M/D, where M and D respectively belong to \overline{k}[P_{1}, …, P_{m}][U] and \overline{k}[P_{1}, …, P_{m}], and M is primitive. Then, we prove that M = M_{g}.

• M_{g}(g) is identically zero on V = V(J). Since J is radical, M_{g}(g) belongs to J, so M_{g}(g) is zero in B. This implies that M_{g} is a multiple of Q_{g} = M/D in K[U], that is, there exists an equality aM = bM_{g}, where a belongs to \overline{k}[P_{1}, …, P_{m}, U] and b belongs to \overline{k}[P_{1}, …, P_{m}]. Since M is primitive in \overline{k}[P_{1}, …, P_{m}][U], it divides M_{g} in \overline{k}[P_{1}, …, P_{m}, U].

• Conversely, M(g) belongs to the ideal J_{K} ∩ \overline{k}[P_{1}, …, P_{m}, X_{1}, …, X_{n}]. Since no prime component of the ideal J contains an element in \overline{k}[P_{1}, …, P_{m}], this intersection is exactly J. Consequently, M(g) belongs to J, so M vanishes on W, that is M_{g} divides M in \overline{k}[P_{1}, …, P_{m}, U].

This implies that M = M_{g} up to a factor in \overline{k}, so that Q_{g} = M/D = M_{g}/D_{g}.

The algebra B is defined by polynomials with coefficients in k, so if g has its coefficients in k then its minimal polynomial Q_{g} belongs to k(P_{1}, …, P_{m})[U], i.e. M and D can be taken with coefficients in k. This proves the proposition. □

Remark 1. The proof shows the following fact: for every p = (p_{1}, …, p_{m}) ∈ \overline{k}^{m} which does not cancel the denominator D_{g}, the polynomial Q_{g}(p, U) ∈ \overline{k}[U] vanishes on the values taken by g in the fiber π^{−1}(p) ∩ V. We use this remark in Subsection 3.4.

3.3 Degree of a parametric resolution

Using well-known techniques [39, 44, 48, 4], we can recover a parametric resolution through the computation of the minimal polynomial of a “generic primitive element” of K → B; this minimal polynomial is also called a u-resultant [12] or a Chow form [4]. As a consequence, the degree bound obtained in the previous proposition will apply for the whole resolution.

These results are summarized in the following proposition, which gives the first part of Theorem 1.

Proposition 2 Let u = ∑ u_{i}x_{i} be a primitive element of K → B with coefficients in \overline{k}, so that the following relations are satisfied in B:

Q_{u}(u) = 0, { Q′_{u}(u)x_{1} = V_{1}(u), ; ⋮ ; Q′_{u}(u)x_{n} = V_{n}(u),

where

• the polynomials Q_{u} and V_{1}, …, V_{n} are in K[U],

• the polynomial Q_{u} is the monic minimal polynomial of u, the polynomials V_{i} have degree less than Q_{u}.

Then Q_{u} has degree deg_{π}.

We recall that K is the field \overline{k}(P_{1}, …, P_{m}). Then all numerators and the least common multiple of the denominators of all the coefficients of the polynomials Q_{u}, V_{1}, …, V_{n} have degree at most deg_{V}.

Furthermore, if u has coefficients in k, then the polynomials Q_{u}, V_{1}, …, V_{n} have their coefficients in the subfield k(P_{1}, …, P_{m}) of K.

Before proving the proposition, we mention the following useful consequences. Given a primitive element u = ∑ u_{i}x_{i}, the separability of its minimal polynomial Q_{u} implies that Q′_{u} is invertible modulo Q_{u}, so the parametrization introduced in the above proposition makes sense. Inverting Q′_{u} modulo Q_{u}, we can write the alternative parametrization, reminiscent of the Shape Lemma [24]:

Q_{u}(u) = 0, { x_{1} = W_{1}(u), ; ⋮ ; x_{n} = W_{n}(u).

Both forms of parametrizations will be used in the sequel, so we give them specific names. Formulae similar to those given in Proposition 2 can be found in Kronecker’s work [39], hence the following denomination.

Definition 1 Let u = ∑ u_{i}x_{i} be a primitive element of K → B with coefficients in \overline{k}.

• We call Kronecker parametrization the vector R_{u} = [Q_{u}, V_{1}, …, V_{n}] defined in Proposition 2.

• We call Shape Lemma parametrization the vector S_{u} = [Q_{u}, W_{1}, …, W_{n}] defined above.

Proposition 2 could be used to give a bound on the degrees of the coefficients in a Shape Lemma parametrization, at best quadratic in the Bézout number, that is, much worse than the bound for the Kronecker form. Practice reflects this point: introducing the normalization Q′_{u} has the effect to lower the size of the parametrization obtained through the Shape Lemma, as was already noticed in [4, 50].

Remark 2. If u = ∑ u_{i}x_{i} is a primitive element of K → B, and if we denote U = ∑ u_{i}X_{i}, then the equality between radical ideals

J_{K} = (Q_{u}(U), X_{1} − W_{1}(U), …, X_{n} − W_{n}(U))

holds in K[X_{1}, …, X_{n}]. This remark is used in the next subsection.

Proof of Proposition 2: introduction of a generic linear form. We deduce Proposition 2 from the study of the minimal (or characteristic) polynomial of a linear combination of the variables X with generic coefficients. To this effect, we extend the base field, and correspondingly extend the notations. The objects introduced below will be considered again in the next subsection.

Let (λ_{1}, …, λ_{n}) be new indeterminates, which will be used as coefficients of the generic linear form. We denote by k_{Λ} the field k(λ_{1}, …, λ_{n}) and proceed to extend all previous constructions to this new base field; they will be denoted by a subscript Λ. Since no confusion can occur, we still use the letters P, X for indeterminates.

Thus, we denote by I_{Λ} the ideal generated in \overline{k_{Λ}}[P_{1}, …, P_{m}, X_{1}, …, X_{n}] by the polynomials f = f_{1}, …, f_{n}, and by J_{Λ} its saturation with respect to the Jacobian determinant jac(f, X). The zero-set of J_{Λ} is denoted by V_{Λ}. It is a routine check that J_{Λ} is the extension of J in \overline{k_{Λ}}[P_{1}, …, P_{m}, X_{1}, …, X_{n}]; using the definition of degree from [33], we deduce that the degree of V_{Λ} equals the degree of V.

We denote by K_{Λ} the field \overline{k_{Λ}}(P_{1}, …, P_{m}), which is the analogous to the field K used up to now. Similarly, J_{K,Λ} denotes the ideal generated in K_{Λ}[X_{1}, …, X_{n}] by the polynomials in J_{Λ}, B_{Λ} is the quotient K_{Λ}[X_{1}, …, X_{n}]/J_{K,Λ}, and x_{i} is the image of X_{i} in B_{Λ}.

Up to now, the variables λ_{1}, …, λ_{n} have played no active role. We now denote by U_{Λ} the generic linear form ∑ λ_{i}X_{i}, and by u_{Λ} = ∑ λ_{i}x_{i} its image in B_{Λ}. Given a new indeterminate U_{Λ}, we denote by χ_{Λ} ∈ K_{Λ}[U_{Λ}] the characteristic polynomial of u_{Λ} in B_{Λ}.

The polynomial χ_{Λ} is the eliminating polynomial for a generic linear form we wanted to introduce. We now present its basic properties, and see how to use it for obtaining a parametric resolution.

Lemma 2 For any u = (u_{1}, …, u_{n}) in \overline{k}^{n}, the characteristic polynomial of ∑ u_{i}x_{i} in B is the specialization χ_{Λ}(., u).

Proof. Any K-basis of B is also a K_{Λ}-basis of B_{Λ}. Consequently, the polynomial χ_{Λ} is the determinant of U_{Λ}I − ∑ λ_{i}M_{x_{i}}, where M_{x_{i}} is the matrix of multiplication by x_{i} in such a basis, and I is the identity matrix. This proves the lemma. □

Lemma 3 χ_{Λ} is an homogeneous polynomial of degree deg_{π} in (U_{Λ}, λ_{1}, …, λ_{n}), monic in U_{Λ}. It can be written χ_{Λ} = Ξ_{Λ}/Θ_{Λ}, where Ξ_{Λ} belongs to k[P_{1}, …, P_{m}, λ_{1}, …, λ_{n}][U_{Λ}], Θ_{Λ} belongs to k[P_{1}, …, P_{m}], both polynomials have degree in (P_{1}, …, P_{m}) bounded by deg_{V}, and Θ_{Λ} is the leading coefficient of Ξ_{Λ}.

Proof. The first point is an easy consequence of the proof of the previous lemma, so we concentrate on the second assertion.

By Corollary 2, there exists a primitive element u = ∑ u_{i}x_{i} of K → B with coefficients in \overline{k}. Then, also by Corollary 2, its characteristic polynomial χ_{u} has no multiple root. Using the specialization property of the previous lemma, this shows that the polynomial χ_{Λ} cannot have multiple roots. Thus it coincides with the minimal polynomial of u_{Λ}. We can then apply Proposition 1 to the variety V_{Λ} defined over the field k_{Λ} to conclude the proof of the lemma. □

Concluding the proof. The polynomial χ_{Λ}(U_{Λ}) belongs to the ideal J_{K,Λ}, so it can be written ∑ g_{j}F_{j}, where (F_{j}) are generators of J in k[P_{1}, …, P_{m}, X_{1}, …, X_{n}] and (g_{j}) are polynomials in K_{Λ}[X_{1}, …, X_{n}]. Since χ_{Λ} is a polynomial in (λ_{1}, …, λ_{n}), the polynomials g_{j} can be taken polynomial in (λ_{1}, …, λ_{n}) too. The derivative of χ_{Λ}(U_{Λ}) with respect to λ_{i} is ∑ \frac{∂g_{j}}{∂λ_{i}} F_{j}; it can also be written (\frac{∂χ_{Λ}}{∂U_{Λ}} X_{i} + \frac{∂χ_{Λ}}{∂λ_{i}})(U_{Λ}).

Since u is primitive, its minimal polynomial Q_{u} coincides with its characteristic polynomial χ_{u}. Lemma 2 then shows that the specialization (λ_{1}, …, λ_{n}) ← (u_{1}, …, u_{n}) in χ_{Λ} and \frac{∂χ_{Λ}}{∂U_{Λ}} are Q_{u} and Q′_{u}. We take for V_{i} the specialization of −\frac{∂χ_{Λ}}{∂λ_{i}} and let U be ∑ u_{i}X_{i}. Since the polynomials \frac{∂g_{j}}{∂λ_{i}} are polynomial in the variables λ_{1}, …, λ_{n}, the specialization (λ_{1}, …, λ_{n}) ← (u_{1}, …, u_{n}) shows that (Q′_{u}X_{i} − V_{i})(U) belongs to J_{K}, so vanishes in B.

We have thus obtained a parametrization. Let us prove that is has the announced degree properties. This will conclude the complexity analysis: the conditions given in Proposition 2 obviously impose uniqueness for V_{1}, …, V_{n}.

First, note that Lemma 3 gives the bound on the degree in U_{Λ} of χ_{Λ} and its derivatives.

The polynomial χ_{Λ} is polynomial in (λ_{1}, …, λ_{n}), so differentiation with respect to λ_{i} does not alter the degrees of the rational functions in (P_{1}, …, P_{m}) that appear. Thus, all numerators in Q_{u}, V_{1}, …, V_{n} have degree at most deg_{V} by Lemma 3. Moreover, all denominators of the coefficients of the polynomials Q_{u}, V_{1}, …, V_{n} divide the denominator of χ_{Λ}, which is Θ_{Λ}, so they have degree at most deg_{V} by Lemma 3 again. Thus the complexity analysis is complete.

We conclude the proof by a trivial remark: if the coefficients (u_{1}, …, u_{n}) are in the base field k, our construction shows that Q_{u}, V_{1}, …, V_{n} have coefficients in k too. □

3.4 Lucky specializations

A parametric resolution can be specialized on an open subset of the parameter space 𝔸^{m}(\overline{k}), to give a description of the simple solutions of the corresponding specialized system. In this section, we define a universal discriminant locus describing the values where the specialization fails, and give a bound on its degree. This bound will be used to control some error probabilities, so it is given in terms of the input data (n, d).

It will be useful to consider both forms of parametrization introduced in the previous subsection, the Kronecker form R_{u} and the Shape Lemma form S_{u}.

Proposition 3 There exists a non-zero polynomial Δ in k[P_{1}, …, P_{m}, λ_{1}, …, λ_{n}] of degree at most d^{n}(2d^{n} + nd + 1) in (P_{1}, …, P_{m}) and 2d^{2n} in (λ_{1}, …, λ_{n}) such that:

• for all p in \overline{k}^{m} and u = (u_{1}, …, u_{n}) in \overline{k}^{n}, if Δ(p, u) is not zero, then u = ∑ u_{i}x_{i} is a primitive element of K → B and p cancels none of the denominators in either R_{u} or S_{u}. Furthermore, the system f(p, .) = 0 has deg_{π} simple solutions, which are described by the specialization of either R_{u} or S_{u} at p.

• for all p in \overline{k}^{m}, if Δ(p, .) is not zero, and U = ∑ u_{i}X_{i} is a linear form that induces a primitive element for the extension corresponding to the simple solutions of the specialized system f(p, .) = 0, then, taking u = (u_{1}, …, u_{n}), Δ(p, u) is not zero, so the same conclusion holds.

Before proving the proposition, we consider its following consequence. If the cardinality of k is greater than 2d^{n}(2d^{n} + nd + 1), then there exists a value (p, u) in k^{m+n} which does not cancel Δ. Then the corresponding element u = ∑ u_{i}x_{i} is a primitive element of K → B. This concludes the proof of Theorem 1.

Proof of Proposition 3. We consider again the variables λ_{1}, …, λ_{n}, and the polynomials χ_{Λ}, Ξ_{Λ} and Θ_{Λ} introduced in the end of the previous subsection, notably Lemma 3.

Let R_{Λ} be the resultant of Ξ_{Λ} and Ξ′_{Λ}, which is a polynomial in k[P_{1}, …, P_{m}, λ_{1}, …, λ_{n}]. From the proof of Lemma 3, we see that this polynomial is non-zero. The bounds given in Lemma 3 show that R_{Λ} has degree at most 2 deg_{π} deg_{V} ≤ 2d^{2n} in (P_{1}, …, P_{m}) and 2 deg_{π}^{2} ≤ 2d^{2n} in (λ_{1}, …, λ_{n}). We now show that this polynomial controls the denominators appearing in a parametric resolution, in either Kronecker or Shape Lemma form.

We recall that the characteristic polynomial of an element u in B is written χ_{u} = Ξ_{u}/Θ_{u}, where Ξ_{u} is primitive in \overline{k}[P_{1}, …, P_{m}][U] and Θ_{u} ∈ \overline{k}[P_{1}, …, P_{m}] is its leading coefficient.

Lemma 4 For any u = (u_{1}, …, u_{n}) in \overline{k}^{n}, the resultant of Ξ_{u} and Ξ′_{u} divides R_{Λ}(., u). If u = ∑ u_{i}x_{i} is a primitive element of K → B, the numerator of the discriminant of Q_{u} divides R_{Λ}(., u).

Proof. Lemma 2 shows that the characteristic polynomial χ_{u} = Ξ_{u}/Θ_{u} is equal to the specialization Ξ_{Λ}(., u)/Θ_{Λ}. Since Ξ_{u} is primitive, it divides Ξ_{Λ}(., u), the possible factor lying in \overline{k}[P_{1}, …, P_{m}]. This implies that the resultant of Ξ_{u} and Ξ′_{u} divides the resultant of Ξ_{Λ}(., u) and Ξ_{Λ}(., u)′. Since by construction the leading term of Ξ_{Λ} does not depend on (λ_{1}, …, λ_{n}), this resultant is the specialization R_{Λ}(., u). This proves the first point.

If u is a primitive element of K → B, its minimal polynomial Q_{u} coincides with its characteristic polynomial Ξ_{u}/Θ_{u}. The numerator of the discriminant of Q_{u} then divides the resultant of Ξ_{u} and Ξ′_{u}, which proves the second point. □

Lemma 5 Let u = (u_{1}, …, u_{n}) be in \overline{k}^{n}. If u = ∑ u_{i}x_{i} is a primitive element of K → B, then all the denominators in R_{u} and S_{u} divide Θ_{Λ}R_{Λ}(., u).

Proof. The proof of Proposition 2 shows that all denominators appearing in R_{u} divide Θ_{Λ}, which proves the assertion for R_{u}. Going from the representation R_{u} to the representation S_{u} requires to invert Q′_{u} modulo Q_{u}. The denominators that appear in the inversion divide the numerator of the discriminant of Q_{u} and Q′_{u}, so they divide R_{Λ}(., u) by the lemma above. This proves the result. □

Concluding the proof. It remains to exclude the possible degenerate points of V. The Jacobian determinant jac(f, X) has degree at most nd. The intersection of its zero-set with V is a variety of degree at most nd deg_{V} ≤ nd^{n+1}, whose image by π is contained in an hypersurface of 𝔸^{m}(\overline{k}). Let S by a polynomial in k[P_{1}, …, P_{m}] of degree at most nd^{n+1} defining such an hypersurface.

We define Δ as the product SΘ_{Λ}R_{Λ}, which has the requested degree in (P_{1}, …, P_{m}) and (λ_{1}, …, λ_{n}). We now prove that Δ fulfills our requirements.

• Let (p, u) be a (m + n)-uple in \overline{k}^{m+n} which does not cancel Δ, and u the linear form ∑ u_{i}x_{i}.

We first briefly describe the fiber π^{−1}(p) ∩ V. Since p does not cancel Δ, the Jacobian determinant of f vanishes on none of the points in the fiber π^{−1}(p) ∩ V ⊂ V(f). Then the Jacobian criterion shows that these points are in finite number. By [33, Proposition 1], this number is at most deg_{π}.

We now turn to the specialization of the parametric resolution. Lemma 4 shows that the resultant of Ξ_{u} and Ξ′_{u}, and thus the discriminant of the characteristic polynomial χ_{u}, cannot be zero, so by Corollary 2, u is a primitive element of K → B. Lemma 5 then shows that p cancels none of the denominators in the corresponding resolutions R_{u} and S_{u}.

Since the numerator of the discriminant of the polynomial Q_{u}(p) ∈ \overline{k}[U] is not zero, Q_{u}(p) has deg_{π} distinct roots. Consequently, the specialization of either R_{u} or S_{u} at p describes deg_{π} distinct points.

Let f be a polynomial in J. From Remark 2, there exist polynomials (g_{1}, …, g_{n+1}) in K[X_{1}, …, X_{n}] such that f = ∑ g_{i}(X_{i} − W_{i}(U)) + g_{n+1}Q_{u}(U), where U = ∑ u_{i}X_{i}. The rewriting process introduces no new denominator; since p cancels none of the denominators in S_{u}, it cancels none of the denominators in this equality. Consequently, the deg_{π} points described by the specialization of S_{u} at p cancel the polynomial f(p, .).

Thus, these points are included in the fiber π^{−1}(p) ∩ V. Using the above upper bound on the cardinality of this fiber, we obtain the first part of the proposition.

• Let now p be in \overline{k}^{m}, and assume that Δ(p, .) is not zero. There exists u′ = (u′_{1}, …, u′_{n}) in \overline{k}^{n} which does not cancel this polynomial, so the previous point shows that the specialized system f(p, .) = 0 admits deg_{π} simple solutions.

Take now U = ∑ u_{i}X_{i} any linear form inducing a primitive element for the extension of \overline{k} generated by the simple solutions of this system. Then U takes deg_{π} distinct values on these points.

We denote by Ξ_{Λ}(p, u) ∈ \overline{k}[U] the polynomial Ξ_{Λ} whose coefficients are specialized at (p, u). Since Θ_{Λ}(p) is not zero, Ξ_{Λ}(p, u) has full degree deg_{π} in its leading variable U, so the resultant of Ξ_{Λ}(p, u) with its derivative is the specialization of the generic resultant R_{Λ} at (p, u). Remark 1 shows that Ξ_{u}(p) vanishes on the deg_{π} distinct values taken by U, and Lemma 2 shows that it is also the case for Ξ_{Λ}(p, u). Consequently, the resultant R_{Λ}(p, u) is not zero. This concludes the proof. □

4 Outlook of the algorithm

Our purpose is now to compute a parametric resolution. To this effect, we propose the following algorithm, reminiscent of both numerical root-finding techniques and Hensel lifting methods.

1. Initial estimation: given p = (p_{1}, …, p_{m}) in k^{m}, compute a geometric resolution of the simple roots of the specialized system f(p, .) = 0.

2. Approximation: starting from this specialized solution, approximate the coefficients of the corresponding parametric resolution in a ring of formal power series.

3. Reconstruction: recover the coefficients of the parametric resolution from their power series expansion at (p_{1}, …, p_{m}).

The core of the main algorithm is given in the next section. Here, in the following subsections, we detail our solutions to the three points above: the resolution of the specialized system, a formal Newton operator, and the reconstruction of a rational function from its power series expansion. These parts are largely independent.

It will be useful to recall the following complexity notations:

• M_{u}(D) denotes the cost of the multiplication of univariate polynomials of degree D.

• M_{s}(D, m) denotes the cost of m-variate series multiplication at precision D.

4.1 Computing the initial resolution

Given a point p in k^{m}, the first task is to compute a geometric resolution of the simple solutions of the specialized system f(p, .) = 0. Such a process will be denoted Resolution(f, p). In our main algorithm, its output will be denoted [u, R_{u}^{0}]. In this case:

• u = (u_{1}, …, u_{n}) is such that ∑ u_{i}X_{i} induces a primitive element for the simple solutions of the system f(p, .) = 0

• R_{u}^{0} is a vector of polynomials [q_{u}, v_{1}, …, v_{n}] in k[U], forming a geometric resolution for these points.

The superscript 0 in R_{u}^{0} indicates that this resolution is thought as the truncation of a generic resolution at precision 0 around (p_{1}, …, p_{m}).

Several tools are available to compute this initial resolution. In the spirit of the present paper, we mention in particular the algorithm of geometric resolution, initially in [26], see also [27, 25]. A simplified and improved version is given in [28], together with the description of its Magma implementation, called Kronecker [40].

This algorithm applies in the same model as ours: the input system is given by a Straight-Line Program of size L. Its complexity depends on a geometric quantity attached to the system, called δ, which is at most d^{n}. With this notation, Theorem 1 in [28] states that the resolution can be computed within O((nL + n^{4})M_{u}(dδ)^{2}) operations in k. The algorithm is of probabilistic nature, and a probability analysis is done for a similar algorithm in [35].

Let us mention other approaches to this question, which were already presented in the introduction. Popular methods rely on Gröbner bases computations [21, 20], either for a lexicographic ordering, or followed by the computation of a Rational Univariate Representation [50]. Other approaches include the computation of u-resultants by means of linear algebra methods, based on generalizations of Sylvester or Bézout matrices [47, 49].

4.2 Lifting the resolution

Knowing the initial estimate, the successive approximations of the parametric resolution are obtained through a formal Newton approximation process. It consists in computing a sequence of resolutions, whose coefficients are the successive Taylor series expansions of the coefficients of the requested parametric resolution. This subsection is devoted to present the details of this Newton lifting operator, and the complexity of an elementary lifting step.

We first present the method in a general setting: we consider the lifting modulo the powers of any ideal of a given coefficient ring. In a second time, we apply this result to our specific problem, where the lifting is done modulo the powers of the maximal ideal of a m-variate power series ring.

Just as the numerical Newton operator doubles the number of digits of accuracy at each step, the formal version doubles the precision of the power series at each step. Thus, the series we compute have successive precisions 1, 2, …, 2^{κ}, … Proposition 5 below gives the complexity of the basic step of this operator: using the complexity notations given in the introduction, lifting the Taylor series from precision 2^{κ} to precision 2^{κ+1} requires

O((nL + n^{4})M_{u}(deg_{π})M_{s}(2^{κ+1}, m))

base field operations. In characteristic zero, using fast arithmetic, this complexity is linear, up to logarithmic factors, in the size of the output.

Our use of the formal Newton operator is in the continuity of notably [27, 26, 25, 28, 34, 35]. In particular, a ready-to-implement formal Newton operator was given in the article [28]. Our method generalizes this algorithm to a wider class of representations, triangular sets representations. We do not use the full generality of this method here; this is the subject of [59]. Yet, we stress the fact that the presentation of the computations has now become both more general and simpler.

Generalist presentation. We temporarily broaden our framework: we consider a commutative ring with unity A, an ideal I of A and some polynomials F = (F_{1}, …, F_{N}) in A[X_{1}, …, X_{N}]. We will describe how to approximate some “solutions” of the system F modulo the powers of I.

To this effect, let t = (t_{1}, …, t_{N}) be polynomials in A[X_{1}, …, X_{N}]. We suppose that t forms what will be called a triangular set: for each j, the polynomial t_{j} is monic in X_{j}, reduced with respect to (t_{j+1}, …, t_{N}), and depends only on the variables X_{j}, …, X_{N}.

This triangular set is meant to represent some solutions of the system F without multiplicities, in the sense that:

• (H1) there exists a N × N matrix A with entries in A[X_{1}, …, X_{N}] such that the equality F = At holds;

• (H2) the Jacobian determinant jac(F) is invertible in A[X_{1}, …, X_{N}]/(I, t_{1}, …, t_{N}).

For any positive integer κ, we denote (t_{1}^{κ}, …, t_{N}^{κ}) the images of the polynomials (t_{1}, …, t_{N}) in A/I^{2^{κ}}[X_{1}, …, X_{N}]; these images are the successive approximations of (t_{1}, …, t_{N}) we are interested in.

Let κ be a fixed positive integer. We suppose that (t_{1}^{κ}, …, t_{N}^{κ}) are known, and propose an algorithm to compute the new approximations (t_{1}^{κ+1}, …, t_{N}^{κ+1}) in A/I^{2^{κ+1}}[X_{1}, …, X_{N}]. This algorithm is based on Proposition 4 below; it finally amounts to linear algebra operations in a suitable quotient ring.

As input, we take any triangular set T_{κ} = (T_{1}^{κ}, …, T_{N}^{κ}) of polynomials in A/I^{2^{κ+1}}[X_{1}, …, X_{N}] such that:

• (H3) T_{j}^{κ} = t_{j}^{κ} modulo I^{2^{κ}} A/I^{2^{κ+1}}[X_{1}, …, X_{N}], for j in 1, …, N.

Stating the proper sequence of computations requires some new notations.

• We denote by H_{κ} the quotient A/I^{2^{κ+1}}[X_{1}, …, X_{N}]/(T_{1}^{κ}, …, T_{N}^{κ}). The canonical image of any α ∈ A[X_{1}, …, X_{N}] in H_{κ} is denoted by α_{κ}; similar notations hold for vector or matrices of polynomials.

• Jac(F) and jac(F) respectively denote the jacobian matrix of the system F and its determinant; Jac(t) denotes the jacobian matrix of t. Similar notation with subscript κ denotes their images in the matrix algebra over H_{κ}, and Jac(T_{κ}) denotes the jacobian matrix of T_{κ}. The identity matrix is denoted by I. F_{κ} denotes the reduction of the vector of polynomials F in H_{κ}.

• Since T_{κ} is a triangular set, the quotient H_{κ} is a free A/I^{2^{κ+1}}-module, which admits for a basis the set of monomials

{X_{1}^{α_{1}} … X_{N}^{α_{N}} , 0 ≤ α_{j} < deg_{X_{j}} T_{j}^{κ}}.

This canonical basis enables to assign to any element h in H_{κ} a canonical preimage in A/I^{2^{κ+1}}[X_{1}, …, X_{N}], denoted by h̃.

With these notations, the following proposition gives the formula for computing the next approximations (t_{1}^{κ+1}, …, t_{N}^{κ+1}).

Proposition 4 The Jacobian matrix Jac(F_{κ}) is invertible. Let δ_{κ} = (δ_{1}^{κ}, …, δ_{N}^{κ}) be the product Jac(T_{κ})Jac(F_{κ})^{−1}F_{κ}. Then, for all j in 1, …, N, the equality t_{j}^{κ+1} = T_{j}^{κ} + δ̃_{j}^{κ} holds in A/I^{2^{κ+1}}[X_{1}, …, X_{N}].

Proof. The Jacobian determinant jac(F) is invertible in A[X_{1}, …, X_{N}]/(I, t_{1}, …, t_{N}), so by Hensel’s Lemma its image is invertible in the quotient h_{κ} = A[X_{1}, …, X_{N}]/(I^{2^{κ}}, t_{1}, …, t_{N}). It is straightforward to check that h_{κ} ≃ H_{κ}/I^{2^{κ}}H_{κ}, so another application of Hensel’s Lemma shows that jac(F) is invertible in H_{κ}/I^{2^{κ+1}}H_{κ}. We note that I^{2^{κ+1}}H_{κ} = 0, so the previous quotient is H_{κ}. This proves the first point.

The second point is proven through the following explicit computations.

The equality F = At implies that the Jacobian matrix Jac(F) is equal to AJac(t)+B, where all entries in the matrix B belong to the ideal (t_{1}, …, t_{N}). The polynomials T_{j}^{κ} are chosen such that the images of the polynomials t in H_{κ} belong to the ideal I^{2^{κ}}H_{κ}. Consequently, in the relation Jac(F_{κ}) = A_{κ}Jac(t_{κ}) + B_{κ} over H_{κ}, all entries of B_{κ} belong to I^{2^{κ}}H_{κ}.

This formula implies that Jac(t_{κ}) is invertible over H_{κ}/I^{2^{κ}}H_{κ} and so, by a new application of Hensel’s lemma, over H_{κ} as well. Consequently, the equality

I = Jac(t_{κ})Jac(F_{κ})^{−1}A_{κ} + Jac(t_{κ})Jac(F_{κ})^{−1}B_{κ}Jac(t_{κ})^{−1}

holds, which can be rewritten

Jac(t_{κ})Jac(F_{κ})^{−1}A_{κ} = I + C_{κ},

where the entries of C_{κ} belong to I^{2^{κ}}H_{κ}. We then deduce the series of equalities

Jac(t_{κ})Jac(F_{κ})^{−1}F_{κ} = Jac(t_{κ})Jac(F_{κ})^{−1}A_{κ}t_{κ} = t_{κ} + C_{κ}t_{κ} = t_{κ},

since all entries in t_{κ} and C_{κ} belong to I^{2^{κ}}H_{κ}, and I^{2^{κ+1}}H_{κ} = 0.

In a similar way, the entries in the matrix Jac(t_{κ}) differ from the entries in Jac(T_{κ}) by elements in I^{2^{κ}}H_{κ}. Since F_{κ} is in I^{2^{κ}}H_{κ}, this implies the equality over H_{κ}:

Jac(T_{κ})Jac(F_{κ})^{−1}F_{κ} = Jac(t_{κ})Jac(F_{κ})^{−1}F_{κ} = t_{κ}.  (1)

Let δ_{κ} be the vector Jac(T_{κ})Jac(F_{κ})^{−1}F_{κ}, computed over H_{κ}, and δ̃_{κ} the vector of the canonical preimages in A/I^{2^{κ+1}}[X_{1}, …, X_{N}] of its entries.

Equation (1) means that δ̃_{j}^{κ} − t_{j}^{κ+1} belongs to (T_{1}^{κ}, …, T_{N}^{κ}). Consequently, δ̃_{j}^{κ} − t_{j}^{κ+1} + T_{j}^{κ} also belongs to this ideal, and has partial degree in every X_{k} less than deg_{X_{k}} T_{k}^{κ}. Since (T_{1}^{κ}, …, T_{N}^{κ}) forms a triangular set, this implies that δ̃_{j}^{κ} − t_{j}^{κ+1} + T_{j}^{κ} is zero. This concludes the proof. □

Algorithm Lift. The previous proposition is turned into the following procedure Lift in a straightforward way; the notations used in this algorithm are the same as in the previous paragraph.

Symbolic Newton lifting

Procedure Lift(F, T)

Input: the system F, a triangular set T = (T_{1}^{κ}, …, T_{N}^{κ}), which satisfy hypotheses (H1), (H2) and (H3). Output: the polynomials (t_{1}, …, t_{N}) modulo I^{2^{κ+1}}.

# The computations are done over A/I^{2^{κ+1}}[X_{1}, …, X_{N}]/(T_{1}^{κ}, …, T_{N}^{κ})

    Jac(F_{κ}) ← JacobianMatrix(F_{κ});

Jac(F_{κ})^{−1} ← Inverse(Jac(F_{κ}));

    Jac(T_{κ}) ← JacobianMatrix(T_{κ});

δ_{κ} ← Jac(T_{κ})Jac(F_{κ})^{−1}F_{κ}; return (T_{1}^{κ} + δ̃_{1}^{κ}, …, T_{N}^{κ} + δ̃_{N}^{κ});

Application to parametric geometric resolutions. To conclude this subsection, we present the application of this method to our problem of parametric resolutions. The notations are those used in the rest of this paper.

Let p be a point in k^{m}, and u = ∑ u_{i}x_{i} a primitive element of K → K[X_{1}, …, X_{n}]/J_{K}. We consider the Shape Lemma parametrization S_{u} = [Q_{u}, W_{1}, …, W_{n}] defined in Section 3.3 and suppose that p cancels no denominator in S_{u}. In this case, we denote by S_{u}^{κ} = [Q_{u}^{κ}, W_{1}^{κ}, …, W_{n}^{κ}] the vector S_{u}, where all coefficients are replaced by their Taylor expansion at p at precision 2^{κ}.

The key point is that a resolution under Shape Lemma form is a particular form of triangular set. Then Proposition 4 enables us to compute S_{u}^{κ+1} from S_{u}^{κ}. The idea is that computing Taylor expansions at p amounts to compute modulo the powers of the maximal ideal of the m-variate power series ring centered at p.

The proposition below justifies this assertion, and estimates the complexity of the process.

From the proof of this proposition, we deduce that computing the new approximation S_{u}^{κ+1} amounts to calling the procedure Lift defined above, with arguments (f_{1}, …, f_{n}, X_{n+1} − ∑ u_{i}X_{i}) and (X_{1} − W_{1}^{κ}(X_{n+1}), …, X_{n} − W_{n}^{κ}(X_{n+1}), Q_{u}^{κ}(X_{n+1})). We will call Lift(f, S_{u}^{κ}, u) this process; its output is S_{u}^{κ+1}.

Proposition 5 Let u = ∑ u_{i}x_{i} be a primitive element of K → K[X_{1}, …, X_{n}]/J_{K} with coefficients in k, and S_{u} = [Q_{u}, W_{1}, …, W_{n}] the vector corresponding to the parametrization

Q_{u}(u) = 0, { x_{1} = W_{1}(u), ; ⋮ ; x_{n} = W_{n}(u).

Let p be a point in k^{m} which cancels none of the denominators in S_{u}, such that the points described by the specialization of S_{u} at p do not cancel the Jacobian determinant jac(f, X). Given the approximation S_{u}^{κ}, S_{u}^{κ+1} can be computed in

O((nL + n^{4})M_{u}(deg_{π})M_{s}(2^{κ+1}, m))

operations in k.

Proof. We get back to the previous setting by introducing X = (X_{1}, …, X_{n}, X_{n+1}), the triangular set t = (X_{1} − W_{1}(X_{n+1}), …, X_{n} − W_{n}(X_{n+1}), Q_{u}(X_{n+1})) and the equations F = (f_{1}, …, f_{n}, X_{n+1} − ∑ u_{i}X_{i}). There exists a (n + 1) × (n + 1)-matrix A with entries in k(P_{1}, …, P_{m})[X_{1}, …, X_{n}, X_{n+1}] such that F = At, where p cancels none of the denominators in this equality. Let (P′_{1}, …, P′_{m}) be new variables, and A the power series ring k[[P′_{1}, …, P′_{m}]]. All the entries in the previous matrix equality can be rewritten in terms of (P′_{1}, …, P′_{m}), letting P_{j} = P′_{j} + p_{j}, for j = 1, …, m. None of the new denominators vanishes at zero, so all entries admit Taylor expansion in A.

We apply Proposition 4 with I = (P′_{1}, …, P′_{m}). Using the notation of the previous paragraph, the quotient H_{κ} is isomorphic to k[[P′_{1}, …, P′_{m}]]/I^{2^{κ+1}}[U]/Q_{u}^{κ}. Since Q_{u} has degree deg_{π}, and the coefficients are power series in m variables at precision 2^{κ+1}, the cost of an operation in H_{κ} is O(M_{u}(deg_{π})M_{s}(2^{κ+1}, m)) operations in k.

The cost of the computations lies in the evaluation of the vectors and matrices involved, and in the linear algebra operations. Evaluating the non-zero terms in the matrix Jac(T_{κ}) takes O(n deg_{π} M_{s}(2^{κ+1}, m)) operations in k. Using Baur-Strassen’s algorithm [6], evaluating the matrix Jac(F_{κ}) and the vector F_{κ} takes O(nL) operations in the quotient H_{κ}. All linear algebra takes O(n^{4}) operations in H_{κ}, using for instance Leverrier’s algorithm [43] for matrix inversion over a ring. All this sums up to O((nL + n^{4})M_{u}(deg_{π})M_{s}(2^{κ+1}, m)) operations in k, which proves the proposition. □

4.3 Recovering the coefficients

Referring to the algorithm sketched in the introduction of this section, the last question to answer is how to recover the coefficients of a parametric resolution from their Taylor expansion at some point p. To this effect, we present an algorithm for the reconstruction of a rational function.

The problem can be stated as follows: let D be a positive integer, p and q two polynomials in k[P_{1}, …, P_{m}] of degrees at most D such that q(0) ≠ 0, and r the Taylor expansion of p/q at precision 2D + 1. Given r as a polynomial of degree 2D, we want to compute p/q.

In the single variable case, i.e. when m = 1, this question is solved using Padé approximants, see [64]. In our general multivariate case, the question can be solved using linear algebra; other solutions based on Gröbner bases computations are presented in [54, 22]. We propose an algorithm with better complexity, which reduces to the usual computation of Padé approximants when m = 1.

The algorithm is probabilistic: it requires to choose m − 1 values in the base field. We indicate the degree of an hypersurface in 𝔸^{m−1}(\overline{k}) that must be avoided to ensure success.

The first paragraph is devoted to present the algorithm. We then estimate the cost of applying it on all the coefficients of a parametric resolution. Due to our approximation process, the origin of the coordinates in the parameter space has moved; to conclude this subsection, we consider the question of restoring the initial coordinates.

4.3.1 Rational reconstruction

The main idea of our algorithm is to get back to an Euclidean situation: we introduce a new variable s, and substitute the variables P_{i} by P_{i}s in r. Then we perform a single-variable Padé approximation with main variable s, from which we recover the fraction p/q.

This algorithm uses the following subroutines, where [P′_{2}, …, P′_{m}] are new variables.

• PadeApproximant(r̃), where r̃ is a polynomial of degree 2D in k[[P′_{2}, …, P′_{m}]][s], with coefficients of precision D.

The function follows the algorithms given in [64, chapters 5.9 and 11] to compute (D, D) Padé approximant of r̃. This is done by applying the fast extended monic Euclidean algorithm to r̃ and s^{2D+1}. The coefficients are power series of fixed precision D; the function raises an error if a division by a series of positive valuation occurs.

• ConstantCoefficient(\overline{q}), where \overline{q} is in k[[P′_{2}, …, P′_{m}]][s].

The function returns the coefficient of degree 0 of \overline{q}, if this coefficient is not zero. Else, it raises an error.

• Homogenization(p̃, P_{1}), where p̃ belongs to k[P_{2}, …, P_{m}][s].

Suppose that for all i the coefficient of s^{i} in p̃ has degree at most i. Then, for all i, this function homogenizes the coefficient of s^{i} in degree i with respect to the variable P_{1}. If the previous assumption is not satisfied, an error is raised.

Here is our algorithm, which chooses m − 1 values in the base field. The subsequent proposition shows that the output is correct for a generic choice, gives a bound on the complexity of the process and the probability of success.

Rational reconstruction

Procedure RationalReconstruction(r, γ) Input: r in k[P_{1}, …, P_{m}] of degree 2D, γ = (γ_{2}, …, γ_{m}) in k^{m−1}. Output: a fraction p/q or “Failure”.

r̃ ← r(s, P_{2}s, …, P_{m}s); # We change the coordinates. r̃ ← subs(P_{2} = P′_{2} + γ_{2}, …, P_{m} = P′_{m} + γ_{m}, r̃); # The computations are done in k[[P′_{2}, …, P′_{m}]][s] # with coefficients truncated at precision D.

\overline{p}, \overline{q} ← PadeApproximant(r̃);
p̃ ← \overline{p}/ConstantCoefficient(\overline{q});
q̃ ← \overline{q}/ConstantCoefficient(\overline{q});

# We change back the coordinates. p̃ ← subs(P′_{2} = P_{2} − γ_{2}, …, P′_{m} = P_{m} − γ_{m}, p̃); q̃ ← subs(P′_{2} = P_{2} − γ_{2}, …, P′_{m} = P_{m} − γ_{m}, q̃);

P̃ ← Homogenization(p̃, P_{1});
Q̃ ← Homogenization(q̃, P_{1});
return subs(s = 1, P̃/Q̃);

Proposition 6 Suppose that there exist (p, q) of degrees at most D, such that r is the Taylor expansion of p/q at precision 2D + 1. For almost all choices of (γ_{2}, …, γ_{m}), the previous algorithm computes p/q using

O_{log} (M_{u}(D) (M_{s}(D, m − 1) + m^{2}\binom{2D + m − 1}{m − 1})) ⊂ O_{log} (m^{2}M_{u}(D)M_{s}(2D, m − 1))

operations in k. The choices of (γ_{2}, …, γ_{m}) that lead to an error belong to an hypersurface of 𝔸^{m−1}(\overline{k}) of degree at most 2D(2D + 1)^{2}.

The proof is divided in several steps. We first show how using the new main variable s reduces the problem to univariate Padé approximant computations over a rational function field. Then we show how to replace rational function computations by power series computations; from this we deduce the complexity and error analysis.

Introduction of a new main variable. Let R̃ be the substitution r(P_{1}s, …, P_{m}s), for a new variable s. R̃ belongs to k(P_{1}, …, P_{m})[s] and has degree 2D; we denote by (\overline{P} , \overline{Q}) its (D, D) Padé approximant, computed in k(P_{1}, …, P_{m})[s] by applying the fast Euclidean algorithm to R̃ and s^{2D+1}.

The polynomials (p(P_{1}s, …, P_{m}s), q(P_{1}s, …, P_{m}s)) also form a (D, D) Padé approximant of R̃, so uniqueness shows that they differ of (\overline{P} , \overline{Q}) by a factor in k(P_{1}, …, P_{m}). Consequently, dividing \overline{P} and \overline{Q} by the constant coefficient of \overline{Q} yields P̃ = p(P_{1}s, …, P_{m}s)/q(0) and Q̃ = q(P_{1}s, …, P_{m}s)/q(0), and the substitution s = 1 in P̃/Q̃ gives the requested output p/q.

Computing the Padé approximant of R̃ will not lead to an algorithm with good complexity: the coefficients that appear during the computations are rational functions, with increasing degrees. We now show how to replace these rational functions by power series, and deduce the proof of the proposition.

Dehomogenization. The coefficients that appear in the course of Euclid’s algorithm applied to R̃ and s^{2D+1} are homogeneous rational functions. Thus, if we want to introduce power series, is necessary to dehomogenize these coefficients. We now proceed to do so.

Let r̃ be r(s, P_{2}s, …, P_{m}s), and (\overline{p}, \overline{q}) the output of the Padé approximant computation applied to r̃. Since the coefficients that occur during Euclid’s algorithm applied to R̃ and s^{2D+1} are homogeneous, (\overline{p}, \overline{q}) coincide with the dehomogenization of (\overline{P} , \overline{Q}) in the variable P_{1}.

As above, we define (p̃, q̃) as (\overline{p}, \overline{q}) divided by the constant coefficient of \overline{q}. Then the previous remarks show that they coincide with the dehomogenization of (P̃, Q̃) in the variable P_{1}. Since for all i, the coefficients of s^{i} in P̃ and Q̃ are homogeneous of degree i, applying the subroutine Homogenization(., P_{1}) to p̃ and q̃ yields P̃ and Q̃.

Consequently, it suffices to compute (\overline{p}, \overline{q}) then (p̃, q̃) and finally (P̃, Q̃) to solve the problem.

Computation with power series. If none of the denominators occurring during Euclid’s algorithm applied to (r̃, s^{2D+1}) vanishes at zero, all the operations on the coefficients can be done at fixed precision D, replacing the rational functions in (P_{2}, …, P_{m}) by their power series expansion at order D.

To get back to this lucky situation, we perform a linear change of variables. Given a value (γ_{2}, …, γ_{m}) in k^{m−1}, we do the computations in the variables (P′_{2}, …, P′_{m}), where P_{i} = P′_{i}+γ_{i} for i = 2, …, m, and get back to the initial variables afterwards.

If A is a polynomial in k[P_{2}, …, P_{m}], rewritten A′ in the variables (P′_{2}, …, P′_{m}), the constant term in A′(P′_{2}, …, P′_{m}) is A(γ_{2}, …, γ_{m}). Consequently, the Padé approximant computations can be done at fixed precision in the new variables if (γ_{2}, …, γ_{m}) cancels none of the denominators that appear in the course of Euclid’s algorithm applied to r̃ and s^{2D+1}.

From [10, 64], this is the case if and only if (γ_{2}, …, γ_{m}) cancels none of the non-zero subresultant coefficients associated to r̃ and s^{2D+1}. These subresultants are minors of the Sylvester matrix associated to r̃ and s^{2D+1} of sizes i = 1, 3, …, 4D + 1, with entries of degrees at most 2D. Each determinant is a polynomial in k[P_{2}, …, P_{m}] of degree at most 2Di, so their product has degree at most 2D(2D + 1)^{2}. This proves the last point of the proposition.

Complexity. We now estimate the complexity. Since the coefficients are power series in m − 1 variables at precision D, Euclid’s algorithm takes O_{log}(M_{u}(D)M_{s}(D, m − 1)) operations in k, see [64, chapter 11].

The changes of variables require the translation by vectors (γ_{2}, …, γ_{m}) and (−γ_{2}, …, −γ_{m}) on the coefficients of r̃, p̃ and q̃. These polynomials have degree at most 2D, and their coefficients are polynomials in (P_{2}, …, P_{m}) of degrees at most 2D. The following lemma shows that the cost of a translation is O_{log}(mM_{u}(D)\binom{2D+m−2}{m−2}), so the sum of these costs is in O_{log}(mDM_{u}(D)\binom{2D+m−2}{m−2}), which is in O_{log}(m^{2}M_{u}(D)\binom{2D+m−1}{m−1}).

The proof of the proposition is now almost complete. We only have to establish the following lemma, which we used above. It gives the cost of translating the variables in a multivariate polynomial.

Lemma 6 Let A a polynomial in k[P_{2}, …, P_{m}] of degree D, and (γ_{2}, …, γ_{m}) a point in k^{m−1}. Then A(P_{2} + γ_{2}, …, P_{m} + γ_{m}) can be computed in O_{log}(mM_{u}(D)\binom{D+m−2}{m−2}) operations in k.

Proof. We move one variable P_{i} at a time; to keep the notation simple, we describe the case i = 2. We consider the polynomial A in k[P_{2}][P_{3}, …, P_{m}]; then the translation is done by shifting all coefficients.

Using the divide-and-conquer algorithm given in [64, chapter 9.2], shifting a single coefficient requires O(M_{u}(D) log D) operations in k. Since there are at most \binom{D+m−2}{m−2} such coefficients, this sums up to O(M_{u}(D) log D\binom{D+m−2}{m−2}) operations in k. Taking all variables P_{i} into account leads to the announced bound. □

4.3.2 Application to parametric resolutions

The rational reconstruction process will be applied to all the coefficients of a parametric resolution. If R = [Q_{u}, V_{1}, …, V_{n}] is a vector of polynomials in k[[P_{1}, …, P_{m}]][U] and γ a point in k^{m−1}, we denote by RationalReconstruction(R, γ) the application of the reconstruction process to all the coefficients of the polynomials in R. The output is a boolean value b which indicates success, and, if possible, a sequence of polynomials, where all coefficients are reconstructed.

Due to our Newton approximation scheme, the Taylor expansions will be given at precisions of the form 2^{κ}. In this short paragraph, we indicate the total cost of the reconstruction under such constraints. Recall that the reconstruction requires to choose m − 1 values in the base field: we also indicate the degree of a hypersurface of 𝔸^{m−1}(\overline{k}) that must be avoided to ensure success.

We suppose that R corresponds to a parametric resolution in Kronecker form. Then we can suppose that all coefficients in the parametric resolution have for maximal degree an integer denoted by deg_{u}. Using Theorem 1, deg_{u} is bounded by the geometric degree deg_{V}, itself bounded by the Bézout number d^{n}. We will also use the hypothesis that the polynomials [Q_{u}, V_{1}, …, V_{n}] have degree at most deg_{π} ≤ d^{n} in their main variable U.

If p/q is a fraction with numerator and denominator of degrees bounded by deg_{u}, then the first power of 2 that permits reconstruction is the first power of 2 greater than 2 deg_{u} +1. If we denote by ⌈x⌉ the first integer greater than or equal to x, this power is 2^{κ_{0}}, where κ_{0} = ⌈log_{2}(2 deg_{u} +1)⌉. Since deg_{u} ≤ d^{n}, then 2^{κ_{0}} ≤ 4d^{n}.

Applying Proposition 6 with 2D = 2^{κ_{0}} shows that a single coefficient can be reconstructed if the change of variables γ avoids an hypersurface of degree at most 4d^{n}(4d^{n} + 1)^{2}, and the reconstruction then takes O_{log} (m^{2}M_{u}(2^{κ_{0}})M_{s}(2^{κ_{0}}, m − 1)) operations in k. Taking all coefficients into account then leads to O_{log}(nm^{2} deg_{π} M_{u}(2^{κ_{0}})M_{s}(2^{κ_{0}}, m − 1)) operations in k.

Since all polynomials in R have degree in U at most deg_{π} and Q_{u} is monic, there are at most (n + 1)d^{n} rational function to recover, and all of them can be recovered if γ avoids the union of all corresponding hypersurfaces. This union has degree at most 4(n + 1)d^{2n}(4d^{n} + 1)^{2}.

Summary. Let us summarize the results we will need in the sequel:

• We assume that the coefficients to reconstruct are rational functions with numerators and denominators of degree at most deg_{u} ≤ d^{n}.

• The first power of 2 that enables the reconstruction is 2^{κ_{0}}, with κ_{0} = ⌈log_{2}(2 deg_{u} +1)⌉.

• The total cost of the reconstruction is within O_{log}(nm^{2} deg_{π} M_{u}(2^{κ_{0}})M_{s}(2^{κ_{0}}, m − 1)) operations in k.

• The hypersurface of 𝔸^{m−1}(\overline{k}) to avoid has degree at most 4(n + 1)d^{2n}(4d^{n} + 1)^{2}.

4.3.3 Going back to the initial coordinates

Suppose that the reconstruction of the parametric resolution is successful. Due to our approximation process, the coordinates in the parameter space are centered at some value (p_{1}, …, p_{m}) in k^{m}, so we must move back to the initial coordinates. Given a resolution R = [Q_{u}, V_{1}, …, V_{n}] as a vector of polynomials in k(P_{1}, …, P_{m})[U], we introduce a subroutine RestoreCoordinates(R) devoted to this operation.

Using the same notation as above, we suppose that the coefficients of the polynomials in R have degree in P_{1}, …, P_{m} at most 2^{κ_{0}−1}. Lemma 6 shows that the cost necessary to move a single coefficient is O_{log}(mM_{u}(2^{κ_{0}−1})\binom{2^{κ_{0}−1}+m−1}{m−1}) operations in k. If we suppose that the polynomials in R have degree at most deg_{π} in U, then there are (n + 1) deg_{π} coefficients to move. The subroutine RestoreCoordinates then induces a total cost of O_{log}(nm deg_{π} M_{u}(2^{κ_{0}−1})\binom{2^{κ_{0}−1}+m−1}{m−1}) operations in k.

5 The main algorithm

The main algorithm consists in a loop organized around the Newton approximation process; tests are performed at each pass to decide whether to stop the computation or not.

Before giving the details of the main algorithm, we present the additional subroutine, denoted StopCriterion, which decides whether to stop the lifting. In a second time, we give the whole algorithm, and work out its complexity and probability of success. Finally, we mention some possible practical improvements.

We constantly switch between the two forms of parametrization, Kronecker and Shape Lemma: the former has better degree properties, so is well suited for the rational reconstruction, whereas Newton’s iterator works using the later representation. In the sequel, we denote KroneckerParametrization(S_{u}) a routine which, given a Shape Lemma parametrization, outputs the corresponding Kronecker parametrization, and ShapeLemmaParametrization(R_{u}) the converse process.

5.1 The stop criterion

For obvious practical reasons, we do not perform the lifting up to the Bézout bound. Instead, we use a probabilistic test, presented in the subroutine StopCriterion: once we have a candidate resolution, the test mainly consists in testing it on a witness point p′; there is a possibility of choosing a bad witness, which will be taken into account in the proof of Proposition 7.

Stop criterion for the lifting process

Procedure StopCriterion(R, f, p′)

Input: a parametric resolution R = [Q_{u}, V_{1}, …, V_{n}], the system f, p′ in k^{m}.

Output: a boolean value.

    if

p′ cancels none of the denominators in R, the specialization of Q_{u} at p′ is squarefree, the points described by the specialization of R at p′ cancel the system f(p′, .), these points do not cancel the Jacobian determinant of f(p′, .),

    then return true
    else return false

Lemma 7 Suppose that the polynomials in R = [Q_{u}, V_{1}, …, V_{n}] have degree at most deg_{π}, and that their coefficients are rational functions of degree at most D. Then the cost of the subroutine StopCriterion is O_{log}(n deg_{π} \binom{D+m}{m} + (nL + n^{4})M_{u}(deg_{π})) operations in k.

Proof. All the (n + 1) deg_{π} coefficients of R have degree at most D in m variables, so their specialization on the point p′ takes less than (n + 1) deg_{π} \binom{D+m}{m} operations in k. Let us denote [q_{u}, v_{1}, …, v_{n}] the specialized resolution, with coefficients in k.

Testing whether the points described by this specialization cancel f requires to switch to the equivalent Shape Lemma Parametrization [q_{u}, w_{1}, …, w_{n}], with coefficients in k, then to evaluate f on the elements [w_{1}, …, w_{n}] modulo q_{u}. The first task requires to invert q′_{u} modulo q_{u}, hence has cost O_{log}(M_{u}(deg_{π})), using the fast Euclidean algorithm [64]. The second task takes O(LM_{u}(deg_{π})) additional operations.

Similarly, the Jacobian matrix can be evaluated in nL operations modulo q_{u} using Baur-Strassen’s algorithm [6] and its determinant can be computed in n^{4} operations modulo q_{u}, which adds O((nL + n^{4})M_{u}(deg_{π})) operations. Its invertibility can be tested within O_{log}(M_{u}(deg_{π})) operations. This yields the overall complexity bound. □

5.2 The detailed algorithm

We are now ready to present the main algorithm ParametricResolution. It chooses 3m−1 values in the base field: the coordinates of the points p and p′, and the change of variables γ used in the rational reconstruction. The following proposition shows that for a generic choice, the output is correct, and quantifies the bad choices. This brings the proof of Theorem 2 restated in the proposition below.

Computing a parametric resolution

Procedure ParametricResolution(f) Input: the system f = (f_{1}, …, f_{n}). Output: a parametric geometric resolution or “Failure”.

p, p′ ← points in k^{m}; γ ← point in k^{m−1}; u, R_{u}^{0} ← Resolution(f, p);

    S_{u}^{0} ← ShapeLemmaParametrization(R_{u}^{0});

MaxSteps ← ⌈log_{2}(2d^{n} + 1)⌉; κ ← 0;

while κ ≤ MaxSteps do

      R_{u}^{κ} ← KroneckerParametrization(S_{u}^{κ});
      b, R_{u}^{κ} ← RationalReconstruction(R_{u}^{κ}, γ);
      if b then

finished ← StopCriterion(R_{u}^{κ}, f, p′);

        if finished then return RestoreCoordinates(R_{u}^{κ});
      end if;

S_{u}^{κ+1} ← Lift(f, S_{u}^{κ}, u);

κ ← κ + 1;

    end while;

return “Failure”; Proposition 7 Let Γ be a subset of k, and suppose that (p, p′, γ) are chosen in Γ^{3m−1}. If n ≥ 2 and d ≥ 2, then the algorithm ParametricResolution computes a parametric resolution of K → K[X_{1}, …, X_{n}]/J_{K} for all choices except at most 110nd^{4n}|Γ|^{3m−2}.

In case of success, the complexity of the lifting step is

O_{log}((nL + n^{4})M_{u}(deg_{π})M_{s}(4 deg_{u}, m) + nm^{2} deg_{π} M_{u}(deg_{u})M_{s}(4 deg_{u}, m − 1))

operations in k, where deg_{u} is the maximum of the degrees in P_{1}, …, P_{m} of the coefficients that appear in the parametric resolution. Else, the algorithm stops after at most ⌈log_{2}(2d^{n} + 1)⌉ lifting steps, and outputs either “Failure” or a wrong answer.

We first show that the algorithm computes the correct answer for generic choices of (p, p′, γ); we return to the quantification of the bad choices in a second time, then estimate the complexity of the process.

Proof of correctness. Let Δ be the polynomial defined in Proposition 3; we assume that the polynomial Δ(p, .) is not zero, and let u, R_{u}^{0} be the output of Resolution(f, p). The second part of Proposition 3 shows that Δ(p, u) is not zero, and that R_{u}^{0} is the specialization of a generic resolution R_{u} at p.

The points p and u satisfy the hypotheses of Proposition 5 so after κ lifting steps, the coefficients of R_{u} are known at precision 2^{κ}. Let deg_{u} be the maximum of the degrees in P_{1}, …, P_{m} of the coefficients that appear in R_{u}; Proposition 2 shows that deg_{u} ≤ d^{n}. We call κ_{0} the number of lifting steps necessary before the reconstruction of all the coefficients is possible; the results of Subsection 4.3.2 show that κ_{0} ≤ ⌈log_{2}(2 deg_{u} +1)⌉ ≤ ⌈log_{2}(2d^{n} + 1)⌉.

We now rule out the possibility that, for some κ ≤ κ_{0} − 1, the rational reconstruction of R_{u}^{κ} ≠ R_{u} is possible and the subroutine StopCriterion(R_{u}^{κ}, p′, f) outputs true; in this case p′ will be called a bad witness.

Since Δ(p, u) is not zero, the polynomial Δ(., u) itself is not zero. We suppose that p′ does not cancel this polynomial; the first point in Proposition 3 then implies that the simple solutions of the specialized system f(p′, .) = 0 are described by the specialization of R_{u} at p′.

The subroutine StopCriterion outputs true at step κ ≤ κ_{0} −1 if the simple solutions of the system f(p′, .) = 0 are described by the specialization of R_{u}^{κ} at p′, i.e. if the specializations of R_{u} and R_{u}^{κ} coincide at p′. Since R_{u} and R_{u}^{κ} are different, at least one of their coefficients differs. These coefficients are rational functions of degrees at most d^{n} and 2^{κ−1}, so the points where their specializations coincide are contained in an hypersurface of 𝔸^{m}(\overline{k}) of degree at most d^{n} + 2^{κ−1}.

Taking all κ < κ_{0} into consideration shows that the point p′ is a “good witness” if it avoids an hypersurface 𝔸^{m}(\overline{k}) of degree at most d^{n}+0+d^{n}+1+⋯+d^{n}+2^{(κ_{0}−1)−1} ≤ d^{n}κ_{0}+2^{κ_{0}−1} ≤ d^{n}(⌈log_{2}(2d^{n} + 1)⌉ + 2). We suppose that this is the case.

The algorithm can then fail only if the reconstruction at step κ_{0} fails, so success is assured if the change of variables γ avoids the hypersurface defined in Subsection 4.3.2; this hypersurface has degree at most 4(n + 1)d^{2n}(4d^{n} + 1)^{2}.

Estimation of probabilities. We now return to the assumptions made on (p, p′, γ) and use Zippel-Schwartz’ lemma [67, 60] to quantify the choices that assure success. Let Γ be a subset of k, and suppose that the values (p, p′, γ) are chosen in Γ^{m} × Γ^{m} × Γ^{m−1}. Besides, we recall that the polynomial Δ has degree at most d^{n}(2d^{n} + nd + 1) in P_{1}, …, P_{m}.

• There at most d^{n}(2d^{n} + nd + 1)|Γ|^{m−1} values of p such that Δ(p, .) = 0; this discriminates at most d^{n}(2d^{n} + nd + 1)|Γ|^{3m−2} choices of (p, p′, γ).

• For all remaining values of p, there are at most 4(n + 1)d^{2n}(4d^{n} + 1)^{2}|Γ|^{m−2} values of γ which prevent the reconstruction; this represents at most 4(n + 1)d^{2n}(4d^{n} + 1)^{2}|Γ|^{3m−2} choices of (p, p′, γ). If m = 1, the reconstruction is deterministic, so this possibility of failure is not taken into account.

• For all remaining values of p, for any value of u, there are at most d^{n}(2d^{n}+nd+1)|Γ|^{m−1} values of p′ such that Δ(p′, u) = 0; this discriminates at most d^{n}(2d^{n} + nd + 1)|Γ|^{3m−2} choices of (p, p′, γ).

• Finally, for any value of p and u, there are at most d^{n}(⌈log_{2}(2d^{n} +1)⌉+2)|Γ|^{m−1} values of p′ which are bad witnesses. This represents at most d^{n}(⌈log_{2}(2d^{n} + 1)⌉ + 2)|Γ|^{3m−2} triples (p, p′, γ).

The number of bad choices is thus at most

d^{n}(4(n + 1)d^{n}(4d^{n} + 1)^{2} + 4d^{n} + 2nd + ⌈log_{2}(2d^{n} + 1)⌉ + 4)|Γ|^{3m−2}.

Using the rough estimates log(1 + x) ≤ x and nd ≤ d^{n} if d ≥ 2, this quantity is seen to be bounded by

d^{2n}(4(n + 1)(16d^{2n} + 8d^{n} + 1) + 12)|Γ|^{3m−2}.

Using n + 1 ≤ 3n/2 for n ≥ 2, we bound this number by nd^{2n}(96d^{2n} + 48d^{n} + 18)|Γ|^{3m−2}, which itself is bounded by 110nd^{4n}|Γ|^{3m−2}.

Complexity. We finally turn to the complexity of the algorithm, and detail the cost of the last call to each subroutine, in terms of operations in k.

• The last call to Lift brings the precision to 2^{κ_{0}}. Proposition 5 shows that its cost is in O((nL + n^{4})M_{u}(deg_{π})M_{s}(2^{κ_{0}}, m)).

• The subroutine KroneckerParametrization requires to multiply all parametrizations by the derivative of the minimal polynomial Q′_{u}, and to reduce them modulo Q_{u}. This takes O(nM_{u}(deg_{π})M_{s}(2^{κ_{0}}, m)) operations in k.

• We saw in Subsection 4.3.2 that the total cost of the subroutine RationalReconstruction is in O_{log}(nm^{2} deg_{π} M_{u}(2^{κ_{0}})M_{s}(2^{κ_{0}}, m − 1)).

• Lemma 7 shows that the complexity of StopCriterion is in O_{log}(n deg_{π} \binom{2^{κ_{0}−1}+m}{m} + (nL + n^{4})M_{u}(deg_{π})). In case of success, Subsection 4.3.3 shows that the subroutine RestoreCoordinates has complexity O_{log}(nm deg_{π} M_{u}(2^{κ_{0}−1})\binom{2^{κ_{0}−1}+m−1}{m−1}).

All these costs sum up to

O_{log}((nL + n^{4})M_{u}(deg_{π})M_{s}(2^{κ_{0}}, m) + nm^{2} deg_{π} M_{u}(2^{κ_{0}})M_{s}(2^{κ_{0}}, m − 1)).

Under our assumption that M_{s}(d, m) ≤ cM_{s}(2d, m) for some universal constant c < 1, the cost of all steps is bounded by 1/(1 − c) times the cost of the last step. Since 2^{κ_{0}} ≤ 4 deg_{u}, this yields the overall complexity bound. □

5.3 Practical strategies

This section presents possible improvements of the main algorithm, which have important practical impact.

Modular arithmetic. Most systems we present as applications are defined over the rational field. To avoid the growth of the intermediate coefficients, it is natural to adopt a strategy based on modular computation: the resolution is first computed modulo some prime number 𝔭, then lifted modulo the successive powers 𝔭^{2^{κ}}, and the rational numbers are recovered when their 𝔭-adic approximation is precise enough.

This process is quite similar to the lifting of the parameters presented here, and is used in practice. All necessary algorithmic tools are given in this paper, except the reconstruction of rational numbers, which is a well-solved problem [62, 15, 64]. Still, we do not give more details on this question: such a strategy induces a variety of new possibilities of failure, whose analysis requires to use arithmetic versions of Bézout’s theorem and of the Nullstellensatz. This is beyond the scope of this paper; such results may be found in the author’s PhD. Thesis [58].

Factorization. Suppose that the minimal polynomial of the specialized system splits into i irreducible factors of degrees (deg_{π}^{(1)}, …, deg_{π}^{(i)}), so that this fiber can be described by i geometric resolutions of smaller degree. 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 RationalReconstruction and StopCriterion.

The term M_{u}(deg_{π}) = M_{u}(∑ deg_{π}^{(i)}) in the complexity of the lifting step is replaced by ∑ M_{u}(deg_{π}^{(i)}). This is most important for a naive multiplication algorithm where M_{u}(D) = O(D^{2}), or Karatsuba’s method, for which M_{u}(D) = O(D^{1.59}). The cost of the recombination of all factors is analyzed in [41], and does not modify our complexity bound.

Computing outside a given hypersurface. Suppose that we want the points V′ ⊂ V lying outside a given hypersurface H ⊂ 𝔸^{m+n}(\overline{k}). It suffices to remove the points of H that intersect the specialized fiber, and perform the lifting on what remains. The correctness of this process requires that the specialization value p avoids the projection of \overline{V}′ ∩H on 𝔸^{m}(\overline{k}); taking this possibility into account does not modify the rough upper bound 110nd^{4n} in the quantification of the probability of success presented above.

6 Applications

This last section gathers some applications which where the initial motivation for the design of our algorithm. All systems are not displayed for lack of space; the equations are available upon request.

The algorithm is implemented in Magma [2]. The Kronecker package developed by G. Lecerf [40] provided many necessary functionalities. We compared our timings with Gröbner Bases computations, for we could easily find a primitive element in each example; we used Magma for these computation, as it allows for such computations on rational function fields.

6.1 Description of the systems

Number of points of a Jacobian. This example, denoted P19 in section 6.2, describes computations that were performed with P. Gaudry and R. Harley for their genus 2 point counting record [23].

More precisely, the framework is the determination of the number of points of the Jacobian of a curve of genus 2 defined over the finite field k = 𝔽_{𝔭}, where 𝔭 is the first prime greater than 10^{19}. Following Schoof’s algorithm for elliptic curves [57], their algorithm is based on explicit computation of divisors of ℓ^{i}-torsion, for various primes ℓ. This part describes the computations for the case ℓ = 2.

We have considered a polynomial system with coefficients in 𝔽_{𝔭} whose resolution gives some divisors D_{i+1} of 2^{i+1}-torsion from the knowledge of a divisor D_{i} of 2^{i}-torsion: this system thus encodes the halving in the Jacobian.

This polynomial system has 4 equations in 6 variables, which split in the 4 coordinates of D_{i+1} plus 2 parameters that are functions of D_{i}. Given a 2-torsion divisor D_{1}, the resolution of this system gives a 4-torsion divisor D_{2}, which is in turn fed to a similar system, and so on. The objective is to go as far as possible, to refine the knowledge of the cardinality.

As i increases, the divisors D_{i} have their coordinates in extensions of k of increasing degrees, so the resolution of the specialized systems gets harder. It is preferable to compute the parametric resolution for once, then to specialize it when needed.

A priori considerations show that this system has generically 64 solutions, and the resolution of a specialized system shows that they are all simple, so our algorithm applies. Our output is only generically valid, but the specialization values caused no problem. We were thus able to compute divisors up to 256-torsion.

Deformation of singular hypersurfaces. The theoretical background for this problem can be found in the article by F. Rouillier, M.-F. Roy and M. Safey el Din [51] 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. This is done by computing the (complex) zero-set of the system

{P (M) = 0, grad_{M}P//AM},

where the last condition is expressed by setting (n − 1) 2 × 2 determinants to zero.

We suppose we have a generic enough point A, so this system is zero-dimensional if the hypersurface H has a finite number of singularities. When H has an infinite number of singularities, the solution proposed in [51] consists in introducing an infinitesimal ε and studying the critical points on the level sets P^{−1}(ε). This amounts to solving the system

{P (M) = ε, grad_{M}P//AM}.

We see this system as parametrized by ε; Sard’s theorem shows that all solutions of this system are generically simple, so our algorithm can be used to compute a parametric solution. We mention that the final step, described in [51], amounts to studying to limits when ε → 0 of the points described by the parametric resolution, so as to solve the initial problem.

As an illustration, in [52], with F. Rouillier and M. Safey el Din, we treated some examples taken from the Birkhoff Interpolation Problem [31], involving some hypersurfaces in 3 variables (t_{2}, t_{3}, t_{4}), which had not be treated automatically before. In Section 6.2, we consider the two examples

P_{3} = 3t_{3}^{6} + 3t_{4}^{6} + 9t_{3}^{2}t_{4}^{4} + 9t_{3}^{4}t_{4}^{2} + t_{2}^{2}t_{4}^{4} + t_{2}^{2}t_{3}^{4} − t_{2}^{4}t_{4}^{2} − t_{2}^{4}t_{3}^{2}+
t_{2}^{6} + 4t_{4}^{4} + 4t_{3}^{4} + 4t_{2}^{4} + 2t_{2}^{2}t_{3}^{2}t_{4}^{2} + 8t_{3}^{2}t_{4}^{2} − 4t_{2}^{2}t_{4}^{2} − 4t_{2}^{2}t_{3}^{2},
P_{10} = −t_{3}^{6} − 3 − 4t_{3}^{4}t_{4}^{2} + t_{2}^{2}t_{3}^{4} − 4t_{2}^{4}t_{4}^{2} − t_{2}^{4}t_{3}^{2} − 3t_{2}^{6} + t_{3}^{4} − 9t_{2}^{4}−
t_{3}^{2} − 9t_{2}^{2} − 4t_{4}^{2} + 4t_{2}^{2}t_{3}^{2}t_{4}^{2} − 2t_{2}^{2}t_{3}^{2} − 8t_{2}^{2}t_{4}^{2} + 4t_{3}^{2}t_{4}^{2}.

Using a randomly chosen point A with integer coordinates, we generate two systems called Birkhoff_{3} and Birkhoff_{10}.

Computing relations. This system was solved to answer a question of I. Bershenko-Kogan [7]. The initial goal is the determination of the conjugacy classes of ℚ[A, B, C] under the action of GL_{3}(ℚ), and notably the study of the orbit of the form A^{n} + B^{n} + C^{n}.

This requires to compute the relation R between three rational functions X/D, Y/D, Z/D in two variables P, Q. This relation is an equation of the image of the corresponding rational application, so it is an irreducible polynomial.

We viewed the system the following way: we introduce two parameters x and y, three variables z, P, Q, and the system (f_{1}, f_{2}, f_{3}) = (Dx − X, Dy − Y, Dz − Z). The relation R is the minimal polynomial of the variable z in ℚ(x, y) → ℚ(x, y)[P, Q, z]/(f_{1}, f_{2}, f_{3}). Once we have a parametric resolution, computing this minimal polynomial is easy: a straightforward way to do so is the computation of the squarefree part of a resultant; a more efficient solution is described in [61].

This system is denoted Bershenko is the sequel. Our output was checked in a direct manner: since the relation R must be irreducible, it suffices to evaluate the relation R on the functions X/D, Y/D, Z/D, and check that we obtain zero.

The system Hawes. This last example is taken from the database SymbolicData [5] (see also J.-C. Faugère’s homepage [1]). Its purpose is to illustrate the behavior of our algorithm with respect to the representation of the input system.

This system has 6 equations plus two inequations in the 8 variables a, b, c, x, y_{1}, z_{1}, y_{2}, z_{2}.

3z_{1}^{2} + y_{1}^{2} + b = 0,
3z_{2}^{2} + y_{2}^{2} + b = 0,
2cy_{1} + x + 5y_{1}^{4} + 3ay_{1}^{2} + 2y_{1}z_{1} = 0,
2cy_{2} + x + 5y_{2}^{4} + 3ay_{2}^{2} + 2y_{2}z_{2} = 0,
−cy_{2}^{2} + cy_{1}^{2} − xy_{2} + xy_{1} + z_{1}^{3} + y_{1}^{2}z_{1} + bz_{1} − y_{2}^{5} − ay_{2}^{3} − y_{2}^{2}z_{2} − z_{2}^{3} − bz_{2} + y_{1}^{5} + ay_{1}^{3} = 0,
2(30y_{1}^{3}z_{1} − y_{1}^{2} + 9ay_{1}z_{1} + 3z_{1}^{2})(−9ay_{2}^{2}z_{2} + 3xz_{2} − 3y_{2}z_{2}^{2} − 45y_{2}^{4}z_{2} + y_{2}^{3} − by_{2})
−2(30y_{2}^{3}z_{2} − y_{2}^{2} + 9ay_{2}z_{2} + 3z_{2}^{2})(−9ay_{1}^{2}z_{1} + 3xz_{1} − 3y_{1}z_{1}^{2} − 45y_{1}^{4}z_{1} + y_{1}^{3} − by_{1})
+6cy_{1}^{2}z_{1}(3xz_{2} − 3y_{2}z_{2}^{2} − 45y_{2}^{4}z_{2} + y_{2}^{3} − by_{2}) − 6cy_{2}^{2}z_{2}(3xz_{1} − 3y_{1}z_{1}^{2} − 45y_{1}^{4}z_{1} + y_{1}^{3} − by_{1}) = 0.

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}”. As in the previous example, this amounts to computing the minimal polynomial with coefficients in ℚ(a, b) of c in a suitably-defined quotient algebra. Our solution was described in the previous paragraph; here, we only consider the preliminary task, the computation of a parametric resolution.

The system was initially given under an expanded form, requiring 92 multiplications. An additional work was necessary to recover the original, more compact, formulation given above, with 55 operations. The impact of this reformulation was important on the computation times, as the following tables will show.

6.2 Computation times

We now present the computation times. As the complexity results are stated in terms of arithmetic operations, we first give the times for the resolution over a finite fields, where arithmetic operations have a constant cost; we choose the field 𝔽_{𝔭}, 𝔭 = 10000000000000000051 being the first prime greater than 10^{19}. In the second table, we present the results of the computations over ℚ, for which we used the strategy described in Section 5.3.

The computations were done on the machines of the UMS MEDICIS [3], on Compaq Alpha EV6 XP/1000 500 Mhz processors with 640 MB of RAM, using Magma v. 2.6.

Here is the legend for Figures 2 and 3:

• The first four lines give the measure of the input: (n, m, d) and the number of multiplications in the Straight-Line Program giving the system;

• The next two lines give the generic number of solutions (deg_{π}) and the degree in the parameters of the coefficients of the output (deg_{u}).

Figure 2: Computation over 𝔽_{10000000000000000051}

System  P19 Birkhoff_{3} Birkhoff_{10} Bershenko  Hawes
Variables  4  3  3  3  6
Parameters  2  1  1  2  2
Degree  10  6  6  12  9
Multiplications  253  91  80  32  55 (92)
Generic Degree  64  34  136  24  120
Degree of the Coefficients 44  5  19  14  46
Probability of Error  0.01 % 2. 10^{−6}% 2. 10^{−6}%  3 %  0.26 %
Time  2.6 h.  6 s.  8 m.  1 min. 6 h. (14 h.)
Magma (Gröbner)  ∞  11 s.  > 5 h.  10 min.  ∞

Figure 3: Computation over ℚ

System  Birkhoff_{3} Birkhoff_{10} Bershenko  Hawes
Time  15 s.  20 min.  7 min. 15 h. (35 h.)
Magma (Gröbner) 5 min.  ∞  ∞  ∞

• The next line gives the probability of error, if the 3m−1 values are chosen uniformly in 𝔽_{10000000000000000051}. We use the precise evaluation of the probability given in Section 5.

• The last lines in Figure 2 compare the running times of our algorithm with the time taken by Magma’s Gröbner engine, on the field k(P_{1}, …, P_{m}), for a lexicographic ordering. We have used the factorization method described in Section 5.3.

• The two lines in Figure 3 give the times of the computation over ℚ.

• In the last column, the figures in parentheses indicate the complexity of evaluation of the developed form of the system, and the corresponding computation times.

Our algorithm behaves very well, and outperforms the Gröbner basis computation proposed by Magma on all these examples. A more extensive list of examples is given in the paper [52] and the author’s PhD. Thesis [58], and confirms its good behavior. For the biggest examples, the size of the output becomes the limiting factor: writing down the output of the system Hawes takes more than 20 MB.

The counterpart is a possibility of failure, but the estimates are quite reasonable: there is never more than a few percents of chance that the algorithm fails. When a verification was possible, it never revealed an error.

Finally, the influence of the complexity of evaluation L is predominant for the running time, as the example Hawes clearly shows: (almost) doubling the number of multiplications has the immediate effect of (almost) doubling the computation time. We stress that, as to no surprise, most applications we have met can be formulated in an easy-to-evaluate form, which is to the advantage of our method. Yet, we have observed that the modelization through a Computer Algebra system may spoil such good behavior, since most systems represent polynomials through the list of their coefficients on a monomial basis.

7 Conclusion

In this article, we have studied the geometry of parametric systems, proposed an elimination procedure adapted to such situations, and demonstrated its good practical behavior. Here are some directions for future work.

• As mentioned earlier, when the base field is ℚ, the first task is to take modular computations into account, and in particular to quantify the new possibilities of degeneracy. This requires to use the arithmetic forms of the geometric results used here; the kind of results we need is given for instance in [38]; the subsequent algorithm is given in the author’s PhD. Thesis [58].

• Our algorithm works only for the components were the Jacobian determinant is generically invertible. Recently, G. Lecerf proposed in [41] an extension of the Newton lifting process for multiple components, whose projection on the parameter space is dominant and generically finite. This new tool will enable a extension of our algorithm to the general case.

• Finally, we do not yet make use of the full strength of our Newton operator, which applies for more general representations than those based on primitive elements. A natural generalization of our algorithm is to use an encoding of the output by triangular sets. We refer to [59], where these aspects are developed.

References

[1] Jean-Charles Faugère’s homepage. http://posso.lip6.fr/ jcf/.

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

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

[4] Alonso, M. E., Becker, E., Roy, M.-F., and Wörmann, T. Zeros, multiplicities and idempotents for zero-dimensional systems. In Proceedings of MEGA’94 (1996), vol. 142 of Progress in Mathematics, Birkhäuser, pp. 1–15.

[5] Bachmann, O., Dengel, M., Faugère, J.-C., Graebe, H.-G., and Schoenemann, H. The SymbolicData Project. http://symbolicdata.medicis.polytechnique.fr/.

[6] Baur, W., and Strassen, V. The complexity of partial derivatives. Theoretical Computer Science 22 (1983), 317–330.

[7] Bershenko-Kogan, I. Inductive approach to Cartan’s Moving Frame Method with applications to classical invariant theory. PhD thesis, University of Minesotta, 2000.

[8] Boulier, F., Lazard, D., Ollivier, F., and Petitot, M. Representation for the radical of a finitely generated differential ideal. In Proceedings of ISSAC’95 (1995), ACM Press, pp. 158–166.

[9] Buchberger, B. Gröbner bases: An algorithmic method in polynomial ideal theory. In Multidimensional System Theory. Reidel, Dordrecht, 1985, pp. 374–383.

[10] Buchberger, B., Collins, G., and Loos, R., Eds. Computer Algebra. Springer Verlag, 1982.

[11] Bürgisser, P., Clausen, M., and Shokrollahi, M. A. Algebraic Complexity Theory. Springer, 1997.

[12] Canny, J. Some algebraic and geometric problems in PSPACE. In Proceedings of the 20th ACM Symposium on Theory of Computation (1988), pp. 460–467.

[13] Cohen, H. Advanced Topics in computational number theory. No. 193 in Graduate Texts in Mathematics. Springer-Verlag, 2000.

[14] Dellière, S. Triangularisation de systèmes constructibles — Application à l’évaluation dynamique. PhD thesis, Université de Limoges, 1999.

[15] Dixon, J. Exact solution of linear equations using p-adic expansions. Numerische Mathematik 40 (1982), 137–141.

[16] Duval, D. Diverses questions relatives au calcul formel avec des nombres algébriques. PhD thesis, Université scientifique, technologique et médicale de Grenoble, 1987.

[17] Eisenbud, D. Commutative Algebra with a view toward Algebraic Geometry. No. 150 in Graduate Texts in Mathematics. Springer, 1996.

[18] Emiris, I., and Mourrain, B. Matrices in elimination theory. Journal of Symbolic Computation 28, 1–2 (1999).

[19] Emiris, I., and Rege, A. Monomial bases and polynomial system solving. In Proceedings of ISSAC’94 (1994), ACM Press, pp. 114–122.

[20] Faugère, J.-C. A new efficient algorithm for computing Gröbner bases (F_{4}). Journal of Pure and Applied Algebra 139, 1–3 (1999), 61–88.

[21] Faugère, J.-C., Gianni, P., Lazard, D., and Mora, T. Efficient computation of zero-dimensional Gröbner bases by change of ordering. Journal of Symbolic Computation 16, 4 (1993), 329–344.

[22] Fitzpatrick, P., and Flynn, J. A Gröbner basis technique for Padé approximation. Journal of Symbolic Computation 13 (1992), 133–138.

[23] Gaudry, P., and Harley, R. Counting points on hyperelliptic curves over finite fields. In Proceedings of ANTS IV (2000), no. 1838 in Lecture Notes in Computer Science, pp. 313–332.

[24] Gianni, P., and Mora, T. Algebraic solution of systems of polynomial equations using Gröbner bases. In Proceedings of AAECC-5 (1989), vol. 356 of Lecture Notes in Computer Science, Springer, pp. 247–257.

[25] Giusti, M., Hägele, K., Heintz, J., Morais, J.-E., Montaña, J.-L., and Pardo, L.-M. Lower bounds for Diophantine approximation. In Proceedings of MEGA’96 (1997), no. 117, 118 in Journal of Pure and Applied Algebra, pp. 277–317.

[26] Giusti, M., Heintz, J., Morais, J.-E., Morgenstern, J., and Pardo, L.-M. Straight-line programs in geometric elimination theory. Journal of Pure and Applied Algebra 124 (1998), 101–146.

[27] Giusti, M., Heintz, J., Morais, J.-E., and Pardo, L.-M. When polynomial equation systems can be solved fast? In Proceedings of AAECC-11 (1995), vol. 948 of Lecture Notes in Computer Science, Springer, pp. 205–231.

[28] Giusti, M., Lecerf, G., and Salvy, B. A Gröbner free alternative for polynomial system solving. Journal of Complexity 17, 1 (2001), 154–211.

[29] Gomez-Diaz, T. Quelques applications de l’évaluation dynamique. PhD thesis, Université de Limoges, 1994.

[30] Gonzalez-Vega, L. Computer Algebra in Science and Engineering. World Scientific Publishing, 1995, ch. Some examples of problem solving by using the symbolic viewpoint when dealing with polynomial systems of equations, pp. 102–116.

[31] Gonzalez-Vega, L. Applying quantifier elimination to the Birkhoff interpolation problem. Journal of Symbolic Computation 22 (1996), 83–104.

[32] Grigoriev, D., and Vorobjov, N. Bounds on numbers of vectors of multiplicities for polynomials which are easy to compute. In Proceedings of ISSAC 2000 (2000), ACM Press, pp. 137–145.

[33] Heintz, J. Definability and fast quantifier elimination in algebraically closed fields. Theoretical Computer Science 24, 3 (1983), 239–277.

[34] Heintz, J., Krick, T., Puddu, S., Sabia, J., and Waissbein, A. Deformation techniques for efficient polynomial equation solving. Journal of Complexity 16, 1 (2000), 70–109.

[35] Heintz, J., Matera, G., and Waissbein, A. On the time-space complexity of geometric elimination procedures. Applicable Algebra in Engineering, Communication and Computing 11 (2001), 239–296.

[36] Heintz, J., Roy, M.-F., and Solernó, P. On the complexity of semi-algebraic sets. In Proc. IFIP 89, San Francisco (1989), North-Holland, pp. 293–298.

[37] Heintz, J., Roy, M.-F., and Solernó, P. On the theoretical and practical complexity of the existential theory of the reals. Comput. J. 36, 5 (1993), 427–431.

[38] Krick, T., Pardo, L. M., and Sombra, M. Sharp estimates for the arithmetic Nullstellensatz. In Proceedings of EACA’99 (1999).

[39] Kronecker, L. Grundzüge einer arithmetischen Theorie der algebraischen Grössen. J. reine angew. Math. 92 (1882), 1–122.

[40] Lecerf, G. Kronecker, a package for Magma for polynomial system solving. UMS MEDICIS, Laboratoire GAGE, http://kronecker.medicis.polytechnique.fr/, 1999.

[41] Lecerf, G. Quadratic Newton iteration for systems with multiplicity. Manuscript, Laboratoire GAGE, École polytechnique, France, April 2001.

[42] Lecerf, G., and Schost, É.. Fast multivariate power series multiplication in characteristic zero. Manuscript, Laboratoire GAGE, École polytechnique, France, April 2001.

[43] Leverrier, U. J. J. Sur les variations séculaires des éléments elliptiques des sept planètes principales: Mercure, Venus, la terre, Mars, Jupiter, Saturne et Uranus. Journal de Mathématiques Pures et Appliquées 4 (1840), 220–254.

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

[45] Montes, A. A new algorithm for discussing Gröbner bases with parameters. Journal of Symbolic Computation 33, 2 (2002), 183–208.

[46] Morrison, S. The differential ideal [P] : M^{∞}. Journal of Symbolic Computation 28 (1999), 631–656.

[47] Mourrain, B., and Pan, V. Multivariate polynomials, duality and strucured matrices. Journal of Complexity 16, 1 (2000).

[48] 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), 255–299.

[49] Rojas, J. Computing complex dimension faster and deterministically. Tech. rep., Department of Mathematics, City University of Hong Kong, 2000.

[50] Rouillier, F. Solving zero-dimensional systems through the Rational Univariate Representation. Applicable Algebra in Engineering, Communication and Computing 9, 5 (1999), 433–461.

[51] Rouillier, F., Roy, M., and Safey el Din, M. Finding at least one point in each connected component of a real algebraic set defined by a single equation. Journal of Complexity 16, 4 (2000), 716–750.

[52] Rouillier, F., Safey el Din, M., and Schost, É. Solving the Birkhoff interpolation problem via the critical point method : an experimental study. In Proceedings of ADG 2000 (2001), vol. 2061 of Lecture Notes in Artificial Intelligence, Springer-Verlag.

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

[54] Sakata, S. Extension of the Berlekamp-Massey algorithm to n dimensions. Information and Computation 84 (1990), 207–239.

[55] Schönhage, A. Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2. Acta Informatica, 7 (1977), 395–398.

[56] Schönhage, A., and Strassen, V. Schnelle Multiplikation großer Zahlen. Computing 7 (1971), 281–292.

[57] Schoof, R. Elliptic curves over finite fields and the computation of square roots mod p. Mathematics of Computation 44 (1985), 483–494.

[58] Schost, É.. Sur la résolution des systèmes polynomiaux à paramètres. PhD thesis, École polytechnique, 2000.

[59] Schost, É. Degree bounds and lifting techniques for triangular sets. In Proceedings of ISSAC 2002 (2002).

[60] Schwartz, J. T. Fast probabilistic algorithms for verification of polynomial identities. Journal of the ACM 27, 4 (1980), 701–717.

[61] Shoup, V. Efficient computation of minimal polynomials in algebraic extensions of finite fields. In Proceedings of ISSAC’99 (1999), ACM Press, pp. 53–58.

[62] Thue, A. Et par andtydninger til en talteoretisk methode. Christiana 7 (1902).

[63] Trinks, W. On improving approximate results of Buchberger’s algorithm by Newton’s method. In Proceedings of EUROCAL’85 (1985), no. 204 in Lecture Notes in Computer Science, Springer-Verlag, pp. 608–611.

[64] von zur Gathen, J., and Gerhard, J. Modern computer algebra. Cambridge University Press, 1999.

[65] Weispfenning, V. Comprehensive Gröbner bases. Journal of Symbolic Computation 14, 1 (1992), 1–30.

[66] Winkler, F. A p-adic approach to the computation of Gröbner bases. Journal of Symbolic Computation 6 (1988), 287–304.

[67] Zippel, R. Probabilistic algorithms for sparse polynomials. In Proceedings of EU-ROSAM’79 (1979), no. 72 in Lecture Notes in Computer Science, Springer, pp. 216– 226.