Generic bivariate multi-point evaluation, interpolation and modular composition with precomputation
Vincent Neiger Johan Rosenkilde Grigory Solomatov
cs.SC






  Univ. Limoges, CNRS, XLIM, UMR 7252
  F-87000 Limoges
  France





  Technical University of Denmark
  Kgs. Lyngby
  Denmark





  Technical University of Denmark
  Kgs. Lyngby
  Denmark




  Suppose  is a large enough field and ⊂^2 is a fixed, generic set of points which is available for precomputation.
  We introduce a technique called reshaping which allows us to design quasi-linear algorithms for both:
  computing the evaluations of an input polynomial f ∈[x,y] at all points of ;
  and computing an interpolant f ∈[x,y] which takes prescribed values on  and satisfies an input y-degree bound.
  Our genericity assumption is explicit and we prove that it holds for most point sets over a large enough field.
  If  violates the assumption, our algorithms still work and the performance degrades smoothly according to a distance from being generic.
  To show that the reshaping technique may have an impact on other related problems, we apply it to modular composition: suppose generic polynomials M ∈[x] and A ∈[x] are available for precomputation, then given an input f ∈[x,y] we show how to compute f(x, A(x))  M(x) in quasi-linear time.








































Generic Bivariate Multi-point Evaluation, Interpolation and Modular Composition with Precomputation
    Grigory Solomatov
    March 2020
===================================================================================================





§ INTRODUCTION





  
Outline
Let  be an effective field.
We consider the three classical problems for bivariate polynomials [x,y] mentioned in the title.
We assume a model where part of the input is given early as preinput which is available for heavier computation, and the primary goal is to keep the complexity of the online phase, once the remaining part of the input is given, to a minimum.

Multi-point evaluation (MPE): with preinput a point set = {(α_i,β_i)}_i=1^n ⊆^2 and input f ∈[x,y], compute (f(α_i,β_i))_i=1^n.
We give two algorithms: the first requires pairwise distinct α_i's and has online complexity f  f + n as long as  is balanced, a notion described below;
  the second accepts repeated x-coordinates with online complexity f( f +  f) + n as long as a certain “shearing” of  is balanced.
“soft-O” ignores logarithmic terms: f(n)(log f(n))^c⊂f(n) for any c ∈_≥ 0.

Interpolation: with preinput a point set  as before, and input values γ⃗∈^n, compute f ∈[x,y] such that (f(α_i,β_i))_i=1^n = γ⃗, satisfying some constraints on the monomial support.
We give an algorithm which preinputs a degree bound d and outputs f such that f < d and f ∈n/d.
The online complexity is n if  and a shearing of  are both balanced; d should exceed the x-valency of , i.e. the maximal number of y-coordinates for any given x-coordinate.

Modular composition: with preinput M, A ∈[x], we input f ∈[x,y] and compute f(x, A)  M.
Our algorithm has online complexity f  f +  A +  M, as long as the bivariate ideal M, y-A is balanced.

We prove that if ⊆^2 is random of fixed cardinality n, and if || ≫ n^2 log(n) then  is balanced with high probability.
Similarly, if M is square-free and A is uniformly random of degree less than M, then M, y-A is balanced with high probability.
Our proof techniques currently do not extend to proving that sheared point sets are balanced.
A few trials we conducted suggest that this may often be the case if the x-valency of  is not too high.
The cost of the second MPE algorithm is not symmetric in the x- and y-degree, so whenever
f <  f one should consider transposing the input, i.e. evaluating f(y, x) on { (β_i, α_i) }_i=1^n.
In this case, the balancedness assumption is on the transposed point set.

Our algorithms are deterministic, and once the preinput has been processed, the user knows whether it is balanced and hence whether the algorithms will perform well.
Further, the performance of our algorithms deteriorates smoothly with how “unbalanced” the preinput is, in the sense of certain polynomials, which depend only on preinput, having sufficiently well behaved degrees.
In a toolbox one might therefore apply our algorithms whenever the preinput turns out to be sufficiently balanced and reverting to other algorithms on very unbalanced preinput.

A typical use of precomputation is if we compute e.g. MPEs on the same point set for many different polynomials.
This occurs in coding theory, where bivariate MPE corresponds to the encoding stage of certain families of codes such as some Reed-Muller codes <cit.> and some algebraic-geometric codes <cit.>: here  is fixed and communication consists of a long series of bivariate MPEs on .
In these applications,  is often not random but chosen carefully, and so our genericity assumptions might not apply.



  
Techniques

We introduce a tool we call reshaping for achieving the following: given an ideal I ⊆[x,y] and f ∈[x,y], compute f̂∈ f + I with smaller y-degree.
For instance in MPE, we let ⊂[x,y] be the ideal of polynomials which vanish on all the points .
Then all elements of f + have the same evaluations on , so we compute a f̂∈ f + of y-degree 0 (it exists if  has distinct x-coordinates), and then apply fast univariate MPE.

An obvious idea to accomplish this iteratively is to find some g ∈ of lower y-degree than f and whose leading y-term is 1, and then compute f̃ = f  g.
The problem is that the x-degree of f̃ may now be as large as f + ( f -  g) g.
Our idea is to seek polynomials g that we call reshapers, which have the form
=2pt
=2pt

    g = y^2d/3 - ĝ ,

where ĝ < d/3 and d =  f+1 (for simplicity, here 3 divides d).
Writing f = f_1 y^2d/3 + f_0 with f_0 < 2d/3,
then f̃ = f_1ĝ + f_0 is easy to compute, has y-degree less than 2d/3, and x-degree only f +  g.
Repeating such a reduction log(d) times with reshapers of progressively smaller y-degree, we eventually reach y-degree 0.

For efficiency, we therefore need the x-degrees of all these reshapers g to be small.
For MPE, stating that g ∈ specifies n linear contraints on the coefficients of ĝ, so we look for g with about n monomials.
Generically, since ĝ≈ d/3, one may expect to find g with g ≈ 3n/d.
Informally,  is balanced if all the reshapers needed in the above process satisfy this degree constraint.

Above, we assumed the point set has distinct x-coordinates.
To handle repetitions, we shear the points by (α, β) ↦ (α + θβ, β), where θ generates an extension field of  of degree 2.
The resulting point set has distinct x-coordinates.
This replaces f(x,y) with f(x - θ y, y), and whenever f <  f we stay within quasi-linear complexity if the sheared point set is balanced.



  
Previous work
Quasi-linear complexity has been achieved for multivariate MPE and interpolation on special point sets and monomial support: Pan <cit.> gave an algorithm on grids, and van der Hoeven and Schost <cit.> (see also <cit.>) generalised this to certain types of subsets of grids, constraining both the points and the monomial support.
See <cit.> for references to earlier work on interpolation, not achieving quasi-linear complexity.

In classical univariate modular composition, we are given f, M, A in [x] and seek f(A)  M.
Brent and Kung's baby-step giant-step algorithm <cit.> performs this operation in n^(ω+1)/2, where ω is the matrix multiplication exponent with best known bound ω < 2.373 <cit.>.
Nüsken and Ziegler <cit.> extended this to a bivariate f, computing f(x, A)  M in complexity f( f)^(ω+1)/2, assuming that A and M have degree at most f  f.

They applied this to solve MPE in the same cost; in this paper,
we use essentially the same link between these problems.
To the best of our knowledge, this is currently the best known cost bound for these problems,
in the algebraic complexity model.

In a breakthough, Kedlaya and Umans <cit.> achieved
“almost linear” time for modular composition and MPE, for specific types of
fields  and in the bit complexity model.
For modular composition, the cost is n^1+ϵ bit operations for any ϵ>0,
while for MPE it is (n + ( f)^2)^1 + ϵ,
assuming f <  f (the algorithm also supports multivariate MPE).
Unfortunately, these algorithms have so far resisted attempts at a practical implementation <cit.>.

Our quasi-linear complexities improve upon the above results (including Kedlaya and Umans' ones since quasi-linear compares favorably to almost linear);
however we stress that none of the latter have the two constraints of our work: allowing precomputation, and genericity assumption.
For modular composition, precomputation on M was suggested in <cit.> to leverage its factorisation structure.
Except for slight benefits of precomputation in Brent and Kung's modular composition (used in the Flint and NTL libraries <cit.>),
we are unaware of previous work focusing on the use of precomputation for MPE, Interpolation, and Modular Composition.

Genericity has recently been used by Villard <cit.>, who showed how to efficiently compute the resultant of two generic bivariate polynomials; a specific case computes, for given univariate M and A, the characteristic polynomial of A in [x] /  M, with direct links to the modular composition f(A)  M <cit.>.
This led to an ongoing work on achieving exponent (ω+2)/3 for modular composition <cit.>.
In that line, the main benefit from genericity is that M, y - A admits bases formed by m polynomials of y-degree <m and x-degree at most (M)/m, for a given parameter 2 ≤ m ≤(M).
Such a basis is represented as an m × m matrix over [x] with all entries of degree at most (M)/m, and one can then rely on fast univariate polynomial matrix algorithms.
In this paper, genericity serves a purpose similar to that in <cit.>: it ensures the existence of such bases for several parameters m, and also of the reshapers g mentioned above; besides we make use of these bases to precompute these reshapers.
Whereas an important contribution of <cit.> is the efficient computation of such bases, here they are only used to find reshapers in the precomputation stage and the speed of computing them is not a main concern.
Once the reshapers are known, our algorithms work without requiring any other genericity property.



  
Organisation After some preliminaries in <ref>, we describe the reshaping strategy for an arbitrary ideal in <ref>.
Then <ref> give algorithms for each of the three problems.
We discuss precomputation in <ref> and genericity in <ref>.



§ PRELIMINARIES



For complexity estimates, we use the algebraic RAM model and count arithmetic operations in .
By (n) we denote the cost of multiplying two univariate polynomials over  of degree at most n;
one may take (n) ∈n log n loglog n⊂n <cit.>.
Division with remainder in [x] also costs (n) <cit.>.
When degrees of a polynomial, say f ∈[x,y], appear in complexity estimates, we abuse notation and let f denote max( f, 1).

It is well-known that univariate interpolation and multi-point evaluation can be done in quasi-linear time <cit.>:
given f ∈[x] and α_1,…,α_n ∈, we may compute (f(α_i))_i=1^n in time ( f + n)log n⊆ f + n;
given α_1,…,α_n and β_1,…,β_n in  with the α_i's pairwise distinct, we may compute the unique corresponding interpolant in time (n)log n⊆n.
We will also use the fact that two bivariate f, g ∈[x,y] can be multiplied in time (d_x d_y)⊂d_x d_y, where d_x = max( f,  g) and d_y = max( f,  g) <cit.>.

For a bivariate polynomial f = ∑_i=0^k f_i(x) y^i ∈[x,y] such that f_k ≠ 0, we define its y-leading coefficient as (f) = f_k ∈[x].

For our genericity results, we will invoke the following staple:

  
  Let f ∈[x_1,…,x_n] be non-zero of total degree d, and ⊆ be finite.
  For α_1,…,α_k ∈ chosen independently and uniformly at random,
  the probability that f(α_1,…,α_k) = 0 is at most d/||.


For a point set ⊆^2, the x-valency of , denoted by (), is the largest number of y-coordinates for any given x-coordinate, i.e.

    () = max_α∈|{β∈| (α,β) ∈}|  .

When () = 1, the x-coordinates of  are pairwise distinct.

The vanishing ideal of  is the bivariate ideal

    () = { f ∈[x,y] | f(α,β) = 0  for all  (α,β) ∈} ,

Hereafter,  stands for the lexicographic order on [x,y] with x  y, and (f) is the -leading term of f∈[x,y].
The following is folklore and follows e.g. from <cit.> and <cit.>.


  
  Let ⊂^2 be a point set of cardinality n and let G = {g_1,…,g_s} be the reduced -basis of (), ordered by .
  Then g_1 ∈[x], and g_s is y-monic with g_s = ().




§ RESHAPE



We first describe our algorithm  which takes f ∈[x,y] and an ideal I and finds f̂∈ f + I whose y-degree is below some target.
This will pass through several intermediate elements of f + I of progressively smaller y-degree.
This sequence of y-degrees has the following form:


  
  We say = (_i)_i=0^k ∈_> 0^k+1 is a (_0, _k)-reshaping sequence if
  _i-1 > _i≥ 2 3 _i-1 for i=1,…,k.
  For I ⊆[x,y] an ideal and = ( _i )_i=0^k a reshaping sequence,
  we say = ( g_i )_i=1^k ∈ I^k is an -reshaper for I if g_i = y^_i + ĝ_i where ĝ_i ≤ 2_i - _i-1, for each i = 1,…,k.


Our algorithms are faster with short reshaping sequences, so we should choose _i ≈2/3_i-1, and hence 2_i - _i-1≈1/3_i.
It is easy to see that for any a,b ∈, there is an (a,b)-reshaping sequence of length less than log_3/2(a) + 2.
Observe that for any (a,b)-reshaping sequence we have η_i ≥2/3 (η_i-1-1) for i=1,…,k and therefore

    2η_i - η_i-1≥η_i-1 - 4/3≥η_i/3 - 1  .

By considering the cases η_i≥ 3 and η_i = 1,2, we get 2η_i - η_i-1≥ 0.

 
  <ref> is correct and has complexity
  
    ∑_i=i_0^k _i ( f + ∑_j=i_0^i   g_j)
       ⊆k  f f + k∑_i=i_0^k _i  g_i ,

  for the smallest i_0 such that _i_0≤ f.



  Let f̂_i, f̂_i,0, f̂_i,1 be the values of f̂, f̂_0, f̂_1 at the end of iteration i.
  First, the iterations for i<i_0 perform no operation and keep f̂_i = f, since _i > f̂_i-1 implies f̂_i,1 = 0 and f̂_i = f̂_i-1.
  In particular, if _i> f for all i then the algorithm is correct and returns f without using any arithmetic operation.
  Now for i≥ i_0, observe that
  f̂_i = f̂_i,1ĝ_i + f̂_i,0
                = f̂_i-1-f̂_i,1 g_i;
  thus in the end f̂∈ f + I since each g_i belongs to I.
  We show the following loop invariants, which imply the degree bounds on the output:
  f̂_i ≤ f + ∑_j=i_0^i  g_j,
     and f̂_i < _i.
  Both are true for i=i_0-1 (just before the loop, if i_0=1).
  For the x-degree, f̂_i = f̂_i-1 - f̂_i,1 g_i yields f̂_i ≤f̂_i-1 +  g_i, and the loop invariant follows.
  For the y-degree, by construction f̂_i,0 < _i and f̂_i,1≤f̂_i-1 - _i hold; the assumption f̂_i-1 < _i-1 then gives f̂_i,1ĝ_i < _i, hence f̂_i < _i.

  For complexity, the only costly step is at <ref> and for iterations i≥ i_0.
  From the above bound f̂_i,1ĝ_i < _i, multiplying f̂_i,1 and ĝ_i costs ((f̂_i,1 + ĝ_i) _i).
  Since ĝ_i =  g_i, since both f̂_i,0 and f̂_i,1 have x-degree at most f̂_i-1, and since f̂_i,0 < _i, the total cost of the ith iteration is in
  
    (f̂_i-1 + ĝ_i)_i⊆( f + ∑_j=i_0^i  g_j)_i.

  Summing over all iterations, we get the first complexity bound in the theorem; the second one follows from it, using the fact that f ≥_i_0 > _i_0+1 > … > _k and i_0≥1.




We now define the balancedness of a point set. In <ref> we prove that this notion captures the expected x-degree of reshapers.

 
  Let ⊆^2 be a point set of cardinality n, and let = ( _i )_i=0^k be a reshaping sequence.
  Then  is -balanced if there exists an -reshaper = (g_i)_i=1^k ∈[x,y]^k for () such that
  g_i ≤n2η_i - η_i-1 + 1 + 1 for i=1,…,k.


The next bound is often used below for deriving complexity estimates; it follows directly from <ref>.

 
  Let = (_i)_i=0^k be a reshaping sequence, ⊆^2 be an -balanced point set of cardinality n, and = (g_i)_i=1^k be an -reshaper for ().
  Then ∑_i=i_0^k _i  g_i ≤ (3n+_i_0)k for 1 ≤ i_0 ≤ k.


We conclude this section with two results about the existence of -reshapers for vanishing ideals of point sets.

 
  Let ⊆^2 be a point set and = ( _i )_i=0^k a reshaping sequence.
  If () ≤min_1≤ i≤ k(2_i - _i-1 + 1), then there exists an -reshaper ∈[x,y]^k for ().


  By <ref>, the reduced -Gröbner basis G of () contains a polynomial with -leading term y^().
  Thus y^ G < () for any , and setting g_i = y^_i - (y^_i G) yields an -reshaper as long as () ≤ 2_i - _i-1 + 1 for all i.



  
  Let ⊆^2 be a point set of cardinality n and a, b ∈ with n > a > b ≥().
  Then there is an (a,b)-reshaping sequence  which satisfies the condition of <ref> and has length k ≤log_3/2(a)+1 ∈ O(log(a)).


  Let v = ()-1 and let ' = ('_0,…,'_k) be any (a-v,b-v)-reshaping sequence with k ≤log_3/2(a-v)+1. Now let = (_0,…,_k) be defined by _i = '_i+v for i = 0,…,k. Then,  is an (a,b)-reshaping sequence. Indeed, clearly the endpoints are correct and _i-1 > _i for i = 1,…,k; moreover,
  
    _i = '_i+v ≥23'_i-1 + v = 23_i-1 + 13 v≥23_i-1 .

  To conclude, we use 2'_i - '_i-1≥ 0 as mentioned above to observe that
  2 _i - _i-1 + 1
    = 2'_i - '_i-1 + v + 1 ≥ v+1 = ().




§ MULTI-POINT EVALUATION



In this section we use reshaping for MPE with precomputation; i.e. given a point set ⊂^2 upon which we are allowed to perform precomputation, and a polynomial f ∈[x,y] which is assumed to be received at online time, compute f(P) for all P ∈.
<ref> deals with the case () = 1, which we reduce to an instance of univariate MPE using .
The cost of <ref> follows directly from <ref>.



  
  <ref> is correct.
  If  is -balanced and  has length in log(n), the complexity is f  f + n.


<ref> can easily be extended to the case where () > 1 by partitioning  into () many subsets, each having x-valency one. This approach also has quasi-linear complexity in the input size as long as () ≪ n, or more precisely if n() ∈n.

When () is large, this strategy is costly, and we proceed instead by shearing the point set, as proposed by Nüsken and Ziegler <cit.>, so that the resulting point set has distinct x-coordinates:
by taking θ∈∖, where  is an extension field of  of degree 2, we apply the map (α,β) ↦ (α + θβ,β) to each element of .
The problem then reduces to evaluating f̅ = f(x - θ y, y) at the sheared points.
To compute f̅, <cit.> provides an algorithm with complexity (d_x(d_x + d_y))log(d_x) using a univariate Taylor shift of f seen as a polynomial in x over the ring [y].
<ref> describes an algorithm for this task which improves the cost on the logarithmic level, by using Taylor shifts of the homogeneous components of f.


a
b

 
  <ref> correctly computes f(  x +   y,y), which has x-degree at most d_x and y-degree at most d_x + d_y,
  at a cost of
  (d_x + d_y)(d_x)log(d_x)⊂d_x(d_x + d_y)
  operations in .


  Observe that y^th_t(x/y) is the homogeneous component of f of degree t, and in particular f = ∑_t = 0^d_x + d_y y^th_t(x/y).
  Thus
  
    f(  x +   y,y)
      = ∑_t=0^d_x + d_yy^t h_t(  x +   y/y)
      = ∑_t=0^d_x + d_yy^t s_t(x/y),

  hence the correctness. The degree bounds on the output are straightforward.
  As for complexity, only <ref> uses arithmetic operations.
  First, scaling h_t(z) ↦ h_t(  z) costs d_x operations in , since h_t ≤ d_x;
  then the Taylor shift h_t(  z) ↦ h_t(  z +  ) costs (d_x)log(d_x) operations in  according to <cit.>.
  Summing over the d_x+d_y iterations yields the claimed bound.



This leads to <ref>, where  may have repeated α_i's.



  <ref> is correct.
  If  is -balanced and  has length in log(n), its complexity is f( f +  f) + n.




§ INTERPOLATION



In this section we use reshaping for the interpolation problem in a similar setting: we input a point set  for precomputation, and input interpolation values at online time.
When  is appropriately balanced, we solve the interpolation problem in quasi-linear time (see <ref>).
The strategy is to first shear the point set to have unique y-coordinates and compute u ∈[y] which interpolates the values on the sheared y-coordinates.
Then we reshape this into r ∈[x,y] with x- and y-degrees roughly √(n).
Shearing back this polynomial to interpolate the original point set is now in quasi-linear time;
a last reshaping allows us to meet the target y-degree.




  <ref> is correct and has complexity
  
    Õ(
          k_1 n + k_2 (√(n) + ∑_j = 1^k_1 g_1,j)^2 + ∑_ℓ=1^2 k_ℓ∑_j=1^k_ℓη_ℓ,k g_ℓ,j).

  If  is _1-balanced and  is _2-balanced, and both _1 and _2 have length in log n, then the complexity is n.


  First note that a reshaping sequence of length O(log n) and satisfying the preinput constraints exists, due to <ref> and the assumption d ≥().
  For correctness, observe that all points in  have pairwise distinct y-coordinates, so computing u makes sense.
  Viewing u as an element of [x,y] with u = 0, we have u(α_i,β̅_i) = _i.
  By <ref> then r has the same evaluations and r < √(n) and r ≤∑_i=1^k_1 g_1,i.

  Then, in both cases ()=1 and ()>1, we have
  
    _i = r(α_i,β̅_i) = s(α_i,β_i) = s_1(α_i,β_i) + θ s_2(α_i,β_i)

  for i = 1,…,n. Since s_1,s_2 ∈[x,y] and all _i's are in , we get s_2(α_i,β_i) = 0 and s_1(α_i,β_i) = γ_i for i = 1,…,n.
  We also then have that s_1 ≤ s < √(n) and
  
    s_1 ≤ s ≤ r +  r ≤√(n) + ∑_j = 1^k_1 g_1,j .

  Thus, by <ref> again, the output f is such that f(α_i,β_i) = _i for i = 1, …, n, and f < d, and
  
    f ≤√(n) + ∑_j=1^k_1 g_1,j + ∑_j=1^k_2 g_2,j .

  The complexity bound gathers the calls to <ref>,
  and the relaxed cost assuming balancedness is due to <ref>.




§ MODULAR COMPOSITION



We now turn to the following modular composition problem: given M, A ∈[x] with n :=  M >  A, and f ∈[x,y], compute

    f(x, A(x))  M(x) ∈[x]  .

We consider the variant of the problem where M and A are available for precomputation.
Computing (<ref>) is tantamount to computing the unique element of (f + I) ∩[x] of degree less than n,
for the ideal I = M, y-A⊆[x,y].
One can thus see this as a reshaping task: given f of some y-degree, reshape it to a polynomial of y-degree 0 while keeping it fixed modulo I: this is formalised as <ref>.



Like for point sets above, if = (_i)_i=0^k is a reshaping sequence, we say that I = M, y-A is -balanced if there exists an -reshaper = (g_i)_i=1^k for I such that g_i ≤n/2_i - _i-1+1 + 1.


  <ref> is correct.
  If M, y-A is -balanced and  has length in log(n), the complexity is f  f + n.




§ PRECOMPUTING RESHAPERS





 §.§ Reshapers for general ideals


Here we describe <ref> for precomputing reshapers for any zero-dimensional ideal I ⊆[x,y], given a -basis of I.
It operates through the [x]-module I_δ := { f ∈ I | f < δ}, so we first expound the relation between this and I as a corollary of Lazard's structure theorem on bivariate -bases <cit.>.


  
  Let G = { b_0,…,b_s }⊂[x,y] be a minimal -basis defining an ideal I = G.
  For δ∈, let I_δ = { f ∈ I | f < δ}, let
  ŝ = max{ i | b_i < δ,  0 ≤ i ≤ s }, let d_i =  b_i for 0 ≤ i ≤ŝ and d_ŝ+1 = δ.
  Then I_δ is a [x]-submodule of [x,y]_ < δ which is free of rank δ - d_0 and admits the basis
  { y^j b_i  |  0 ≤ j < d_i+1 - d_i, 0 ≤ i ≤ŝ}.

A proof is given in appendix.
We will use the following [x]-module isomorphism which converts between bivariate polynomials of bounded y-degree and vectors over [x]:
for any δ∈_>0,

    _δ: f = ∑_j=0^δ-1 f_j(x)y^j ∈[x,y] ↦ [f_0,…,f_δ-1] ∈[x]^1 ×δ .

If I is zero-dimensional then in <ref> we have d_0 = 0 and I_δ has rank δ.
Any basis B of I_δ can be represented as a nonsingular matrix M_B ∈[x]^δ×δ whose rows are _δ(B).
Then, (I_δ) := (M_B) does not depend on the choice of B since all bases of I_δ have the same determinant up to scalar multiplication.

In this section, we use the Popov form <cit.>, which can be defined for any matrix and with “shifts”; here we only need the unshifted, nonsingular square case.

  For any row vector v⃗∈[x]^1 ×δ its row degree denoted v⃗ is the maximal degree among its entries.
  The pivot of v⃗ is the rightmost entry of v⃗ with degree v⃗.
  A nonsingular matrix P = [ p_ij ] ∈[x]^δ×δ is in Popov form if p_ii is the pivot of the ith row, is monic, and p_ii >  p_ji for any j ≠ i.


For a (free) [x]-submodule ⊂[x]^1 ×δ of rank δ, we identify a basis of  as the rows of a nonsingular matrix in [x]^δ×δ.
Any such  has a unique basis P ∈[x]^δ×δ in Popov form, which we call the Popov basis of .
It has minimal row degrees in the following sense:
if N ∈[x]^δ×δ is another basis of , there is a bijection ψ from the rows of P to the rows of N such that p⃗≤ψ(p⃗) for any row p⃗ of P.
The Popov basis satisfies () = (P) = |(P)|, using the following notation:
the sum of the entries of a tuple t⃗∈_≥ 0^δ is denoted |t⃗|;
the column degree of a matrix B ∈[x]^δ×δ is (B) = (d_i)_i=1^δ∈_≥ 0^δ, with d_i the largest degree in the ith column of B (for a zero column, d_i =0).

The next result allows us to compute Popov forms efficiently.


  
  There is an algorithm which inputs a nonsingular matrix B ∈[x]^δ×δ and outputs the Popov basis of the [x]-row space of B using δ^-1 |(B)| operations in , assuming that δ∈|(B)|.


Since Popov forms are “column reduced”, they are well suited for matrix division with remainder <cit.>:
if P ∈[x]^δ×δ is the Popov basis of , then for any v⃗∈[x]^1 ×δ there is a unique u⃗∈v⃗ + such that (u⃗) < (P) entrywise; we denote u⃗ = v⃗ P.
Furthermore, u⃗ has minimal row degree among all vectors in v⃗ +.
Such remainders can be computed efficiently:


  
  There is an algorithm which inputs a Popov form P ∈[x]^δ×δ and v⃗∈[x]^1 ×δ such that (v⃗) < (P) + ((P),…,(P)) entrywise, and outputs v⃗ P using δ^-1(P) operations in , assuming that δ∈(P).





  
  <ref> is correct.
  Assuming ∈(I_δ), it costs δ^ω-1(I_δ) +  s  b_0 operations in .



  Since G is a -basis, if y^+I contains a polynomial of y-degree less than δ, then (y^η G) ≤δ and the algorithm does not fail at <ref>.

  For correctness of the output, observe that y^η - R ∈ I so satisfactory g = y^η - g̃ all have g̃∈ R + I_δ.
  Now, ĝ of <ref> is clearly in R + I_δ since P is the Popov basis of I_δ, but also ĝ has minimal x-degree in the coset R + I_δ.
  Hence among all g of the correct form, the algorithm returns that of minimal x-degree.

  For complexity, work is done in <ref>.
  Since G is reduced, b_0 > … >  b_s.
  Therefore the diagonal entries in B are dominant in their columns and | B| = (B) = (P) = (I_δ).
  For <ref>, we use the algorithm of <cit.> with cost η s  b_0, see <ref>.
  <ref> costs δ^ω-1| B| by <ref> and <ref> costs δ^ω-1(P) since R <  b_0 < (P).




 §.§ Reshapers for the considered problems



We turn to obtaining the reduced -basis of ().
We will consider the [x]-submodule _m() = () ∩[x,y]_ < m which by <ref> is free and of rank m.
To obtain a -basis, our approach is to first compute the Hermite basis of _m().
This is the unique basis whose corresponding matrix H ⊂[x]^m × m is lower triangular, with each diagonal entry monic and strictly dominating the degrees in its column.


  
  For any point set ⊆^2 and any m > (), we have () = _m() and (_m()) = ||.


  By <ref> the elements of the reduced -Gröbner basis of () have y-degree at most (), implying the first claim.
  Further, this means the quotient [x,y] / () is isomorphic to the quotient of modules [x,y]_ < m / _m().
  It is a basic property of zero-dimensional varieties that the -dimension of the former is the number of points in , which is hence also the -dimension of the latter.
  This dimension is (_m()) by <cit.>.



  
  There is an algorithm which inputs ⊂^2 and outputs the reduced -basis of () and has complexity ()^-1 ||.


  Let = (), _m = _m(), and m = () + 1.
  We first compute the Hermite basis H of _m() in time m^-1 || using (a special case of) <cit.>, in which taking s⃗ = (0,n,…,(m-1)n) ensures that the s⃗-Popov basis P of _m is the Hermite basis.

  Let G = { g_0,…,g_m-1}⊂[x,y] be given as the _m^-1-image of the rows of H.
  By <ref> and since H is lower triangular, G is a -basis of  but not necessarily minimal.
  Construct G' ⊆ G from G by excluding the elements g ∈ G such that there is g' ∈ G with g' <  g and ((g')) ≤((g)), i.e. (g') divides (g).
  This makes G' a minimal -Gröbner basis of  <cit.>, and we claim it is the reduced one.
  Indeed, since H is in Hermite form, the selection criteria for G' ensures that for any g ≠ g' in G' and any term x^iy^j in g', we have i < ((g)) or j <  g, and hence G' is reduced.
  Obtaining G' from H costs no arithmetic operations.



  
  Given a point set ⊆^2 of cardinality n and a reshaping sequence =(_i)_i=0^k with n ≥_k and satisfying the condition of <ref>, then we can determine if  is -balanced and compute an -reshaper = (g_i)_i=1^k for  where each element has minimal possible x-degree in complexity k_0^ω-1 n + _0  n k.


  By <ref>, computing a reduced -basis G = (b_i)_i=0^ of () costs ν_x^-1 n⊂_0^-1 n.
  We then run <ref> on input = _i and δ_i = 2_i - _i-1 + 1 > for i=1,…,k.
  <ref> ensures (_δ()) = n for any δ >,
  thus the cost of each call to <ref> becomes _0^ω-1 n + _0  n.



  Given M, A ∈[x] with n :=  M >  A and a reshaping sequence =(_i)_i=0^k with n ≥_k, then we can determine if I := M, y - A is -balanced and compute an -reshaper = (g_i)_i=1^k for  where each element has minimal possible x-degree in complexity k _0^ω-1 n.


  For any δ, and using the notation of <ref>, the basis B of I_δ is lower triangular with diagonal entries (M,1,…,1).
  Hence (B) = (I_δ) = n.
  Using s = 1 and b_0 =  M = n, the cost follows from <ref>.




§ GENERICITY



Now we show that on random input our algorithms usually have quasi-linear complexity, i.e. that random point sets are balanced and that M, y-A is balanced for random univariate A, M.


  
  Let α_1,…,α_n ∈ be distinct, let y_1,…,y_n be new indeterminates,
  and consider for s ∈ the matrix
  
    A_s = [ V_s | D V_s |…| D^m-1 V_s ] ∈[y_1,…,y_n]^n × ms

  where D is the diagonal matrix with entries (y_1,…,y_n), and V_s = [ α_i^j-1 ]_1≤ i≤ n,1 ≤ j≤ s∈^n × s.
  Then A_s has rank min(n, ms).


  Note that by rank of a matrix over [y_1,…,y_n], we mean the rank of that matrix seen as over the field of fractions (y_1,…,y_n).
  If we specialise y_i to α_i^s for i = 1, …, n, we obtain the Vandermonde matrix Â_s = [ α_i^j-1 ]_1 ≤ i≤ n,1 ≤ j ≤ ms∈^n× ms of the points α_1,…,α_n.
  Since these points are distinct, Â_s has full rank min(n, ms).
  Hence A_s must also have full rank.


The columns of A_s can be identified to monomials x^i y^j with i<s and j<m.
In particular, if p ∈() is a bivariate polynomial with x-degree less than s and y-degree less than m which vanishes on a point set = { (α_i,β_i) }_i=1^n ⊂^2 with distinct α_i's, then the coefficients of p form a vector in the right kernel of the matrix Â_s = (A_s)_|y_i →β_i∈^n × ms specializing y_i to β_i.

The next lemma determines the exact row degrees of the Popov basis P ∈[x]^m × m of _m(_m()) for a “random” point set , where _m() = () ∩[x,y]_ < m as in <ref>.


  
  Let α_1,…,α_n ∈ be distinct, let ⊆ be a finite subset, and let λ: ^n →^n be an affine map.
  For γ_1,…,γ_n ∈ chosen independently and uniformly at random, set = { (α_i,β_i) }_i=1^n where (β_1,…,β_n) = λ(γ_1,…,γ_n).
  Let m ∈ with () < m ≤ n and let (d,t) = (n, m).
  With probability at least 1 - 2nm/||, the Popov basis P ∈[x]^m × m of _m(_m()) has exactly m-t rows of degree d and t rows of degree d+1 and
  in particular _x P ≤ d+1.


  Let p_1,…,p_m ∈[x,y] be the polynomials defined by the rows of P.
  <ref> shows (P) = n = ∑_i=1^m _x p_i.

  For any s ∈, let A_s ∈[y_1,…,y_n]^n × ms be as in <ref>, hence (A_s) = min(n, ms).
  Let Â_s = (A_s)_|y_i →β_i∈^n × ms.
  Taking s = d, as mentioned above, if _x p_i < d for some i, then the coefficient vector of p_i is in the right kernel of Â_d, and so (Â_d) < (A_d) = md ≤ n.
  Thus, letting M ∈[y_1,…,y_n] be a non-zero md × md minor of A_d then M(β_1,…,β_n) = M(λ(γ_1,…,γ_n)) = 0;
  M has degree at most m-1 in each variable, so the total degree of M is less than nm, and since λ is affine the composition M(λ(z_1,…,z_n)) also has total degree less than nm.
  Then, by <ref> the probability that M(λ(γ_1,…,γ_n)) = 0 is at most nm/||.

  Assume now that all rows of P have degree at least d.
  For each i such that p_i = d, the coefficients of p_i form a vector in the right kernel of Â_d+1∈^n × m(d+1).
  By <ref>, A_d+1 has a right kernel (over the fractions) of dimension m(d+1) - n = m-t.
  Since the rows of P are linearly independent over [x], and therefore also over , whenever (Â_d+1) = (A_d+1) at most m-t rows of P have x-degree d.
  We thus consider N ∈[y_1,…,y_n] a non-zero n × n minor of A_d+1.
  Again N has total degree less than nm and the probability that N(β_1,…,β_n) = N(λ(γ_1,…,γ_n)) = 0 is at most nm/||, bounding the probability that (Â_d+1) < (A_d+1).

  Hence, with probability at least 1 -2nd/||, P has all rows of degree at least d and j rows of degree exactly d with j ≤ m-t.
  Each of the remaining m-j rows has degree at least d+1, while their degrees must sum to
  n - jd = md + t - jd = (m-j)d + t ≤ (m-j)(d+1).
  Hence each of them has degree exactly d+1.


<ref> for computing reshapers outputs a g = y^η - ĝ with ĝ < δ satisfying ĝ≤ P, where P is the Popov basis of _δ().
<ref> states that generically we can expect P ≤n/δ + 1, and so when δ = 2η_i - η_i-1 + 1 in a reshaping sequence, this matches the definition of -balanced.


  
  Let α_1,…,α_n ∈ be distinct, let ⊆ a finite subset, and let λ: ^n →^n be an affine map.
  For γ_1,…,γ_n ∈ chosen independently and uniformly at random, set = { (α_i,β_i) }_i=1^n where (β_1,…,β_n) = λ(γ_1,…,γ_n).
  Let = ( _i )_i=0^k be a reshaping sequence with _0 ≤ n and satisfying the constraint of <ref>.
  Then  is -balanced with probability at least 1 - n^2k/||.


The above proposition directly applies to both our MPE and interpolation algorithms on random point sets with unique x-coordinates.
Note that in the case of interpolation, where the point set is sheared if its y-valency is greater than one, the property of being -balanced is not inherited a priori by the sheared point set. The probability of being -balanced, however, is preserved, since the shearing acts as an affine transformation on the y-coordinates.
There are many formulations depending on the type of randomness one needs over the point sets; the following is a simple example over finite fields:


  
  Let d, n ∈ with d ≤ n and q be a finite field with q elements, and let = { (α_i,β_i) }_i=1^n ⊆ q^2 be chosen uniformly at random among point sets with cardinality n.
  Then with probability of at least (1 - n^2/q)(1 - 3n^2(log_3/2(n)+1)/q) over the choice of  the following two problems can be solved with cost n:
  
    
  * Input polynomial f ∈ q[x,y] such that f < n/d and f < d, and output ( f(α_i, β_i) )_i=1^n ∈ q^n.
    
  * Input interpolation values γ⃗= ( γ_i )_i=1^n ∈ q^n, and output f ∈ q[x,y] satisfying f(α_i, β_i) = γ_i for i=1,…,n, as well as _y f < d and _x f ≤ cn for some constant c which depends only on n and d.
  


  The probability simply bounds the probability that  has unique x-coordinates and that it is balanced in all the necessary ways.
  By <ref> there is an appropriate reshaping sequence of length at most log_3/2(n) + 2.


We do not make a claim about the genericity in <ref>: due to the shearing in that algorithm, the arguments of this section do not immediately apply.
Lastly, we turn to modular composition.

  Let M ∈[x] be square-free of degree n and let  be a (d, 1)-reshaping sequence of length k with 0 < d ≤ n.
  Let ⊆ be a finite subset, and let A = ∑_i=0^n-1 a_i x^i-1∈[x] where a_0,…,a_n-1 are chosen independently and uniformly at random from .
  Then M, y - A is -balanced with probability at least 1 - n^2k/||.


  Let  be the splitting field of M, so M = ∏_i=1^n (x-α_i) for some pairwise distinct α_1,…,α_n ∈.
  Define the stochastic variables β_i = A(α_i) for i=1,…,n;
  the map λ(a_0,…,a_n-1) = (β_1,…,β_n) is -linear.
  Consider = { (α_i, β_i) }_i=1^n ⊆^2.
  Then <ref> implies that  is -balanced with probability at least 1 - n^2k/||.
  In this case, for each i there exists g_i = y^_i + ĝ_i ∈ I_ where ĝ_i < 2_i - _i-1 and ĝ_i ≤n/2_i - _i-1+1 + 1, and where I_ = M, y-A⊗_.
  Let {1,ζ,…,ζ^s-1}⊂ be a basis of : and write g_i = g_i,0 + ζ  g_i,1 + … + ζ^s-1  g_i,s-1 with g_i,j∈[x,y].
  Then g_i ∈ I_ implies that g_i,0∈ I, and by the shape of g_i then g_i,0 = y^η_i +  ĝ_i,0 where the x- and y-degree of ĝ_i,0 satisfy the same bounds as ĝ_̂î.
  Then the tuple g⃗_0 = (  g_1,0, …, g_k,0 ) ∈[x,y]^k forms a balanced -reshaper for I.


ACM-Reference-Format

31


















#1     #1#1#1     #1     #1     #1     #1#1          #1#1           






[Assmus and KeyAssmus and Key1992]
        assmus_designs_1992
authorpersonE. F. Assmus and personJ. D.
  Key. year1992.
booktitleDesigns and Their Codes.
publisherCambridge University Press.
978-0-521-41361-9

<https://doi.org/10.1017/CBO9781316529836>



[Brent and KungBrent and Kung1978]
        brent-kung-1978
authorpersonR. P. Brent and personH. T.
  Kung. year1978.
Fast Algorithms for Manipulating Formal Power
  Series.
journalJ. ACM volume25,
  number4 (year1978), pages581–595.
0004-5411

<https://doi.org/10.1145/322092.322099>



[Cantor and KaltofenCantor and
  Kaltofen1991]
        cantor_fast_1991
authorpersonD. G. Cantor and personE.
  Kaltofen. year1991.
On fast multiplication of polynomials over
  arbitrary algebras.
journalActa Informatica volume28,
  number7 (year1991), pages693–701.
0001-5903, 1432-0525

<https://doi.org/10.1007/BF01178683>



[Cox, Little, and O'SheaCox
  et al2015]
        CoxLittleOShea2015
authorpersonD. A. Cox, personJ. Little,
  and personD. O'Shea. year2015.
booktitleIdeals, Varieties, and Algorithms
  (edition4th ed.).
publisherSpringer.

<https://doi.org/10.1007/978-3-319-16721-3>



[CoxonCoxon2018]
        coxon_fast_2018
authorpersonN. Coxon. year2018.
Fast systematic encoding of multiplicity codes.
journalJ. Symb. Comput. (year2018).
0747-7171

<https://doi.org/10.1016/j.jsc.2018.08.005>



[DahanDahan2009]
        Dahan2009
authorpersonX. Dahan. year2009.
Size of Coefficients of Lexicographical Gröbner
  Bases: The Zero-Dimensional, Radical and Bivariate Case. In
  booktitleProceedings ISSAC 2009.
  pages119–126.

<https://doi.org/10.1145/1576702.1576721>



[DeMillo and LiptonDeMillo and
  Lipton1978]
        DeMilloLipton78
authorpersonR. A. DeMillo and personR. J.
  Lipton. year1978.
A Probabilistic Remark on Algebraic Program
  Testing.
journalInf. Process. Lett. volume7,
  number4 (year1978), pages193–195.

<https://doi.org/10.1016/0020-0190(78)90067-4>



[Hart, Johansson, and PancratzHart
  et al2015]
        flint
authorpersonW. Hart, personF. Johansson,
  and personS. Pancratz. year2015.
titleFLINT: Fast Library for Number Theory.
Version 2.5.2, <http://flintlib.org>.


[Jeannerod, Neiger, Schost, and
  VillardJeannerod et al2016]
        jeannerod_fast_2016
authorpersonC.-P. Jeannerod, personV.
  Neiger, personÉ. Schost, and personG. Villard.
  year2016.
Fast Computation of Minimal Interpolation
  Bases in Popov Form for Arbitrary Shifts. In
  booktitleProceedings ISSAC 2016.
  pages295–302.
978-1-4503-4380-0

<https://doi.org/10.1145/2930889.2930928>



[KailathKailath1980]
        kailath_linear_1980
authorpersonT Kailath.
  year1980.
booktitleLinear Systems.
publisherPrentice-Hall.



[Kedlaya and UmansKedlaya and Umans2011]
        kedlaya_fast_2011
authorpersonK. Kedlaya and personC.
  Umans. year2011.
Fast Polynomial Factorization and Modular
  Composition.
journalSIAM J. Comput. volume40,
  number6 (dateJan. year2011),
  pages1767–1802.
0097-5397

<https://doi.org/10.1137/08073408X>



[LazardLazard1985]
        lazard_ideal_1985
authorpersonD. Lazard.
  year1985.
Ideal bases and primary decomposition: case of two
  variables.
journalJ. Symb. Comput. volume1,
  number3 (year1985).

<https://doi.org/10.1016/S0747-7171(85)80035-3>



[Le GallLe Gall2014]
        le_gall_powers_2014
authorpersonF. Le Gall.
  year2014.
Powers of tensors and fast matrix multiplication.
  In booktitleProceedings ISSAC 2014.
  publisherACM, pages296–303.

<https://doi.org/10.1145/2608628.2608664>



[MiuraMiura1993]
        miura_algebraic_1993
authorpersonS. Miura. year1993.
Algebraic geometric codes on certain plane curves.
journalElectronics and Communications in Japan (Part
  III: Fundamental Electronic Science) volume76,
  number12 (year1993), pages1–13.
1520-6440

<https://doi.org/10.1002/ecjc.4430761201>



[Neiger, Salvy, Schost, and VillardNeiger
  et al2020]
        NeigerSalvySchostVillard
authorpersonV. Neiger, personB. Salvy,
  personÉ. Schost, and personG. Villard.
  year2020.
titleFaster modular composition (work in progress).



[Neiger and VuNeiger and Vu2017]
        neiger_computing_2017
authorpersonV. Neiger and personT. X.
  Vu. year2017.
Computing Canonical Bases of Modules of
  Univariate Relations. In booktitleProceedings ISSAC
  2017.

<https://doi.org/10.1145/3087604.3087656>



[Nüsken and ZieglerNüsken and
  Ziegler2004]
        nusken-ziegler-2004
authorpersonM. Nüsken and personM.
  Ziegler. year2004.
Fast Multipoint Evaluation of Bivariate
  Polynomials. In booktitleProceedings ESA 2004.
978-3-540-30140-0

<https://doi.org/10.1007/978-3-540-30140-0_49>



[PanPan1994]
        pan_simple_1994
authorpersonV. Y. Pan.
  year1994.
Simple Multivariate Polynomial
  Multiplication.
journalJ. Symb. Comput. volume18,
  number3 (year1994), pages183–186.
0747-7171

<https://doi.org/10.1006/jsco.1994.1042>



[Paterson and StockmeyerPaterson and
  Stockmeyer1973]
        paterson_number_1973
authorpersonM. S. Paterson and personL. J.
  Stockmeyer. year1973.
On the number of nonscalar multiplications
  necessary to evaluate polynomials.
journalSIAM J. Comput. volume2,
  number1 (year1973), pages60–66.

<https://doi.org/10.1137/0202007>



[PopovPopov1970]
        popov_properties_1970
authorpersonV Popov. year1970.
Some properties of the control systems with
  irreducible matrix-transfer functions. In booktitleSeminar
  on Diff. Eq. and Dyn. Sys., II. pages169–180.

<https://doi.org/10.1007/BFb0059934>



[SchwartzSchwartz1980]
        Schwartz80
authorpersonJ. T. Schwartz.
  year1980.
Fast Probabilistic Algorithms for Verification of
  Polynomial Identities.
journalJ. ACM volume27,
  number4 (year1980), pages701–717.

<https://doi.org/10.1145/322217.322225>



[ShoupShoup2020]
        ShoupNTL
authorpersonV. Shoup. year2020.
titleNTL: A Library for doing Number Theory, version
  11.4.3.
howpublished.



[van der Hoevenvan der Hoeven2015]
        van_der_hoeven_complexity_2015
authorpersonJ. van der Hoeven.
  year2015.
On the complexity of multivariate polynomial
  division. In booktitleProceedings ACA 2015.
  pages447–458.

<https://doi.org/10.1007/978-3-319-56932-1_28>



[van der Hoeven and Lecerfvan der Hoeven and
  Lecerf2018]
        van_der_hoeven_modular_2018
authorpersonJ. van der Hoeven and personG.
  Lecerf. year2018.
Modular composition via factorization.
journalJ. Complexity  volume48
  (year2018), pages36–68.
0885-064X

<https://doi.org/10.1016/j.jco.2018.05.002>



[van der Hoeven and Lecerfvan der Hoeven and
  Lecerf2019]
        van_der_hoeven_fast_2019
authorpersonJ. van der Hoeven and personG.
  Lecerf. year2019.
Fast multivariate multi-point evaluation
  revisited.
journalJ. Complexity (year2019).
0885-064X

<https://doi.org/10.1016/j.jco.2019.04.001>



[van der Hoeven and Schostvan der Hoeven and
  Schost2013]
        van_der_hoeven_multi-point_2013
authorpersonJ. van der Hoeven and
  personÉ. Schost. year2013.
Multi-point evaluation in higher dimensions.
journalAppl. Algebra Eng. Commun. Comput.
  volume24, number1 (year2013),
  pages37–52.

<https://doi.org/10.1007/s00200-012-0179-3>



[VillardVillard2018a]
        villard_computing_2018
authorpersonG. Villard.
  year2018a.
On computing the resultant of generic bivariate
  polynomials. In booktitleProceedings ISSAC 2018.
  pages391–398.

<https://doi.org/10.1145/3208976.3209020>



[VillardVillard2018b]
        villard_computing_2018_slides
authorpersonG. Villard.
  year2018b.
titleOn computing the resultant of generic bivariate
  polynomials. Presentation at ISSAC 2018.

<http://www.issac-conference.org/2018/slides/villard-computingresultant.pdf>




[von zur Gathenvon zur Gathen1990]
        Gathen1990
authorpersonJ. von zur Gathen.
  year1990.
Functional decomposition of polynomials: the tame
  case.
journalJ. Symb. Comput. volume9,
  number3 (year1990).
0747-7171

<https://doi.org/10.1016/S0747-7171(08)80014-4>



[von zur Gathen and Gerhardvon zur Gathen and
  Gerhard2013]
        von_zur_gathen_modern_2012
authorpersonJ. von zur Gathen and personJ.
  Gerhard. year2013.
booktitleModern Computer Algebra
  (edition3rd ed.).
publisherCambridge University Press.

<https://doi.org/10.1017/CBO9781139856065>



[ZippelZippel1979]
        Zippel79
authorpersonR. Zippel.
  year1979.
Probabilistic algorithms for sparse polynomials.
  In booktitleProceedings EUROSAM'79.
  pages216–226.

<https://doi.org/10.1007/3-540-09519-5_73>








§ APPENDIX




  
  Let G = { b_0,…,b_s}⊂[x,y] be a minimal -basis, sorted according to .
  Then
  
    
  * b_0 < … <  b_s; and
    
  * (b_s) |(b_s-1) |⋯|(b_0).
  



  Since I is an ideal of [x,y] and I_δ = I ∩[x,y]_<δ, then I_δ is a [x]-submodule of [x,y]_<δ.
  Let ℬ denote the (claimed) basis in the corollary. Clearly ℬ⊆ I_δ, and the elements of ℬ all have different y-degree and so are [x]-linearly independent.
  Also |ℬ| = δ - d_0, so if ℬ generates I_δ then it is a basis of it and the rank of I_δ is δ - d_0.
  It remains to show that ℬ generates I_δ, so take some f ∈ I_δ.
  Since f ∈ I the multivariate division algorithm using G and the order  results in q_0,…,q_s ∈[x,y] such that f = q_0 b_0 + … + q_s b_s with q_i ≤ f -  b_i.
  Since f < δ this means q_ŝ+1 = … = q_s = 0.
  Say that in each iteration of the division algorithm, we use the greatest index i for which (b_i) divides the leading term of the current remainder.
  Thus no term of q_i b_i is divisible by (b_i+1) for any i < s.
  But by <ref> then (b_i+1) divides (b_i), and so if (q_i b_i) ≥ b_i+1 then (b_i+1) |(q_i b_i).
  Consequently q_i <  b_i+1 -  b_i, and therefore f is a [x]-linear combination of the elements of ℬ.



  
  There is an algorithm which inputs a -basis G = [ b_0,…,b_s] ⊆[x,y] with b_0 = 0, and a polynomial f ∈[x,y], and outputs f  G in time |G|d_x( f), where d_x = max( f,  b_0).


  ^∘_x
  ^∘_y
  This is a special case of <cit.>: the multivariate division algorithm computes q_0,…,q_s,R ∈[x,y] such that f = q_0 b_0 + … + q_s b_s + R with R = f  G, and the cost of the algorithm can be bounded as
  
    ∑_i=0^s (q_i b_i)(q_i b_i) + (R)(R)  ,

  where (·) denotes an a priori upper bound on the x-degree, and similarly for (·).
  Since G is a -basis, then (q_i b_i) ≤ f and (R) ≤ f.
  For the x-degrees, note that in an iteration of the division algorithm where b_i, i > 0 is used, then R̃ <  b_0, where R̃ is the current remainder, since otherwise the algorithm would have reduced by b_0 as b_0 = 0.
  Hence (q_i) ≤(q_i (b_i)) <  b_0 and so (q_i b_i) ≤ 2 b_0.
  Similarly, (R) <  b_0.
  Left is only (q_0b_0): since q_0 b_0 = f - q_1 b_1 - … - q_s b_s - R, then (q_0 b_0) ≤max( f, 2 b_0).


