On sparse interpolation of rational functions and gcds Joris van der Hoevena, Grégoire Lecerfb CNRS, École polytechnique, Institut Polytechnique de Paris Laboratoire d'informatique de l'École polytechnique (LIX, UMR 7161) 1, rue Honoré d'Estienne d'Orves Bâtiment Alan Turing, CS35003 91120 Palaiseau, France a. Email: vdhoeven@lix.polytechnique.fr b. Email: lecerf@lix.polytechnique.fr Preliminary version of November 4, 2023

1.Introduction

Let be a rational function over an effective field , presented as a blackbox that can be evaluated at points in . The problem of sparse interpolation is to recover in its usual representation, as a quotient of the form

where and for , where and for , where , and where for any .

Here, “sparse” means that we use a sparse representation for multivariate polynomials: only the set of the non-zero terms is stored and each term is represented by a pair made of a coefficient and an exponent vector. When using a sparse representation, the natural complexity measures are the number of non-zero terms and the maximum bit size of the exponents. On the contrary, a dense representation of a polynomial , say of total degree , is the vector of all its coefficients up to degree for a prescribed monomial ordering.

The sparse interpolation of a blackbox polynomial has been widely studied in computer algebra, and very efficient algorithms are known. In this paper, we will use this polynomial case as a basic building block, with an abstract specification. We refer to [1, 4, 8, 9, 13, 17] for state of the art algorithms for polynomial sparse interpolation and further historical references.

Our main problem here is to find an efficient reduction of the general sparse interpolation problem for rational functions to the special case of polynomials. We focus on the case when the total degrees and are “modest” with respect to the total number of terms . This is typically the case for applications in computer algebra. We also recall that algorithms in this area are usually probabilistic, of Monte Carlo type; our new algorithms are no exception. In a similar vein, we assume that the cardinality of is “sufficiently large”, so that random points in are “generic” with high probability.

One efficient reduction from sparse rational function interpolation to sparse polynomial interpolation was proposed by Cuyt and Lee [2]. Our new method is a variant of their algorithm; through the use of projective coordinates, we are able to prove a better complexity bound: see section 3.

Computations with rational functions are intimately related to gcd computations. In section 4, we show that our approach also applies to gcd computations: given a blackbox function for the computation of two sparse polynomials and , we present a new complexity bound for the sparse interpolation of . We also present an alternative heuristic algorithm for the case when and are given explicitly, in sparse representation.

2.Preliminaries

For simplicity, all time complexities are measured in terms of the required number of operations in . Since our algorithms are probabilistic, we assume that a random element in can be produced in unit time. Our complexity bounds are valid in several complexity models. For definiteness, one may assume a probabilistic RAM model that supports operations in .

2.1.On the complexity of sparse polynomial interpolation

The algorithms in this paper all rely on a given base algorithm for the interpolation of sparse polynomials. For a polynomial in variables of total degree with terms, given by a blackbox function , we assume that this algorithm requires evaluations of , as well as additional operations for the actual interpolation from the latter values of . We recall that this algorithm is allowed to be probabilistic, of Monte Carlo type.

As a particular example of interest, consider the case when , and let be the cost to multiply two polynomials of degree in . In a suitable asymptotic region where and do not grow too fast, Prony's traditional probabilistic geometric progression technique yields

For FFT-primes , one may also apply the tangent Graeffe method, which yields

We refer to [9] for a survey of recent complexity bounds for sparse polynomial interpolation. In what follows, we always assume that and are non-decreasing in .

2.2.Genericity and random dilatations

Consider a polynomial or rational blackbox function . We assume that a blackbox function over only involves field operations in (which includes equality testing and branching, if needed), and that the number of these operations is uniformly bounded.

Efficient algorithms for the sparse interpolation of typically evaluate at a specific sequence of points, such as a geometric progression. For the algorithms in this paper, it is convenient to systematically make the assumption that we only evaluate at points with for . This assumption holds in particular for points in geometric progressions.

One problem with the evaluation of a blackbox function that involves divisions is that the evaluation may fail whenever a division by zero occurs. An evaluation point for which this happens is colloquially called a bad point. Such bad points must be discarded and the same holds for sequences of evaluation points that contain at least one bad point. Of course, it would in principle suffice to choose another sequence of evaluation points. However, in general, sufficiently generic sequences are not readily at our disposal with high probability. This problem may be overcome using randomization, as follows.

Since our blackbox function only involves a bounded number of field operations in , the set of bad evaluation points is a constructible subset of . We assume that this subset is not Zariski dense. Consequently, there exists a polynomial such that whenever is a bad point. Given , the probability that a point satisfies

is at least

by the Schwartz–Zippel lemma [5, Lemma 6.44], where denotes the cardinality of .

Now assume that we are given a method for the sparse interpolation of any blackbox function with the same support as . The evaluation points for such a method form a fixed finite subset of . If is sufficiently large, then it follows with high probability that for all , when taking at random. We may thus interpolate

instead of , with high probability. Since are non-zero, we may next recover the sparse interpolation of itself.

3.Rational functions

This section concerns the sparse interpolation of a rational function given as a blackbox.

3.1.Related work

The natural idea to interpolate a univariate rational function given as a blackbox is to evaluate it at many points and then use Cauchy interpolation [5, Corollary 5.18] in combination with the fast Euclidean algorithm [5, chapter 11]; this idea appears in early works on the topic [14]. If the numerator and denominator of have respective degrees and , then the total cost decomposes into evaluations of and

operations in for the Cauchy interpolation. Here is the usual shorthand for , as in [5].

Interpolating rational functions seems more difficult than polynomials. In high degree, it is still a challenging problem to design efficient algorithms in terms of the sparse size of the function. An algorithm with exponential time complexity is known from [6].

In the multivariate case, the interpolation problem can be handled as a dense one with respect to one of the variables and as a sparse one with respect to the other variables. Although no general polynomial complexity bound is known in terms of the sparse size, this approach is extremely useful in practice, as soon as one of the variables has small degree. Roughly speaking, one proceeds as follows: for a fixed point we reconstruct the univariate function with the dense representation in . Then we regard the numerator and the denominator as polynomials in , and it remains to interpolate their sparse coefficients using sufficiently many evaluation points .

The difficulty with the latter interpolation is that the values of and are only determined up to a scalar in . In order to compute this scalar, a normalization assumption is necessary. For instance we may require to be monic in , which can be achieved through a random linear change of coordinates [14]. Another kind of normalization has been proposed by Cuyt and Lee [2]: they introduce a new auxiliary variable to decompose and into their homogeneous components. Their approach is sketched below. Let us finally mention [12], which addresses the special case using a specific strategy.

3.2.The algorithm by Cuyt and Lee

An efficient reduction of sparse interpolation from rational functions to polynomials was proposed by Cuyt and Lee in [2]. Their approach starts by introducing a homogenizing variable and

With high probability, they next reduce to the case when is non-zero, by choosing a random point and performing a shift:

For a fixed point , we may consider as a univariate rational function in , which can generically be interpolated from evaluations for different values of .

The first main idea from [2] is to normalize the resulting rational function in such a way that the constant coefficient of the denominator is equal to one:

This yields blackbox functions that take as input and produce

as output. We may now use sparse polynomial interpolation for each of these blackbox functions and determine the sparse expressions for the above polynomials. We finally obtain using a simple shift.

One complication here is that shifts

destroy the sparsity. The second main observation by Cuyt and Lee is that this drawback only occurs for the non-leading terms: the coefficient (resp. ) coincides with (resp. ). Using this observation recursively in a clever way, it is possible to obtain the non-leading coefficients as well, using at most

blackbox evaluations [2, section 2.2]. However, even though sparse interpolations of shifted polynomials can be avoided, the resulting algorithm still requires evaluations of shifted polynomials in their interpolated form: see the remarks at the end of [2, section 2.2]. In fact, the focus of [2] is on minimizing the number of blackbox evaluations; the overall complexity of the algorithm is not analyzed in detail.

3.3.A variant based on projective coordinates

Instead of using the recursive method from [2, section 2.2] for determining the non-leading coefficients, we propose to avoid the mere existence of such terms altogether. At the start of our algorithm, we assume that the total degrees and of the numerator and the denominator of are known. If we are only given a rough upper bound on these degrees, then we may reduce to this case as follows: take a point at random, and interpolate as a univariate rational function in , from evaluations at different values of . With high probability, and coincide with the degrees of the numerator and denominator of .

Our key observation is the following: if and are homogeneous polynomials then and , so Cuyt and Lee's method mostly reduces to two sparse polynomial interpolations. Our method is therefore to force the polynomials and to become homogeneous by introducing a new projective coordinate .

More precisely, we introduce the rational function

and we write

where

The numerator and denominator of are both homogeneous and coprime. One evaluation of requires one evaluation of plus operations in .

For some random shift , we next consider

One evaluation of requires one evaluation of plus operations in . Before addressing the interpolation algorithm for , let us show that is the canonical representative of .

Proof. We consider the following weighted degree:

For this weight, and are quasi-homogeneous, whence so are their respective irreducible factors and therefore their gcd, written . Consequently, there exists a non-negative integer such that

By construction, and are coprime, whence is a non-zero constant and . Since , we conclude that .

Theorem 2. Let be a rational blackbox function whose evaluation requires operations in . Let , where and are bounds for the number of terms of and . Let , where and are the total degrees of and . Using a probabilistic algorithm of Monte Carlo type, the sparse interpolation of takes at most

operations in , provided that the cardinality of is sufficiently large.

Proof. First of all, a point is taken at random, so holds with high probability. Let us describe a blackbox for evaluating and at a point :

• We perform evaluations of for different values of , which give rise to at most

  

  operations in .

• If the latter evaluations do not fail, then we perform the univariate Cauchy interpolation of , which takes operations in . If this Cauchy interpolation fails, then the evaluation of and aborts.

The failure of the Cauchy interpolation of depends on the actual procedure being used, but we may design it to succeed on a Zariski dense subset of points in . For instance, with the aforementioned method of [5, Corollary 5.18], the reconstruction succeeds whenever the degrees of the reconstructed numerator and denominator are and , respectively. This holds when , where

Since is primitive as an element of , the gcd of and in can be obtained as the primitive representative in of the gcd of and in ; see [15, chapter IV, Theorem 2.3], for instance. Using Lemma 1, we deduce that is not identically zero. Consequently, the evaluation of and at a point does not fail on a Zariski dense subset of .

In order to ensure the sparse polynomial interpolation to succeed with high probability, we appeal to the discussion from section 2.2 for the evaluations of and . More precisely, by taking at random, we can interpolate and with high probability. This random dilatation further increases the cost of each blackbox evaluation by operations in .

The required number of evaluation points is

Once all evaluations are done, we need

further operations to reconstruct the sparse representations of and . The backward dilatation to recover and takes

further operations in .

To conclude the proof, we observe that and coincide with and . But since are are homogeneous, this means that and actually coincide with and .

Remark 3. One nice aspect of Cuyt and Lee's algorithm is that the number of blackbox evaluations can be of the order of

if the terms of and are more or less equally distributed over their homogeneous components. With our method, we always pay the maximal possible overhead with respect to sparse interpolation of polynomials.

4.Sparse gcds

Let be polynomials that are given through a common blackbox. In this section we show how to adapt the technique from the previous section to the computation of .

4.1.Related work

A classical approach for multivariate gcds relies on evaluation/interpolation for all but one of the variables, say , and thus on the computation of many univariate gcds; see for instance [5, chapter 6]. This turns out to be similar to rational function interpolation: the univariate case is handled using dense representations and a suitable normalization is necessary to recombine the specialized gcds.

In general and without generic coordinates, the gcd does not admit a natural normalization. This issue has been studied in [3]. Basically one piece of their solution is to consider and as polynomials in and to force to be monic in . This leads to interpolate rational functions instead of polynomials, and the technique of the previous section may be used. But it is even more efficient to adapt our main idea for rational functions to gcds.

The sparse gcd problem was first investigated by Zippel, who proposed a variable by variable approach [18]. Then, subresultant polynomials turned out to be useful to avoid rational function interpolations. For historical references we refer the reader to [3, 14]. Software implementations and a recent state of the art can be found in [11].

Let us mention that probabilistic sparse gcd algorithms may be turned into Las Vegas type algorithms by interpolating the full Bézout relation, at the price of an increased computational cost; see for instance [5, chapter 6]. In this paper we content ourselves with Monte Carlo algorithms and we focus on gcd computations only.

4.2.A variant based on projective coordinates

Let be polynomials that are given through a common blackbox, and assume that contains terms. We let , , . Using the technique from the previous section we first reduce to the case when and are homogeneous (without increasing their total degrees), so we let:

For , we let

Example 4. With , , , , we have and . We note that , and then that . However .

The latter example shows that we need to take care of the relationship between and . This is the purpose of the following lemma, which extends Lemma 1.

Proof. Let . In the proof of Lemma 1 we have shown that

for some . By construction, we have , so divides . Conversely, divides , so there exists a non-zero constant with

It follows that

Since or , we finally have , so that .

From now on we assume that or , so that or is primitive as an element in . By Lemma 5 and [15, chapter IV, Theorem 2.3], their gcd in can be obtained as the primitive representative in of the gcd of and in . We further observe that , whence

(1)

This allows us to obtain

as the gcd of and in , and to recover via

This approach yields the following complexity bound for the sparse interpolation of .

Theorem 6. Consider a blackbox function that computes the two polynomials and in using operations in . Let , , and be the respective total degrees of , , and . Let . Using a probabilistic algorithm of Monte Carlo type, we may interpolate using

operations in , provided that the cardinality of is sufficiently large.

Proof. We first take at random, so or with high probability. In order to evaluate at a point we proceed as follows:

• We evaluate and for different values of . The evaluation of , , takes

  

  operations in . The evaluation of involves more operations. Therefore, one evaluation of and amounts to

  

  operations in .

• We interpolate and , and then compute

  

  This requires operations in . By (1), the normalization coincides with whenever , , and the subresultant coefficient of degree of and is non-zero: see [5, chapter 6] or [16], for instance.

In this way, we have built a blackbox that evaluates on a Zariski dense subset of . In order to ensure that the sparse interpolation of succeeds with high probability, we appeal to the discussion from section 2.2. Precisely, we take at random and we interpolate . This random dilatation further increases the cost of each blackbox evaluation by operations in .

We need evaluation points, together with further operations in , for the sparse interpolation of . Finally, we recover the requested gcd of and from , using extra operations in .

4.3.Another heuristic approach

Let us now focus on the case when and are presented as sparse polynomials instead of blackbox functions (more generally, one may consider the model from [7, section 5.3]). The number of terms of (resp. of ) is written (resp. ). The above method has the disadvantage that we need to evaluate and at many shifted points that are generally not (e.g.) in a geometric progression. Such evaluations are typically more expensive, when using a general algorithm for multi-point evaluation, so it is better to avoid shifting altogether in this case.

Now consider the particular case when

and contains a single term . Then we may use an alternative normalization where we divide by instead of , and multiply by to keep things polynomial. For a random point , we may check whether we are in this case (with high probability) by testing whether

Once has been interpolated, we recover as

where stands for the gcd of all monomials occurring in and similarly for and . With these modifications, the complexity bound of Theorem 6 becomes

The same approach can be used if contains a single term .

In the case when contains several terms, we consider

for a suitable vector of weights. Whenever contains a single term, we may use the above method to compute and then . However, the degree of is potentially times larger than the degree of , where . This leads to the question how small we can take while ensuring that contains a single term with probability at least . We conjecture that this is the case when taking at random with . In practice, the idea is to try a sequence of random weights for which follows a slow geometric progression.

Let us provide some partial evidence for our conjecture. In dimension , consider the Newton polygon at of a polynomial of degree made of the maximal number of slopes with . For each integer , let be the number of slopes with , , and . We have

Under this constraint, the sum is maximized by considering Newton polygons that first exhaust all slopes for which is minimal. Now for all . Taking minimal with , it follows that

Now , so and

On the other hand, for any , the number of rational numbers with is at least

Now let be the smallest integer with . For a random weight with , we then conclude that has a unique term with probability at least .

As a second piece of evidence, let us show that there always exists a weight vector for which contains a unique term. Let be an element of the support of for which is maximal. Then is included in the ball with center and radius , so the tangent space to this ball at intersects only at . Taking to be the normal vector to this tangent space, it follows that for some .

Bibliography