> <\body> ||<\doc-note> Most of this work was achieved in 2019 and presented in the preprint. During that period, we note that Joris van der Hoeven was working at the #3069 of the CNRS. |||<\author-affiliation> CNRS, École polytechnique, Institut Polytechnique de Paris Laboratoire d'informatique de l'École polytechnique (LIX, UMR 7161) Bâtiment Alan Turing, CS35003 1, rue Honoré d'Estienne d'Orves 91120 Palaiseau, France |>>||<\author-affiliation> CNRS, École polytechnique, Institut Polytechnique de Paris Laboratoire d'informatique de l'École polytechnique (LIX, UMR 7161) Bâtiment Alan Turing, CS35003 1, rue Honoré d'Estienne d'Orves 91120 Palaiseau, France |>>||>> Consider a multivariate polynomial K,\,x|]>> over a field , which is given through a black box capable of evaluating at points in >, or possibly at points in> for any -algebra . The problem of sparse interpolation is to express in its usual form with respect to the monomial basis. We analyze the complexity of various old and new algorithms for this task in terms of bounds and for the total degree of and its number of terms. We mainly focus on the case when is a finite field and explore possible speed-ups. > Consider a polynomial function \K> over a field given through a black box capable of evaluating at points in >. The problem of sparse interpolation is to recover the representation of K,\,x|]>> in its usual form, as a linear combination <\equation> f=,\,i>f,\,i>*x>*\*x> of monomials. The aim of this paper is to analyze various approaches for solving this problem, with our primary focus on the case when is a finite field. We will survey and synthesize known methods, but we will also present a few new algorithms, together with improved complexity bounds for some important special cases. Efficient algorithms for the task of sparse interpolation go back as far as to the eighteen's century and the work of Prony. The first modern version of the algorithm is due to Ben Or and Tiwari. This method was swiftly embraced in computer algebra; for early implementations, we refer to. There has been aregain of interest for the problem during the last decade, both from atheoretical perspective and from the practical point of view. We also mention the survey paper by Roche on the more general topic of computations with sparse polynomials and the recent PhD thesis by Perret du Cray for the currently best known algorithm for univariate polynomials with integer coefficients; see3.18> or1.2>. Putting the most efficient algorithms into practice constitutes an important challenge for computer algebra systems. The ultimate implementation should be able to interpolate sparse polynomials almost as fast as dense ones of the same size. This motivates us to explore various methods under heuristic conditions that we expect to fairly reflect average behavior in practice. Sometimes we will prefer a relaxed and intuitive style of exposition to mathematically precise theorems with rigorous proofs. Throughout this paper stands for the total degree of and for the number of nonzero terms in(). Whenever available, the uppercase characters and represent upper bounds for and . We also write for the number of ring or field operations in that are required in order to evaluate at a point. The complexity analysis of sparse interpolation has to be carried out with a lot of care, due to the large variety of cases that can occur: <\itemize> What kind of complexity/evaluation model do we use? <\itemize> Do we count the number operations in or the number of bit operations? Are we interested in theoretic (asymptotic) or practical complexity? Are divisions allowed for the evaluation of and how do we count them? Are we only allowed to evaluate at points in > or also at points in > for certain extension rings or fields >? What kind of coefficient field do we use? <\itemize> A field from analysis such as >. A discrete field such as > or a finite field >. Fields with roots of unity > of large smooth order in . The univariate case () the multivariate case (2>). Informally speaking, there are three levels of \Psparsity\Q: <\itemize> Weakly sparse: total degrees of the order >; Normally sparse: total degrees of the order >>; Super sparse: total degrees of order with >. We also notice that almost all general algorithms for sparse interpolation are probabilistic of Monte Carlo type. Indeed, without further knowledge about , such as its support or its number of terms, the mere knowledge of a finite number of evaluations of only allows us to guess plausible expressions for . In this paper, we will be mostly interested in the practical bit complexity of sparse interpolation over finite fields >. Sparse interpolation over the rational numbers can often be reduced to this case as well, in which case is a well chosen prime number that fits into 32 or 64 bits and such that admits a large smooth divisor; see section. We analyze the complexities of specializations of existing algorithms to the finite field case and also present afew new algorithms and tricks for this specific setting. Due to the large number of cases that can occur, we will not prove detailed complexity bounds for every single case, but rather outline how various ideas may be used and combined to reduce the practical complexity. From our practical perspective, it is important to take into account logarithmic factors in complexity bounds, but it will be convenient to ignore sub-logarithmic factors. For this reason, we use the notation <\equation*> \=O>|)>|def>\=O*|)>>*|)>>|)> for any functions ,\>. We will also write > for the bit cost of multiplying two polynomials of degree over and abbreviate \>>. For instance, the naive multiplication algorithm yields =O>*log q|)>>. For our complexity analyses we will give priority to the asymptotic complexity point of view and use the well known bound =O>>. Many of the challenges concerning sparse interpolation already arise in the univariate case when . As we will see in section, the multivariate case can actually be reduced to the univariate one using the technique called \PKronecker segmentation\Q, even though other approaches may be more efficient. For this reason, a large part of the paper is devoted to methods for interpolating a univariate black box function >. We distinguish three major types of algorithms: <\itemize> Cyclic extension methods (section); Geometric progression methods (section); FFT based methods (section). For the first two types of methods we mostly review existing algorithms, although we do propose some new variants and optimizations. The third FFT based method is new, as far as we are aware. For each of the three methods, an important is to evaluate > modulo -1> for one or more suitable orders , after which we reconstruct > from these modular projections. Cyclic extension methods directly evaluate over the cyclic extension ring /-1|)>>>. This has the advantage that can be freely chosen in a suitable range. However, the evaluation of over such a large cyclic extension induces a non-trivial overhead in the dependency of the complexity on . Geometric progression methods rather evaluate at a sequence ,\,\> of pairwise distinct elements in (or inside an extension of of modest degree ). If > is a finite field, then > necessarily has an order that divides (or -1> when working over an extension of degree ). Although the evaluations of become more efficient using this approach, the recovery of > modulo -1> from ,f|)>,\,f|)>> requires extra work. The cost of this extra work highly depends on the kind of orders that can be taken as divisors of (or -1> for small ). Theoretically speaking, the existence of suitable orders is a delicate issue; in practice, they always tend to exist as long as >>; see sections, and for some empirical evidence. Geometric progression methods allow us to take much larger than , but they involve anon-trivial cost for recovering > modulo -1> from ,f|)>,\,f|)>>. If |)>>>, then this cost may even dominate the cost of the evaluations of . In such situations, an alternative approach is to evaluate at ,\,\> and to recover > modulo-1> using one inverse DFT of length . However, this puts an even larger strain on the choice of , since it is important that r=O> for this approach to be efficient. Moreover, the recovery of > from its reductions modulo -1> for several orders of this type is more delicate and based on a probabilistic analysis. Yet, suitable orders again tend to exist in practice as long as >>. The expected complexities of the best versions of the three approaches are summarized in Table. These bounds rely on two types of heuristics: <\description-compact> >For d>, the exponents of are essentially randomly distributed modulo . >For >, the number -1> admits a \Plarge smooth divisor\Q. We will present a more precise version of> in section. The heuristic> will be made precise in sections and and numeric evidence is provided in sections, and. The last section is dedicated to the interpolation of multivariate polynomials. We start with the well known strategies based on Kronecker segmentation (section) and prime factorization(section). For sparse polynomials in many variables, but with amodest total degree , we also recall the inductive approach by packets of coordinates in section. If q>, then geometric progression and FFT based methods should be particularly favorable in combination with this inductive approach, since one can often avoid working over extensions of in this case. We conclude section with a few algorithms for special situations.|||||>||>|)>>>|> and >>>||>|)>*T*|)>*log|)>>>|>, >, and >>>>||>*T*|)>*log|)>>>|>, >, and >>>>>>>> Heuristic complexity bounds for the three main approaches to sparse interpolation. > We wish to thank the anonymous referee for helpful corrections.> One remarkable feature of the finite field > with elements is that every \> satisfies the equation =a>. In particular, for any sparse polynomial as in() and with coefficients in >, the polynomial takes the same values as <\equation*> ,\,i>f,\,i>*x rem >*\*x rem >, for ,\,x\\>, where \P\Q stands for the remainder of a Euclidean division. In other words, the exponents of are only determined modulo , so we may assume without loss of generality that they all lie in ,q-2|}>> and that the total degree of satisfies n*>. On the other hand, in the case where our black box function can be evaluated not only over >, but also over field extensions >> (this typically holds if is given by an expression or a directed acyclic graph (dag)), then the exponents in the expression() can be general non-negative integers, but the above remark shows that we will crucially need to evaluate over extensions fields >> with 1> in order to recover exponents that exceed . More generally, if we choose to evaluate only at points ,\,a|)>\\>> such that ,\,a>> are -th roots of unity, then we can only hope to determine the exponents modulo in the expansion(). In that case, must divide the order -1> of the multiplicative group of >>. In addition, as recalled in sections and below, several important tools such as polynomial root finding and discrete logarithms admit faster implementations if we can take of the form *r> with =O> and where > is smooth. Sometimes, primitive roots of unity of such orders already exist in >. If not, then we need to search them in extension fields >> with 1> as small as possible. Let us briefly investigate the prime factorization of -1> for various prime numbers and small . As observed in4.1 and4.6>, the number -1> typically admits many small prime divisors when is itself smooth. This phenomenon is illustrated in Table for small values of . Given , it seems easy in practice to find a small and divisors *r\-1|)>>> such that \T> and > is smooth.|||>|-1>>|-1>>|-1>>>|>||>|>>>|>|>|>>|\3>>>|>|>|13>>|\31>>>|>|5>>|\5>>|\3\13>>>|>|\7>>|\7\13>>|\3\7\31>>>|>|5\17>>|\5\41>>|\3\13\313>>>|>|\5\7\13>>|\5\7\13\73>>|\3\7\13\31\601>>>|>|5\17\257>>|\5\17\41\193>>|\3\13\17\313\11489>>>|>|\5\7\13\17\241>>|\5\7\13\41\73\6481>>|\3\7\13\31\313\601\390001>>>||\7\11\31\151\331>>|\7\11\13\31\61\271\4561>>|\3\7\11\31\61\71\181\521*\>>>|>|\5\7\13\19\37\73\109>>|\5\7\13\19\37\73\757\530713>>|\3\7\13\19\31\37\601\829*\>>>>>>>> Prime factorizations of -1>, -1>, and -1> for various small smooth values of . > For larger , we may need to resort to larger extension degrees in order to find appropriate orders . A natural question is whether -1> is guaranteed to have a non-trivial smooth divisor for large and a fixed value of . This leads us to introduce the following guaranteed lower bound: <\equation> a\lim\\> gcd-1\q\q|)>, where -1\q\q|)>> denotes the of the integers -1> where runs over the prime numbers greater than or equal to>. In Table, we listed the prime factorizations of -1> and > for and various small smooth values of. Here was chosen such that /2> is also prime. For the practical applications in this paper, the table suggests that it remains likely that suitable orders can still be found whenever needed, and that> is usually quite pessimistic, even for large values of . Let us finally mention that the sequence> coincides with Sloane's integer sequence ; see >.|||>|-1>>|>>>|>|649871>>|>>|>|\3\4513\649871>>|\3>>>|>|7\13\397\649871\6680137>>|>>|>|\3\5\193\4349\4513\40253\649871>>|\3\5>>>|>|\3\7\13\31\397\4513\649871\6680137\18164844799>>|\3\7>>>|>|\3\5\17\73\193\241\4349\4513\40253\649871\298090889\941485217>>|\3\5>>>|>|\3\7\13\31\193\397\4349\4513\40253\649871\6680137\387205657*\>>|\3\5\7\13>>>|>|\3\5\17\73\193\241\4349\4513\40253\649871\3955153\298090889*\>>|\3\5\17>>>|>|\3\5\7\13\17\31\73\193\241\397\4349\4513\40253\649871*\>>|\3\5\7\13>>>|>|\3\7\11\13\31\61\71\271\397\701\881\1171\2411\4513\649871*\*>>|\3\7\11\31>>>|>|\3\5\7\13\19\31\37\109\193\397\757\1657\4349\4513\40253\649871*\>>|\3\5\7\13\19\37>>>>>>>> Prime factorizations of -1> and > for and various values of . > As already mentioned in the introduction, most algorithms for sparse interpolation are probabilistic of Monte Carlo type. We notice that it is easy to check (with high probability) whether a candidate sparse interpolation >> of is correct. Indeed, it suffices to evaluate >> at random sample points and check whether the result vanishes. Deterministic algorithms exist but with higher complexities; see for instance. Many algorithms for sparse interpolation require extra information, such as boundst> and d> for the number of terms and the total degree of. Furthermore, several algorithms are only able to guess some of the terms of with high probability, but not all of them. In this section, using ideas from, we show how to turn such \Ppartial\Q algorithms for sparse interpolation into full-blown ones. Provided that the characteristic of is zero or sufficiently large, we also show how to upgrade an interpolation method modulo -1,\,x-1|)>> into a general algorithm, following. Assume that we have an algorithm for \Papproximate\Q (or \Ppartial\Q) sparse interpolation that takes ablack box for as input, together with bounds and for and . The algorithm should always return a sparse polynomial >>, of total degree at most and with at most terms. Moreover, for some fixed constant \\1>, if T> and D>, then >>> should contain at most *t> terms, with high probability. If T> or D>, then we allow>> to be essentially arbitrary under the above constraint that >> has at most terms of degreeD>. Then we may use the following strategy for \Pfull\Q sparse interpolation: <\quote-env> <\specified-algorithm> a polynomial black box function ,\,x|)>>. the sparse interpolation >> of . <|specified-algorithm> <\enumerate> Let >\0> be an initial approximation of . Set initial bounds 1> and 1> for the number of terms and total degree. Determine the approximate sparse interpolation >> to \f-f>> using and as bounds for the number of terms and the total degree of >. Set >\f>+\>>. If >> with high probability, then return >>. Whenever appropriate, increase and/or , and reset >\0>. Return to step3. In step 6, the precise policy for increasing and may depend on the application. We typically double when is suspected to exceed and we double when is suspected to exceed . In this way, the bounds and are at most twice as large as the actual values and , and the running time is essentially a constant times the running time of the approximate sparse interpolation with bounds and . However, for this \Pmeta complexity bound\Q to hold, it is crucial in step3 that the sparse approximation >> can be evaluated efficiently at the sample points used during the sparse interpolation (the naive evaluation of a polynomial with terms at > points would take time |)>>, which is too much). Fortunately, this is the case for all sparse interpolation strategies that will be considered in this paper. When do we suspect or to be too low in step6? In the case of , a natural idea is to test whether the number of terms in >> or >> exceeds a fixed constant portion of . This strategy assumes that >> be essentially random when is too small (if the number of terms of >> is much smaller than whenever T>, then we might need more than > iterations before the algorithm converges). The handling of exponents and degree bounds is more involved. When interpolating over a finite field, all non-zero evaluation points are roots of unity, so the exponents can only be determined modulo a certain integer (or even modulo a submodule of>). If the characteristic of is sufficiently large, then the exponents can be computed directly: see the next subsection. Otherwise, we need to recombine reductions with respect to several moduli: see section below. This also provides a natural test in order to check whether D>. Indeed, it suffices to compute the sparse interpolation of for one or more extra moduli and check whether the results agree with our candidate interpolation. Assume that we have an algorithm that allows us to compute the sparse interpolation of modulo =-1,\,x-1|)>> for given moduli . Assume also that we have access to the program that computes , typically in terms of a straight-line program. If or max > f,\,deg> f|)>>, then let us show how to derive an algorithm for the sparse evaluation of . It is well known that the technique of \Pautomatic differentiation\Q (due to Baur and Strassen) allows us to evaluate the gradient f/\ x,\,\ f/\ x|)>> using at most operations in . Using operations, this provides an algorithm for the simultaneously evaluation of ,\,g> with =* f/\ x|)>>> for ,n>. With a small overhead, this next allows us to jointly compute the sparse interpolations of ,\,g> modulo >. Now assume that >*\*x>> is a term of such that for any other term *x>*\*x>> of , we have mod r,\,e mod r|)>\ mod r,\,e mod r|)>>; we say that this term \Pdoes not collide\Q modulo . Then the sparse interpolation of modulo> contains the non-zero term mod r>*\*x mod r>>. Moreover, given ,n|}>> with \0>, the sparse interpolation of > modulo > also contains the non-zero term *c*x mod r>*\*x mod r>>. This allows us to retrieve > through one simple division *c|)>/c>. Furthermore, if the modulus was picked at random with t>, then there is a high probability that a fixed non-zero proportion of terms in do not collide modulo . Combined with Algorithm, this yields an algorithm for obtaining the sparse interpolation of. This strategy for sparse interpolation was first exploited by Huang. <\remark*> For simplicity, we consider sparse interpolation of polynomials over fields in this paper. In fact, the algorithms also work for vectors of polynomials in ,\,x|]>>>>, by considering them as polynomials with coefficients in >>. We implicitly used this fact above when saying that we \Pjointly\Q compute the sparse interpolation of ,\,g>> modulo >. In summary, we have shown how to reduce the general problem of sparse interpolation to the case when <\enumerate> we have bounds for the number of terms and the total degree, and we only require an approximate sparse interpolation (in the sense of section). One general approach for sparse interpolation of univariate polynomials over general base fields was initiated by Garg and Schost. It assumes that the black box function can be evaluated over any cyclic extension of the form /-1|)>>. The evaluationof <\equation> f=c*x>+\+c*x> at inside such an extension simply yields <\equation*> f>\c*x rem r>+\+c*x rem r>. In absence of \Pcollisions\Q =e> modulo for j>, this both yields the coefficients of and its exponents modulo . By combining the evaluations for various moduli, it is possible to reconstruct the actual exponents using Chinese remaindering. Throughout this section, we assume that we are given bounds and for the degree and the number of terms of. Garg and Schost's original algorithm was deterministic under these assumptions. However, their algorithm was not designed to be efficient in practice. In the past decade, many variants have been proposed. Roughly speaking, they all follow the same recipe that is summarized in Algorithm below. The variants mainly differ in the precise way recombinations are done in step3. <\quote-env> <\specified-algorithm> a black box polynomial >, a degree bound , a sparsity bound . an approximate sparse interpolation >> of . <|specified-algorithm> <\enumerate> Determine suitable moduli ,\,r\T> with ,\,r|)>\D>. Evaluate at in />-1|)>> for ,l>. Determine matching terms in the expansions of |]>>,\,f|]>>> that are likely to be the reductions of the same term *x>> in the expansion of . Return the sum >> of all terms *x>> as above. For all matching strategies that have been considered so far, the cost of steps1 and3 is dominated by the cost of step2. If the evaluation of only involves ring operations, then the running time of Algorithm is therefore bounded by |)>+\+|)>|)>*L|)>>. The moduli> are usually all of the same order of magnitude \T>> for some small \1> that depends on the matching strategy. Then we may take >, so the cost simplifies to >|)>*log D/log T|)>>. For finite fields >>, this cost becomes >>*log D*log q|)>>. For the design of matching strategies, it is therefore important that we can take > as small as possible. <\remark*> The above analysis can be refined by maintaining separate counts >, >, and> for the numbers of additions (or subtraction), multiplications, and divisions that are necessary for one evaluation of . Then the cost of Algorithm over> becomes >>*log D*log q*+L*log T+L*|)>|)>>. <\remark*> The complexity analysis may need to be adjusted somewhat if is so large that we run out of suitable moduli >. If our matching strategy requires prime numbers of the order of >, then this happens when exceeds approximately the same order>. In that case, we need to replace by an appropriate power of in our complexity bounds. Alternatively, if the characteristic of is sufficiently large, then we may fall back on the strategy from section. Garg and Schost's original algorithm from uses prime numbers for the moduli>. Assuming that >> admits terms, the algorithm is based on the observation that the projection of the polynomial |)>*\*|)>> modulo coincides with rem p|)>>*\* mod p|)>>>. This allows for the recovery of |)>*\*|)>> through Chinese remaindering, by working with a sufficient number of primes. It then suffices to determine the zeros ,\,e>> of and to recover > as the coefficient of rem p>> in >> for ,t>. However, this strategy is very sensitive to collisions, and requires T> in order to ensure with high probability that >> admits exactly terms. In other words, it forces us to take \2> in the complexity analysis. Garg and Schost's original algorithm is actually deterministic and uses (L*T* )> operations in. The derandomization is achieved by using *log D|)>> different primes. Another matching strategy for step3 of Algorithm has been proposed by Arnold, Giesbrecht, and Roche. The idea is to pick =p*p> for ,l>, where ,\,p>> are primes with \T> and ,\,p\T>> for some \0> (so that =1+\>). For ,t> and ,l>>, there then exists a fixed non-zero probability such that the term *x rem p>> of|]>>> matches a term *x rem r>> of |]>>>. Let > be the set of indices for which we have a match. For some fixed constant 0>, we then have \k\\|)>\T*l>> with high probability. By taking log D/*log T|)>> in step 1, this implies \k\\|)>\D>. With high probability, this allows us to reconstruct those terms *x>> such that \e>> modulo> for all i>. The sum of these terms gives the desired approximation >> of for which a fixed nonzero proportion of terms are likely to be correct. Giesbrecht and Roche proposed yet another matching strategy which is based on the concept of \Pdiversification\Q. The polynomial is said to be if its coefficients> are pairwise distinct. Assuming that is sufficiently large, it is shown in that the polynomial *x|)>> is diversified with high probability for a random choice of \K>>. Without loss of generality, this allows us to assume that is diversified. In step3 of Algorithm, we then match a term > of |]>>> with a term *x>> of >|]>>> if and only if >>. Giesbrecht and Roche's original algorithm uses > moduli ,\,r> of size *log D|)>>. Consequently, their algorithm for sparse interpolation uses *|)>> operations in >. As we will see below, their probabilistic analysis is actually quite pessimistic: in practice, it is possible to take ,\,r=> as long as >>. Let us now focus on the design and analysis of a probabilistic algorithm that exploits the idea of diversification even more than Giesbrecht and Roche's original method from. Given a diversified polynomial , together with bounds and for its degree and its number of terms, our aim is to compute > terms of , with high probability. Our algorithm uses the following parameters: <\itemize> A constant \1>. Orders \\\r> that are pairwise coprime, with \\*T> for ,l>. The minimal number \,l|}>> such that ,\,r>|)>=r*\*r>>\D>. The precise choice of > and will be detailed below; the parameter > and the ratio>> should be sufficiently small for the algorithm to be efficient, but > and > should also be sufficiently large for our algorithm to succeed with high probability. We now use Algorithm below in order to compute an approximate sparse interpolation of. It is a variant of Algorithm with a matching strategy that has been detailed in steps2, 3, and4. Each individual term > of is reconstructed from only asubset of its reductions modulo >-1,\,x>-1>. <\quote-env> <\specified-algorithm> a diversified black box polynomial > and parameters as above. an approximate sparse interpolation >> of . <|specified-algorithm> <\enumerate> Compute |]>>=f rem >-1|)>> for ,l>. Let >\0>. Let > be the set of all coefficients that occur once in |]>>>, for ,l>. For each C\\\C> do: <\indent> If \c\C|}>> is such that \k\\|)>\D>, then: <\indent> Determine the unique exponent D> such that >> occurs in|]>>> for every \>, and set >\f>+c*x>. Return >>. How to ensure that a non-zero portion of the terms of >> can be expected to be correct? In order to answer this question, we make the following heuristic hypothesis: <\description-compact> >>>For ,l>, the modular reductions of exponents rem r> for ,t> are uniformly distributed in ,r-1|}>>. The distribution associated to > is independent of the one associated to >> whenever k>. Such a heuristic is customary in computer science, typically when using hash tables. According to , the probability that a fixed term > does not collide with another term modulo > is <\equation*> |)>\*T|)>|)>\*T|)>|)>\\*T|)>|)>>\\>. Setting \1-\>>, this probability tends to > for large . The probability that > collides with another term modulo > for exactly > values of is therefore bounded by <\equation*> >*\>*|)>> and this bound is sharp for large . Consequently, the probability that we cannot recover aterm >*x>> in step4 from its reductions modulo >-1> for ,l>> is at most <\equation> P;l,1-\|)>=\>*\*|)>, and this bound is sharp for large . The probability() has been well studied; see for a survey. Whenever <\equation*> a\\/l\1-\, Chernoff's inequality1> gives us <\equation*> P;l,1-\|)>\exp>|)>+*log>|)>|)>|)>. Let \1> be a positive real parameter. In order to ensure ;l,1-\|)>\\> it suffices to have <\equation*> l*>|)>+*log>|)>|)>\log|)>. Now thanks to2(a)> we have <\equation*> a*log>|)>+*log>|)>\2*|)>|)>, so it suffices to ensure that <\equation> -l*|)>|)>\|)>|2>*l. Now let us take =c*|)>*l> with 1>, after which is equivalent to <\equation> \|)>|2*l*|)>>. For fixed > and large ( for large ), it follows that we may take arbitrarily close to. Summarizing, we may thus take \/|)>> in order to recover an average of at least |)>*t>> correct terms, where > can be taken arbitrarily close to : <\proposition> Assume > and let >, >, , ,\,r>, and > be parameters as above. Given \1>, assume that =c*|)>*l>, where satisfies)>. Then Algorithm returns at least |)>*t>> correct terms of on average. Let us now analyze the running time of Algorithm. Taking \/|)>>, and since *T|)>>\D>, the cost of step1 is <\equation*> \\T*\*|1-\>=T*\*\*\>=O*\*\>|)>, and *\>> reaches its minimum value > at =1>. This means that the total complexity is best when > is close to . In other words, this prompts us to take =1>, \\\r\T>, \log D/log T>, and \/|)>>. For this choice of parameters, we obtain the following heuristic complexity bound: <\proposition> Assume that and >. Given \\1> and a diversified polynomial of degree D> and with T> terms, there exists a Monte Carlo probabilistic algorithm which computes at least |)>*t> terms of in time <\equation*> O>. <\proof> We take > to be the -th smallest prime numbers that is larger than , so that =O> is the smallest number with *\*r>\D>. We also take =1>, =>>, and let > be smallest such that() is satisfied for /|)>|)>>. Combining and, we may compute ,\,r> in time >>. Now the running time of step1 of Algorithm is <\equation*> O|)>+\+|)>|)>*L|)>. With >, this cost simplifies to <\equation*> O>*log D/log T|)>=O>. Step3 may be obtained by sorting the coefficients of the |]>>>, in time <\equation*> O=O. Using fast Chinese remaindering, step4 takes time *log D|)>=O>>. <\remark*> If T>, then *x|)>> is diversified with high probability for a random choice of \\\>: see. In the range where > and |)>>, it is possible to work with a slightly weaker assumption: we say that is if ,\,c|}>> is of size>. If >, then the polynomial *x|)>> is still weakly diversified, for a random choice of \\\>. If is only weakly diversified and \i\t\j\i,c\c>|}>>, then our analysis can be adapted to show that Algorithm returns about |)>*t>> correct terms of on average. Finally, in the range where >, we need to work over a field extension >> with \T>, which leads to an additional arithmetic overhead of>>>. <\remark*> Let us show that with high probability, the polynomial >> returned by Algorithm only contains correct terms of (although it does not necessarily contain all terms). For this, we make the additional hypothesis that the coefficients of are essentially random non-zero values in > (which is typically the case after a change of variables \f*x|)>>, where \\\> is random). Now assume that some coefficient in step 4 gives rise to a term > that is not in . Then for every \>, there should be at least two terms in that collide modulo > and for which the sum of the corresponding coefficients equals . The probability that this happens for a fixed \> is bounded by > and the probability that this happens for all \> is bounded by >>, where \\> is minimal with +1>*\*r\D>. For the algorithms in this section, we assumed that a bound for was given. It turns out that a variant of our probabilistic analysis can also be used for the efficient computation of a rough estimate for . This yields an interesting alternative to the naive doubling strategy described in section. Let us still assume that holds. We will also assume that colliding terms rarely cancel out (which holds with high probability if is sufficiently large). This time, we compute rem -1|)>> for \*t>, where \1> is to be determined, and let be the number of terms in this remainder. When randomly distributingballs overboxes, the probability that none of the balls lands in a particular box is >. Consequently, the expected number of boxes with no ball is *B>, whence <\equation*> B-N\|)>*B\\>>*B. It follows that <\equation*> >\log|)>, and thus <\equation> t=>\B*log|)>. By doubling until >, we may then use the approximation() as a good candidate for . Notice that we have > when 2*t/log t>. Cyclic extension methods for sparse interpolation are attractive due to their generality and the possibility to derive deterministic complexity bounds. However, even their most efficient probabilistic versions suffer from the overhead of arithmetic operations in cyclic extension algebras /-1|)>>. The matching strategy based on diversification leads to the best practical complexity bounds, as shown in section. Assuming, >, and >, we have given a Monte Carlo algorithm for sparse interpolation of complexity >.> The case when > can be reduced to this case using a field extension of degree >. Assuming only and >, we thus obtain a probabilistic algorithm that runs in time <\equation> O>|)>. Prony's method is one of the oldest and most celebrated algorithms for sparse interpolation of univariate polynomials. It is based on the evaluation of at points in a geometric progression. Since there are many variants, we keep our presentation as general as possible. As in the previous section, assume that <\equation> f=c*x>+\+c*x> and that we know bounds and for the degree and the number of terms of . <\quote-env> <\specified-algorithm> a black box polynomial >, a degree bound , a sparsity bound . the sparse interpolation >> of . <|specified-algorithm> <\enumerate> Find a suitable element \K> with multiplicative order max>>, while replacing by an extension if necessary. Evaluate ,f|)>,f|)>,\,f|)>>. Compute a minimal T>, a monic \K> of degree , and an K> of degreet>>, such that the following identity holds modulo |)>>: <\equation> \>f|)>*z=i\t>|1-\>*z>=|\>. \; Find the roots >> of \z+\*z+\+\\K>. Compute the discrete logarithms of the roots to base > to discover the exponents ,\,e> of as in(). Compute ,\,c> from ,\,f|)>> and ,\,e>, using linear algebra. It is well known that steps3 and6 can be performed in time *log T|)>>, through fast Padé approximation in the case of step3, and using atransposed version of fast multipoint interpolation for step6. If >, then this bound becomes >*log q|)>>. The efficiency of steps4 and5 highly depends on the coefficient field. In the remainder of this section, we will discuss this issue in detail in the case when > is a finite field. Finding the roots of a univariate polynomial over a finite field is a well-known and highly studied problem in computer algebra. The most efficient general purpose algorithm for this task is due to Cantor and Zassenhaus. It is probabilistic and computes the roots of> in time >*|)>>. In, several alternative methods were designed for the case when *r> with =O> and \T>> smooth (in the sense that =O>> for each prime factor> of>). The idea is to proceed in three steps: <\enumerate> We first compute the >-th Graeffe transform >|)>> of >, whose roots are the >-th powers of the roots of >. This step can be done in time >*log q|)>> by5>. We next compute the roots of >|)>> through an exhaustive evaluation at all >th roots of unity. This step takes time >>. We finally lift these roots back up to the roots of >. This can be done in time >*log q*log r|)>=O>*log q|)>>, using g.c.d. computations. Altogether, this yields a sparse interpolation method of cost >*log q|)>>. The back-lifting of single roots can be accelerated using so-called \Ptangent Graeffe transforms\Q. The idea is to work over the ring K|]>/|)>> instead of. Then \K> is a root of apolynomial \K> if and only if +\\R> is a root of the polynomial |)>\R>>. Now if we know a single root +\|)>>=\>+r*\-1>*\> of >|)>|)>>, then we may retrieve> using one division of >>by*\-1>> and one multiplication by > (note that > is invertible in > since >divides>). In other words, the back-lifting step can be done in time >>, using> operations in>. However, this method only works if >> is a single root of >|)>>. When replacing > by|)>>> for a randomly chosen \\>, the polynomial >|)>|)>> can be forced to admit anontrivial proportion of single roots with high probability. However, these roots are no longer powers of >, unless we took . Assuming that and using several shifts >, it can be shown12> that the tangent Graeffe method yields a sparse interpolation method of complexity >|)>>. The discrete logarithm problem in abelian groups is a well-known problem in computational number theory. If is smooth, then the Pohlig\UHellman algorithm provides an efficient solution; it allows step5 of Algorithm to be performed in time >>. Under the assumption that we may take >>, this bound reduces to >>. Again, the same bound still holds if *r> with =O> and > smooth. Indeed, in that case, we may tabulate the powers ,\>,\>,\,\-1|)>*r>>. This allows us to efficiently determine the discrete logarithms of *r>,\,\*r>> with respect to >>, which yields the exponents ,\,e> modulo >. We next use the Pohlig\UHellman algorithm to compute,\,e>>. If > exceeds (or if admits no suitable factors that allows us to apply the above methods), then we may need to work over an extension field >> of >. Notice that this requires our black box representation of to accept inputs in such extension fields. Since evaluations over >> are at least times more expensive, it is important to keep as small as possible. If q>, then we must necessarily have \2*T+1>, whence log/log q>. In general, we want to take |log T/log q|\>|)>>. Since we also needD>> in step 1, this leads to the constraint >>. Under this hypothesis and using the observations from section, it is likely that a suitable extension order and divisor -1|)>> can indeed be found. Still denoting by > the cost of multiplication of polynomials of degree over >, the total cost of sparse interpolation then becomes <\equation> O>|)>*T*|)>. An interesting question is whether we can reduce the number of evaluation points when working over an extension field >>. Indeed, if > is the Frobenius map of >> over >, then |)>=\|)>> for all \>>. If >=\>, then this allows us to obtain the evaluations at the distinct points ,\,\> using a single evaluation at . In step2 of Algorithm, we can therefore avoid the evaluations at > for > and gain aconstant factor > for large . Similarly, we may compute all values ,f|)>,\,f|)>> using approximately evaluations of only; whenever > is small, this allows us to gain a factor . It would be interesting to know whether it is possible to do better and regain a factor > in general. Besides working in an extension field, it may also be interesting to perform part of the computations over a subfield of >. Indeed, the main aim of steps4 and5 of Algorithm is to find the exponents of . Now assume that > admits a subfield >> with \q> and let \\>> be the corresponding trace function. Then and are likely to admit approximately the same exponents. Taking \\>> in step1, we may thus replace ,\,f|)>> by their traces after step2, and determine the exponents of instead of . Although this does not allow us to speed up steps2 and6, we do gain a factor of at least > in steps4 and5. > Once the order has been fixed, Algorithm essentially allows us to interpolate > modulo -1>. With respect to the cyclic extension approach, the main difference is that one expensive evaluation of over the extension ring /-1|)>> is replaced by cheap evaluations over plus >|)>> scalar operations. If T>, then we may also evaluate > modulo -1> for different moduli ,\,r>> and recombine the results using one of the matching strategies from section. However, in the present context, we are not completely free to pick our moduli, since we need corresponding primitive roots of unity \\>>,\,\\>>>> of orders ,\,r> in small extensions >>,\,\>>> of >. Let us focus more specifically on Algorithm, which requires in particular that T>. We need ,\,r\2*T> to be suitable for steps4 and5, so =r*r> with =O> and > is smooth. On the other hand, we may relax the conditions on ,\,r>. In this case, the complexity does not depend on >, so it is better to choose the > much larger than, preferably of the order of >. It is also not necessary that ,\,r> be pairwise coprime: it suffices that >,\,r>>|)>\D> for any k\\\k>\l>. Ideally speaking, we should have \\log D/log T>. Although there is no reason for suitable ,\,r> of this kind to exist, we expect that this will often be the case in practice, as long as >. Evidence in this direction will be presented in sections and below, under even stronger constraints on the orders ,\,r>. Assuming that we are indeed able to find suitable ,\,r>, the overall runtime complexity becomes <\equation*> O>|)>*T*|)>+\+|)>|)>|)>. When using naive arithmetic with =O>*log q|)>> and assuming that =O>, this complexity bound simplifies into <\equation*> O>|)>*T*|)>*log q|)>. In summary, the efficiency of the geometric progression approach over a finite field rests on our ability to find suitable divisors of -1> for small values of . If T> and q>>, then we essentially need an order -1|)>> of the type *r> with \T>, >smooth, and small. By what has been said in section, it is very likely that such and always exist. If >> and q>, then it remains likely that various divisors of this type can be combined, as explained in section. If T>, then we first need to replace by asuitable extension. Assuming that >> and that suitable orders as above can indeed be found, the cost of the geometric progression approach is bounded by <\equation> O>|)>*T*|)>*log |)>. In favorable cases when q=T>> and is smooth, we obtain the complexity bound <\equation> O>*T*log q|)>, instead of(), by using the tangent Graeffe method. In comparison with algorithms based on cyclic extensions, the main advantage of the algorithms in this section is that we avoid expensive evaluations over cyclic extension rings. On the other hand, the cost of the root finding step may dominate the cost of the evaluations of whenever |)>>; in that case, cyclic extension methods may become competitive. Methods based on geometric progressions are also less suited for the super sparse case. <\remark> For practical implementations, it is not convenient to use an comparison in order to determine which of the or > terms dominates the cost in(). Since we usually try to interpolate for increasingly large bounds, it is better to test whether step1 of Algorithm requires more running time than step3 as we increase. Whenever the cost of step3 starts to dominate, we may switch to acyclic extension style approach (or an FFT based approach, to be presented next). Geometric progression style algorithms for sparse interpolation admit the big advantage that they have a sharp complexity with respect to . However, if evaluations of are cheap, then the> term in() may dominate . In this section, we will investigate a strategy to reduce the dependency in , at the price of > more evaluations. The idea is to more aggressively exploit the observations from sections and. Alternatively, we can regard our proposal as a careful adaptation of the cyclic extension approach that allows us to replace evaluations over cyclic extensions by evaluations over itself.\ -1>> For a fixed parameter > that we control and a modulus close to *T>, let us first study how to evaluate > efficiently modulo -1>>. Instead of evaluating at only points ,\,\>> for a primitive -th root of unity > (as in step 2 of Algorithm), we rather evaluate at all -th roots of unity ,\,\>. The advantage is that we may then use FFT-based techniques as a replacement for the remaining steps of Algorithm. Moreover, if > and > lives in an extension =\>> of , then it suffices to evaluate at only r/s> points in order to determine all values ,\,f|)>> using the Frobeniusmap. Recovering > modulo -1> from these values can also be done efficiently using the inverse Frobenius FFT. <\quote-env> <\specified-algorithm> a black box polynomial > over > and \\>> of order \>. a polynomial >\\> of degree r> with >> modulo -1>. <|specified-algorithm> <\enumerate> Let > be a section of ,\,\|}>> under the Frobenius map > that is suitable for the Frobenius DFT. Compute |)>> for each \\>. Compute the coefficients of >> by applying one inverse Frobenius DFT of order to |)>|)>\\>> and return >>. Assuming for simplicity that |\|>=\> and that the computation of an inverse Frobenius DFT is > times faster than the computation of a full DFT, we note that the cost of one run of Algorithm is bounded by <\equation> O>*T*\*|s>|)>. If T> and T>, then we wish to apply a similar strategy as in section and recombine the values of rem -1|)>> for various moduli ,\,r>. This time, in order to perform step1 of Algorithm using Algorithm, we need the orders ,\,r> to be close to . With\1> as in section, we thus assume =\*T,\,r=\*T> with ,\,\>\\\1>. As in section, we also impose the condition that >,\,r>>|)>\D> for any k\\\k>\l>. Ideally speaking, we have \log D/log T> and \/>|)>>. Under these conditions, thanks to, the overall running time of Algorithm is bounded by <\equation> O>*T**|)>|s>+\+\*|)>|s>|)>|)>. When using naive arithmetic with =O>*log q|)>> and assuming that =O>, this bound simplifies into <\equation*> O>*T*|)>*log q|)>. Due to our hypothesis that \/>|)>> and the analysis from section, we can still expect the algorithm to return about half of the terms of . <\remark> If T>, then we need to replace by an extension of degree at least before being able to diversify , without being able to profit from the Frobenius map. This leads to an incompressible overhead of . Nevertheless, in the special case when the exponents ,\,e> are already known, the matching step4 of Algorithm can be simplified, and working over an extension can be avoided. >> In order to apply Algorithm as in the previous section, we need primitive roots of unity of suitable orders ,\,r>> in algebraic extensions of>. Although we have no general theory, it is instructive to study a few practical examples in order to gain insight what we may reasonably hope for. In this section, we consider the case when , >, and>. First of all, we need to replace > by a sufficiently large extension >> in order to diversify (at least in the weak sense from the end of section). Since it is favorable to take <\equation*> |log q>\|log 2>\20 as smooth as possible, we opt for =q>. For >>, we next take our orders \10> with \>-1|)>> and > as small as possible, and such that <\equation*> \\lcm,\,r|)> is as large as possible: <\equation*> |||||||||||||||||||||||||||||>|||>|>||||151\331>|>|>|>|10>>|>||||>||||61\1321>||>|>|10>>|>||||>||||23311>||>|>|10>>|>||||>||||\11\19\631>||>|>|10>>|>||||>||||61681>||>|>|10>>|>||||>||||\7\41\241>||>|>|10>>|>||||>||||4051>||>|>|10>>|>||||>||||1801>||>|>|10>>|>||||>||||100801>||>|>|10>>|>||||>||||7\54001>||>|>|10>>>>> Taking =1>, we obtain \18/6=3> and \*\>. To be on the safe side, we take . The minimum least common multiple of three ( four) orders among ,\,r> is <\equation*> lcm,r,r|)>\4.5\10 ( ,r,r,r|)>\1.1\10>), so we have =4>. Notice that some of the =r/T> are quite a bit larger than >, with an average of =+\+\|)>/10\1.26>. For and *l\4.52\\> , we therefore consider it likely that the condition \k\\|)>\D> in step 4 of Algorithm is satisfied with probability 1/2> (a rigorous analysis would be welcome). >> As our second example, consider the case when (as in section), >, and >. This time > contains at least elements, so we may directly work over >. Proceeding as in the previous subsection, we may now take: <\equation*> ||||||||||||||||||||||||||||||>|||>|>||||649871>|>|>|>|10>>|>||||>||||\3\4513>||>|>|10>>|>||||>||||\40253>||>|>|10>>|>||||>||||193\4349>||>|>|10>>|>||||>||||2411>||>|>|10>>|>||||>||||13\31\397>||>|>|10>>|>||||>||||\17\73\241>||>|>|10>>|>||||>||||3\7\37\109>||>|>|10>>|>||||>||||11\71\881>||>|>|10>>|>||||>||||148853>||>|>|10>>>>> With =1>, , and =4> as before, the minimum least common multiple of > orders among ,\,r> is ,r,r,r|)>\1.2\10>. Again, the \Peffective\Q mean value =+\+\|)>/10\1.5> satisfies >*l\5\\>. Notice that > approximately grows like \1.1\k>. On the first example from the previous subsection, we rather have \0.7\k>. This suggests that the growth of > might often be somewhat above*k>. As an optimization, we also notice that it should be better to take =r/T>> somewhat larger for small values of , so as to balance the terms *|)>/s> in the sum(). For instance, the last number > is quite a lot larger than ; it should be more efficient to resort to a larger extension degree =13> and take =1220444=2\305111>. Another option would be to include =3> and =6680137>. In addition, by increasing > and >, it might also be possible to take . As mentioned in the introduction, one important application of sparse interpolation over finite fields is sparse interpolation over the rational numbers. In that case, we typically combine sparse interpolation over different prime fields > for which is divisible by a high power of two. We usually proceed in two stages. We first determine the exponents ,\,e> of (typically using a few primes only). Once the exponents are known, the computation of the coefficients modulo further primes (typically many ones) becomes more efficient. We finally use the Chinese remainder theorem and rational number reconstruction to determine the actual coefficients. Let> denote the Euler totient function. For any coprime positive integers and we we write > for the number of prime integers x> such that . Roughly speaking we have <\equation> \\*log x> for sufficiently large values of ; recent advances in this research area can be found in. Unfortunately, for practical purpose, the bounds presented in require to be exponentially larger than . Assuming the Generalized Riemann Hypothesis (GRH), known estimates of > are more favorable: in fact, the explicit bound proved in raises this requirement. The techniques from this section may help to accelerate both stages. If is normally sparse, then we may directly apply the algorithms from above to find the exponents. If is super sparse, then we will now present a multi-modular variant of the strategy from section. Given a number T>, the aim is to determine with high probability those exponents> such that \e> for all i> and use this to reconstruct an approximate sparse interpolation of . We next apply Algorithm. So let T> be fixed and consider successive prime numbers ,q,\> for which -1> is divisible by . Let be minimal such that *\*q\D>. Under GRH, the estimate yields |log/log r|\>>, since *\*q\r>. Hence =O> and =O>, for ,l>. It follows that *\*q|)>=O>. Instead of relying on GRH for bounding the >, we might have used the Bombieri\UVinogradov theorem, by adapting3.4>, that is based on the explicit bound proved in1.5>. With high probability, this yields a suitable at random among the primes in, say, the range from to . We first compute ,r|]>>\|)> rem -1|)>>> for ,l>> using Algorithm. We next reconstruct *\*q|)> rem -1|)>> using Chinese remaindering. In a similar way, we compute mod q*\*q|)> rem -1|)>>>. For any coefficient *x>> of such that \c> modulo for all i>>, we thus obtain both coefficients > and *e> modulo *\*q>. Whenever *\*q,c|)>>=1>>, this allows us to retrieve >. By considering a few more primes ,\,q>>> for sanity checking, we may thus compute > with high probability. The running time of this algorithm is bounded by >*T*log D|)>>. If the exponents ,\,e> are known, then we may also use the techniques from this section to further speed up the computation of the coefficients >,\,f>>. For this, let r\r\\>> be an increasing sequence of pairwise coprime numbers. Assume that the numerators and denominators of the coefficients >> are bounded by with >>, and assume that =O> for the first > terms of the sequence |)>1>>. For each, let > be the smallest prime number such that > divides -1>. Starting with =\=\=\>, we now repeat the following loop for > until \>q\B> for all ,t>: <\itemize> We compute |)> rem >-1|)>> using Algorithm. For ,t>, if \e> modulo > for all i>, then set \\\>. At the end of this loop, for ,t>, we may recover the coefficient >> from the coefficients of rem r>> in |)> rem >-1|)>>, where runs over >. With high probability, this algorithm runs in time >*t*log B|)>>, modulo the heuristic. For comparison, the best known complexity bound for univariate polynomials with integer coefficients is |)>>, unconditionally1.2>. The FFT based technique from this section is intended to be used in combination with aprobabilistic recombination method such as Algorithm or the technique from the previous subsection. In the most favorable case when we are always able to pick close to, we have seen that approximately > terms of are expected not collide modulo -1>. This means that we have evaluate at least *T> times in order to determine its sparse interpolation. Recall that only evaluations were required withAlgorithm. An interesting question is whether there exist any compromises or improvements. One idea is to exploit colliding terms better. For instance, assume that we are computing |]>>\f mod >-1|)>> for various moduli > and that |]>>> contains a term that comes from the sum of exactly two terms *x>> and*x>> with \e> modulo >. Now suppose that we are able to determine *x>> from the reductions |]>>> for a few other moduli \r>. Then we might reuse |]>>> to compute the term *x>> using a simple subtraction. Another idea is to optimize the (non-tangent) Graeffe transform variant of the geometric progression method. With the notations from section, we have seen how to achieve a complexity of >*log r|)>*log q|)>>. When taking =T>>, this complexity bound reduces to >|)>*log q|)>>, whereas only about |)>> terms of collide modulo -1>. For functions that can be evaluated \Pfast\Q (more specifically, this means that the > term dominates the term in the complexity bound() for geometric progression style algorithms), it may be interesting to switch to the FFT-based approach from this section. As a bonus, one may exploit the Frobenius map whenever we need to work over an extension field of . Under the assumption that >> and that suitable orders ,\,r>> as in section indeed exist, this leads to the complexity bound <\equation> O>*T*|)>*log |)>. We recall that this bound hides an increased constant factor with respect to(), so geometric progression based algorithms should remain faster whenever *log D>>. Under the assumption that >>, we also notice that() is always better than(). However, cyclic extension methods should remain faster in the super sparse case. One noticeable exception concerns sparse interpolation over the rationals, in which case we may use the approach from section. In this section we focus on the case when is a multivariate polynomial in ,\,x|]>> given as a black box function. In sections and, we first survey the well-known technique of Kronecker segmentation and known multivariate variants of Algorithm. In section, we recall a technique that we introduced in and in section we list a few generally useful ideas. In sections, , and we propose a few new ideas for particular situations. One first method for interpolating is to reduce to the univariate case using Kronecker segmentation. This reduction assumes that we have strict bounds <\equation*> deg> f\D,\,deg> f\D for the partial degrees of . Then we may consider the univariate polynomial <\equation*> g=f>,\,z*\*d>|)> with D*\*D>. We now compute the univariate sparse interpolation of and retrieve the sparse interpolation of by mapping each term *z> of to the term *x>*x|)> rem d>*\*x*\*d>> of . Another well known method for multivariate sparse interpolation is to adapt Algorithm. Using the notation =x>*\*x>> for -tuples \>, we may still write <\equation*> f=c*x>+\+c*x>. For a suitable vector \K>, we again evaluate at =\*\*\> for ,2*T-1>, and determine the roots >,\,\>> of > as in steps3 and4 of Algorithm for the univariate case. The only subtlety is that we need to chose > in a clever way, so that we can recover the exponents ,\,e> from >,\,\>>. In the case when has characteristic zero, Ben-Or and Tiwari proposed to pick =,\,p|)>>, where =2>, =3>, =5>, > is the sequence of prime numbers; see also. This allows to quickly recover each > from >> through smooth prime factorization. The approach still works for finite fields of characteristic p>. With respect to Kronecker segmentation, the approach is interesting if the total degree is much smaller than the sum +\+D> of the partial degrees. Note that we really need the weighted total degree <\equation*> D\max |)>*log p+\+|)>*log p\1\i\t|}> to be bounded by \log char K>. If is a finite field, then we may still use a variant of this idea, which partially goes back to. Modulo going to asuitable field extension, we first ensure that is of the form >>, where n> and D>. For some primitive element \>> over > and pairwise distinct points ,\,a>\\> we then take =,\,u-a|)>>. For sparse interpolations on GPUs, modular arithmetic is often most efficient for prime numbers of at most or bits. This may be insufficient for the interpolation of very large polynomials with 2> and thereby forces us to use >> instead. In that case, we may take =,\,p,u-a,\,u-a|)>> and ,\,a\\>>. We next require the total degree of in ,\,x> to be bounded by with \q> and we require to be linear in ,\,x>. If one of the above techniques applies, then the complexity analysis is similar to the one for the univariate case. Indeed, after the evaluation step2 of Algorithm, the values ,f|)>,\,f|)>> are in , so steps3, 4, and6 can be applied without changes, whereas the discrete logarithm computations in step5 are replaced by cheaper factorizations. Under the assumption that we can find a suitable divisor *r> of with \T> and > smooth, this leads to a complexity bound of the form <\equation> O>|)>*T*log |)>. For ,n-1|}>> and a random point \K>, consider the specialization <\equation*> \,\,x|)>\f,\,\,x,\,x|)>. With high probability on ,\,\>, the support of> coincides with the support of in,\,x>. The knowledge of the support of > may sometimes help with the computation of the support of, which in turn facilitates its sparse interpolation. Since > has less variables, we may recursively use a similar strategy for determining the support of >. This idea essentially goes back to Zippel in the case when and was applied successfully in in combination with more modern algorithms for sparse interpolation. Let us show how to implement this strategy in the case when we know a bound m>> for the total degree of with respect to ,\,x>. We first pick =,\,\|)>> in asimilar way as > in subsection: if or pm>>>, then we take \,\,p|)>>>; otherwise, we take \,\,u-a|)>> with the notations from there. The auxiliary function <\equation*> g,\,x|)>=f,\,x,\*x,\,\*x|)> admits the property that the terms *x>> of are in one to one correspondence with the terms *\>*\*\>*x>> of (where \|)>>). Moreover, knowing both > and *\>*\*\>> allows us to compute m>\,\,e|)>> through factorization. For a suitable element \K> of order T> and suitable integers ,\,\\0>, the next idea is to compute <\eqnarray*> >|>|>,\,z>|)> rem -1|)>>>|>|>|>,\,z>|)> rem -1|)>>>|>|>|>,\,z>|)> rem -1|)>.>>>> Mimicking Algorithm, we take *r\2*T> as large as possible under the conditions that =O> and > be smooth (ideally speaking, \r\T>). We compute the polynomial >, its roots, and the exponents \e,\,\\e> modulo as in Algorithm, which allows us to recover >. We compute > in a similar way (while exploiting the fact that the polynomial> is the same in both cases with high probability). Now consider some ,t|}>> such that \e\\e>> modulo and m>\em>\\m>\em>> modulo for all i>. Then we may recoverm>> from the terms *z\e rem r>> and *\m>>*z\e rem r>> of > and>. Moreover, \em> rem r=\e-\\em>|)> rem r> is an exponent of > (with high probability). Since > and > are both known, this allows us to recover m>>. Instead of a geometric progression method, we may also opt for the approach from section. In that case, we rather take an order T> close to , and we compute > and > using Algorithm. With high probability, this only leads to a constant portion of the exponents of . By varying the choices of and >, we may increase this portion until we find all exponents. This way of interpolating multivariate polynomials by packets of variables is particularly interesting if > admits significantly less terms than (so that the interpolation of> only requires a small portion of the total time) and m>> is small. In that case, we can typically avoid expensive field extensions as considered in sections, , and. If > is a finite field with T>, this leads to a complexity of >|)>*T*log q|)>> for the geometric progression method and >*T*log q|)>> when using FFTs (but with a larger constant factor for the dependency on ). <\remark> The algorithms in this subsection rely on the prime factoring approach from section. It is also possible to use Kronecker segmentation for the next packet of variables: see. Notice that Kronecker segmentation is mainly useful in combination with geometric progression style algorithms. It combines less well with the FFT based approach, except when and the coefficients of with respect to ,\,x> are dense polynomials in >. In the previous subsections, besides a bound on the total degree , we have also used bounds on the partial degrees ,\,D> and the total degree m>> with respect to a packet of variables ,\,x|}>>. More generally, for the sparse interpolation of a large multivariate polynomial , it is a good idea to first gather some general information about . This can often be done quickly and the information can be used to select a fast interpolation strategy for the specific polynomial at hand. For instance, in the previous subsection, it may help finding a packet of variables ,\,x|}>> that leads to optimal performance. Let us now describe a few probabilistic methods to gain information about . We may test whether through one evaluation at a random point \K>. We may test whether K> through one more evaluation at a random point \K> with \\>, by checking whether |)>=f|)>>. We may test whether is linear through a third evaluation at the point =|)>*\+\*\> for a random \K\>, by checking whether |)>=|)>*f|)>+\*f|)>>. More generally, if the total degree of is small, then we may determine this degree using evaluations of : we first set \f*z,\,\*z|)>> for random ,\,\>\K\>> so that . For \K>> of order d>, we then compute |)>> for > until |)>=g|)>> for the unique polynomial > of degree with |)>=G|)>> for ,i-1>. At that point, we have . The above method generalizes to the case when we wish to determine the degree of with respect to all variables in a subset of ,\,x|}>>. Indeed, modulo renaming variables, we may assume ,\,x|}>>. We then take \f*z,\,\*z>,\,\,\|)>>. We can obtain a rough estimate for the number of terms of as follows: for orders that increase geometrically and random integers ,\,\\0>, we compute \f>,\,z>|)> rem -1|)>>. As soon as exceeds , it becomes likely that contains less than terms. This may be used to determine a rough estimate of using > instead of > evaluations. We may recursively compute the set of variables > that occur in as follows. If is constant, then we return >. If is not constant and , then we return|}>>. Otherwise, we pick a random \K> and consider ,\,x|)>\f,\,x>,,\,\>|)>> with |n/2|\>>, as well as ,\,x|)>\f,\,\>,,\,x>|)>>. We recursively apply the method to and , and return the union of the output sets. This method uses at most |log n|\>*s> evaluations of for an output set of size . A is a subset ,\,x|}>> of variables such that is linear with respect to . Is it possible to quickly determine a maximal linear clique? Setting \deg f>\1|}>>> we must clearly have X> and we may always put all inactive variables inside . Although we do not have a fast algorithm to compute a maximal linear clique, the following randomized algorithm should be able to produce reasonably large ones for practical purposes (especially if we run it several times). Starting with \>, we iterate over all elements \X> in a random order, and add a new > to whenever remains linear with respect to |}>>. Notice that this algorithm only requires > evaluations of. More generally, we may use a variant of the linear clique algorithm to find large packets of variables ,\,x|}>> with respect to which admits a small total degree. This is the kind of packets that we need in section (by taking ,\,x|}>>\S> after a suitable permutation of indices). More precisely, for a finite field > of size , assume that we wish to compute apacket of size (as large as possible) such that the total degree > of with respect to satisfies >\q>. During the algorithm, we maintain a table that associates adegree> to each variable >. The table is initialized with \deg> f> for each . Starting with \>> we now repeat the following steps. We reorder indices such that ,\,x|}>> and \\\\>>. For >, we replace \deg|}>> f> and stop as soon as or \\> for some ,i|}>>. We then determine the index ,i|}>> for which > is minimal. If >\q>, then we add > to and continue; otherwise, we abort the main loop and return . During the loop, we may also remove all variables > for which > exceeds >. With this optimization, the number of evaluations of always remains bounded by*log q|)>>. Let us now consider the interesting special case when there exists a non-trivial linear clique ,\,x|}>> and assume that we know the support of with respect to ,\,x>>. Then the derivation strategy from section yields an efficient interpolation method for, even if is small. Indeed, in that case, we may write <\equation*> f=g,\,x|)>+g,\,x|)>*x+\+g,\,x|)>*x, where =\ f/\ x> for ,n-m>. The idea is to jointly interpolate ,g,\,g> using the fact that the vector ,g,\,g|)>> can be evaluated fast using automatic differentiation. As usual, we may either rely on a geometric progression style approach or on more directFFTs. When using a geometric progression style approach, we assume that \q> and let \K> denote aprimitive >-th root of unity. We pick random integers ,\,\\0>> and set \v>,\,z>|)>>. With high probability, the known exponents ,\,e> of are pairwise distinct modulo . We now compute ,\,V|)>> and interpolate the coefficients of using transposed multipoint interpolation. This yields and in time <\equation> O>|)>*u*log q|)>. When using FFTs, we rather take an -th primitive root of unity with > close to and only obtain a constant proportion of the coefficients of ; applying the same method for different values of , this leads to a runtime of >*u*log q|)>>>, with a higher constant factor with respect to . <\example> Consider the sparse interpolation of the determinant of an n> matrix with generic coefficients |)>i,j\n>>. Then the variables ,\,x> form a (maximal) linear clique. When specializing these variables by random elements in , we obtain a new function > for which the variables ,\,x> form a (maximal) linear clique. And so on for the successive rows of the matrix. Notice that the specialization of all variables > with k> has a support of size . Using the geometric progression strategy of complexity() in a recursive manner, it follows that we can compute the sparse interpolation of in time <\equation*> O>k\n>+|)>|)>*|)>*log q|)>=O>+|)>*n!*log q|)>=O>*n!*log q|)>. Notice that the traditional geometric progression method (of complexity()) takes time <\equation*> O>+|)>*n!*log q|)>=O>**n!*log q|)>. Instead of evaluating \f>,\,z>|)> rem -1|)>> by using Algorithm for random integers ,\,\\0> and a suitable modulus >, the multivariate setting also allows us to consider multivariate FFTs. To do so, we use moduli ,\,r\> and an integer matrix |)>n,j\m>>, and then take <\equation*> F,\,z|)>\f>*\*z>,\,z>*\*z>|)> rem >-1,\,z>-1|)>. Moreover, the > and > should be chosen in such a way that r*\*r> is close to , the linear map <\eqnarray*> >||\>>|/|r*\|\>\\\\/|r*\|\>>>|,\,e|)>>|>|*\+\+e*\,\,e*\+\+e*\|)>>>>> is surjective, and the exponents in the support of are mapped to essentially random elements in >. This generalization admits the advantage that the > do not need to be coprime, which may help to avoid working over large extensions >> if T>. It would be interesting to generalize the Frobenius FFT to this multivariate setting. Is it possible to speed up the interpolation if admits symmetries? As a simple example, let us consider the trivariate case with ,x,x|)>=f,x,x|)>>. It is easy to detect such a symmetry by checking whether ,\,\|)>=f,\,\|)>> holds for a random point \K>. When using multivariate FFTs as in the previous subsection, the idea is to pick parameters > and > of the form <\equation*> F,z,z|)>\f*z>,z*z>,z>*z>*z>|)> rem >-1,z>-1,z>-1|)>, so that is still symmetric with respect to > and >. Now let ,\\K> be primitive roots of unity of orders > and >. In order to evaluate of at all points >,\>,\>|)>> with ,i>=0,\,r-1> and =0,\,r-1>, it suffices to consider the points with \i>, thereby saving a factor of almost two. We may recover the coefficients of using an inverse DFT. This can again be done almost twice as efficiently as for non-symmetric polynomials, using a method from. Overall, we save a factor of almost two when computing sparse evaluations using the FFT method. This approach should generalize to more general types of symmetries, as the ones considered in. In order to interpolate a multivariate polynomial , it is recommended to first gather useful information about , as described in section. With this information at hand, we may then opt for a most appropriate interpolation strategy: <\itemize> If admits special properties, then we may wish to apply a dedicated algorithm, such as the linear clique strategy from section or the symmetric interpolation strategy from section. Whenever *\*D\q-1>, Kronecker segmentation can be combined with any of the univariate interpolation methods in order to obtain an efficient algorithm for the interpolation of . If *\*D\q-1>, but admits a reasonably small total degree, then we may use a variant of the geometric progression style algorithm from section in order to interpolate. If is normally sparse, but too large for one of the above strategies to apply, then interpolation by packets is a powerful tool. It typically allows us to reduce to the most favorable cases from sections and when we do not need any field extensions of > (except when T> or when admits not enough small prime divisors). As in the univariate case, we may then opt for a geometric progression style approach if =O> and for an FFT based approach if |)>>. If is super sparse, then we may use the strategy from section over the rationals and fall back on a cyclic extension style algorithm when is a given finite field. <\bibliography|bib|tm-plain|> <\bib-list|55> M.Agrawal, N.KayalN.Saxena. PRIMES is in P. , 160(2):781\U793, 2004. A.Arnold, M.GiesbrechtD.S.Roche. Sparse interpolation over finite fields via low-order roots of unity. , 27\U34. ACM Press, 2014. A.Arnold, M.GiesbrechtD.S.Roche. Faster sparse multivariate polynomial interpolation of straight-line programs. Symb. Comput.>, 75:4\U24, 2016. A.ArnoldD.S.Roche. Multivariate sparse interpolation using randomized Kronecker substitutions. , 35\U42. ACM Press, 2014. R.ArratiaL.Gordon. Tutorial on large deviations for the binomial distribution. , 51(1):125\U131, 1989. R.C.Baker, G.HarmanJ.Pintz. The difference between consecutive primes, II. , 83(3):532\U562, 2001. W.BaurV.Strassen. The complexity of partial derivatives. , 22:317\U330, 1983. M.Ben-OrP.Tiwari. A deterministic algorithm for sparse multivariate polynomial interpolation. , 301\U309. ACM Press, 1988. M.A.Bennett, G.Martin, K.O'BryantA.Rechnitzer. Explicit bounds for primes in arithmetic progressions. Math.>, 62(1\U4):427\U532, 2018. A.Bostan, G.LecerfÉ.Schost. Tellegen's principle into practice. , 37\U44. ACM Press, 2003. R.P.Brent, F.G.GustavsonD.Y.Y.Yun. Fast solution of Toeplitz systems of equations and computation of Padé approximants. Algorithms>, 1(3):259\U295, 1980. J.Canny, E.KaltofenY.Lakshman. Solving systems of non-linear polynomial equations faster. , 121\U128. ACM Press, 1989. D.G.CantorH.Zassenhaus. A new algorithm for factoring polynomials over finite fields. , 36(154):587\U592, 1981. A.DíazE.Kaltofen. FOXFOX: a system for manipulating symbolic objects in black box representation. , 30\U37. ACM Press, 1998. A.-M.Ernvall-HytönenN.Palojärvi. Explicit bound for the number of primes in arithmetic progressions assuming the Generalized Riemann Hypothesis. , 91:1317\U1365, 2022. T.S.Freeman, G.M.Imirzian, E.KaltofenY.Lakshman. DAGWOOD: a system for manipulating polynomials given by straight-line programs. , 14:218\U240, 1988. S.GargÉ.Schost. Interpolation of polynomials given by straight-line programs. , 410(27-29):2659\U2662, 2009. M.GiesbrechtD.S.Roche. Diversification improves interpolation. , 123\U130. ACM Press, 2011. P.Giorgi, B.GrenetA.Perret du Cray. Essentially optimal sparse polynomial multiplication. , ISSAC '20, 202\U209. New York, NY, USA, 2020. ACM Press. P.Giorgi, B.Grenet, A.Perret du CrayD.S.Roche. Random primes in arithmetic progressions. , arXiv, 2022. . P.Giorgi, B.Grenet, A.Perret du CrayD.S.Roche. Sparse polynomial interpolation and division in soft-linear time. , ISSAC '22, 459\U468. New York, NY, USA, 2022. ACM Press. B.Grenet, J.vander HoevenG.Lecerf. Randomized root finding over finite fields using tangent Graeffe transforms. , 197\U204. New York, NY, USA, 2015. ACM Press. B.Grenet, J.vander HoevenG.Lecerf. Deterministic root finding over finite fields using Graeffe transforms. , 27(3):237\U257, 2016. D.Y.GrigorievM.Karpinski. The matching problem for bipartite graphs with polynomially bounded permanents is in NC. , 166\U172. IEEE Computer Society, 1987. D.HarveyJ.vander Hoeven. Polynomial multiplication over finite fields in time >. ACM>, 69(2):1\U40, 2022. Article12. D.Harvey, J.vander HoevenG.Lecerf. Faster polynomial multiplication over finite fields. ACM>, 63(6), 2017. Article52. J.HeintzC.-P.Schnorr. Testing polynomials which are easy to compute. , 30, 237\U254. Geneva, 1982. Univ. Genève. J.vander HoevenR.Larrieu. The Frobenius FFT. , 437\U444. New York, NY, USA, 2017. ACM Press. J.vander Hoeven, R.LebretonÉ.Schost. Structured FFT and TFT: symmetric and lattice polynomials. , 355\U362. New York, NY, USA, 2013. ACM Press. J.vander HoevenG.Lecerf. On the bit-complexity of sparse polynomial and series multiplication. Symb. Comput.>, 50:227\U254, 2013. J.vander HoevenG.Lecerf. Sparse polynomial interpolation in practice. , 48(3/4):187\U191, 2015. J.vander HoevenG.Lecerf. Sparse polynomial interpolation. Exploring fast heuristic algorithms over finite fields. , HAL, 2019. . J.HuM.B.Monagan. A fast parallel sparse polynomial GCD algorithm. S.A.Abramov, E.V.ZimaX.-S.Gao, , 271\U278. New York, NY, USA, 2016. ACM Press. M.A.HuangA.J.Rao. Interpolation of sparse multivariate polynomials over large finite fields with applications. , 508\U517. Philadelphia, PA, USA, 1996. Society for Industrial and Applied Mathematics. M.-D.A.HuangA.J.Rao . Interpolation of sparse multivariate polynomials over large finite fields with applications. Algorithms>, 33(2):204\U228, 1999. Q.L.HuangX.S.Gao. Sparse Polynomial Interpolation with Finitely Many Values for the Coefficients. V.Gerdt, V.Koepf, W.SeilerE.Vorozhtsov, , 10490 Springer, Cham, 2017. Q.-L.Huang. Sparse polynomial interpolation over fields with large or zero characteristic. , 219\U226. New York, NY, USA, 2019. ACM Press. Q.L.HuangX.S.Gao. Faster interpolation algorithms for sparse multivariate polynomials given by straight-line programs. , 101:367\U386, 2020. M.JavadiM.Monagan. Parallel sparse polynomial interpolation over finite fields. M.Moreno MazaJ.-L.Roch, , 160\U168. New York, NY, USA, 2010. ACM Press. E.Kaltofen. Computing with polynomials given by straight-line programs I: greatest common divisors. , 131\U142. ACM Press, 1985. E.L.Kaltofen. Fifteen years after DSC and WLSS2 what parallel computations I do today: invited lecture at PASCO 2010. , 10\U17. ACM Press, 2010. E.Kaltofen, Y.N.LakshmanJ.-M.Wiley. Modular rational sparse multivariate polynomial interpolation. , 135\U139. New York, NY, USA, 1990. ACM Press. E.KaltofenW.Lee. Early termination in sparse interpolation algorithms. Symb. Comput.>, 36(3):365\U400, 2003. E.KaltofenB.M.Trager. Computing with polynomials given by black boxes for their evaluations: greatest common divisors, factorization, separation of numerators and denominators. Symb. Comput.>, 9(3):301\U320, 1990. E.KaltofenL.Yagati. Improved sparse multivariate polynomial interpolation algorithms. , 467\U474. Springer Verlag, 1988. R.Moenck. Fast computation of gcds. , 142\U171. New York, 1973. ACM Press. H.MuraoT.Fujise. Modular algorithm for sparse multivariate polynomial interpolation and its parallel implementation. Symb. Comput.>, 21:377\U396, 1996. M.Okamoto. Some inequalities relating to the partial sum of binomial probabilities. >, 10(1):29\U35, 1959. A.Perret du Cray. : interpolation, arithmétique, test d'identité>. , Université de Montpellier (France), 2023. R.Prony. Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastiques et sur celles de la force expansive de la vapeur de l'eau et de la vapeur de l'alkool, à différentes températures. de l'École Polytechnique Floréal et Plairial, an III>, 1(cahier 22):24\U76, 1795. J.-J.RislerF.Ronga. Testing polynomials. Symb. Comput.>, 10(1):1\U5, 1990. D.S.Roche. What can (and can't) we do with sparse polynomials? C.Arreche, , ISSAC '18, 25\U30. New York, NY, USA, 2018. ACM Press. A.Schönhage. Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2. , 7:395\U398, 1977. A.Sedunova. A partial Bombieri\UVinogradov theorem with explicit constants. , 101\U110, 2018. R.Zippel. Probabilistic algorithms for sparse polynomials. E.W.Ng, , 72, 216\U226. Springer-Verlag, 1979. <\initial> <\collection> <\attachments> <\collection> <\associate|bib-biblio> <\db-entry|+if2lDqI1qmMqAaI|inproceedings|CKL89> <|db-entry> E. Y. > <\db-entry|+2X2sGxMCqQBxuWk|inproceedings|BoLeSc:2003:tellegen> <|db-entry> G. É. > <\db-entry|+1LCfUIVm228oQhYU|inproceedings|GR11> <|db-entry> D. S. > <\associate|bib-bibliography> <\db-entry|+if2lDqI1qmMqAaI|inproceedings|CKL89> <|db-entry> E. Y. > <\db-entry|+1LCfUIVm228oQhYU|inproceedings|GR11> <|db-entry> D. S. > <\db-entry|+2X2sGxMCqQBxuWk|inproceedings|BoLeSc:2003:tellegen> <|db-entry> G. É. > <\db-entry|+tgCudvXqpzK6MJ|techreport|vdH:ffsparse> <|db-entry> G. > > <\db-entry|+2DPRky2cs1xL3Pp|article|Pro1795> <|db-entry> > de l'École Polytechnique Floréal et Plairial, an III> <\db-entry|+if2lDqI1qmMqAaW|inproceedings|BenOrTiwari1988> <|db-entry> P. > <\db-entry|+8lDHqURSvmeX09|inproceedings|HuangRao1996> <|db-entry> A. J. > <\db-entry|+1CdAZXTV2814h167|inproceedings|Kaltofen1985:stoc> <|db-entry> > <\db-entry|+2MmazzPzwkzLNWg|inproceedings|KaltofenLakshmanWiley1990> <|db-entry> Y. N. J.-M. > <\db-entry|+1QNR1I4u1FLRwI5l|article|KaltofenTrager1990> <|db-entry> B. M. > Symb. Comput.> <\db-entry|+2MmazzPzwkzLNWi|inproceedings|KaltofenYagati1988> <|db-entry> L. > <\db-entry|+1QNR1I4u1FLRwI5m|article|MF96> <|db-entry> T. > Symb. Comput.> <\db-entry|+1zXxMKe2u4Uc2Nj|inproceedings|DiazKaltofen1998> <|db-entry> E. > <\db-entry|+aBpkOnY2VUnXTuN|article|FrImKaLa88> <|db-entry> G. M. E. Y. > <\db-entry|+2Ge6dVXegcIOERX|inproceedings|ArnoldGiesbrechtRoche2014> <|db-entry> M. D. S. > <\db-entry|+tgCudvXqpzK6Ma|article|AGR16> <|db-entry> M. D. S. > Symb. Comput.> <\db-entry|+1CdAZXTV2814h166|inproceedings|ArnoldRoche2014> <|db-entry> D. S. > <\db-entry|+2MmazzPzwkzLNWZ|article|GargSchost2009> <|db-entry> É. > <\db-entry|+1QNR1I4u1FLRwI5Y|inproceedings|GiorgiGrenetPerretduCray2020> <|db-entry> B. A. > <\db-entry|+tgCudvXqpzK6N1|inproceedings|GiorgiGrenetPerretduCray2022> <|db-entry> B. A. D. S. > <\db-entry|+tgCudvXqpzK6NI|article|HuangRao1999> <|db-entry> A. J. Rao> Algorithms> <\db-entry|+1QNR1I4u1FLRwI5g|inproceedings|HuangGao2017> <|db-entry> X. S. > V. W. E. > <\db-entry|+1QNR1I4u1FLRwI5e|inproceedings|Huang19> <|db-entry> > <\db-entry|+1QNR1I4u1FLRwI5i|article|HuangGao2020> <|db-entry> L. X.S. > <\db-entry|+2Ge6dVXegcIOERe|inproceedings|Kaltofen2010:pasco> <|db-entry> > <\db-entry|+1QNR1I4u1FLRwI5h|article|HoevenLecerf2013> <|db-entry> G. > Symb. Comput.> <\db-entry|+8lDHqURSvmeXA0|article|vdH:spinpol> <|db-entry> G. > <\db-entry|+1QNR1I4u1FLRwI5d|inproceedings|HM16> <|db-entry> M. B. > E. V. X.-S. > <\db-entry|+1QNR1I4u1FLRwI5f|inproceedings|JaMo10> <|db-entry> M. > J.-L. > <\db-entry|+1QNR1I4u1FLRwI5j|article|KaltofenLee2003> <|db-entry> W. > Symb. Comput.> <\db-entry|+tgCudvXqpzK6N3|inproceedings|Roche2018> <|db-entry> > > <\db-entry|+2NThGgQo9xz5mXG|phdthesis|PerretduCray2023> <|db-entry> > : interpolation, arithmétique, test d'identité> <\db-entry|+tgCudvXqpzK6Mt|article|vdH:ffnlogn> <|db-entry> J. van der > >> ACM> 12> <\db-entry|+cJSkMhbTYXNoYI|article|Sch77> <|db-entry> > <\db-entry|+tgCudvXqpzK6Mm|article|vdH:ffmul> <|db-entry> J. van der G. > ACM> 52> <\db-entry|+11kKTCYk2Q6o1vKG|inproceedings|HeintzSchnorr1980> <|db-entry> C.-P. > <\db-entry|+1QNR1I4u1FLRwI5n|article|RislerRonga1990> <|db-entry> F. > Symb. Comput.> <\db-entry|+if2lDqI1qmMqAaV|article|BS83> <|db-entry> V. > <\db-entry|+tgCudvXqpzK6Mc|article|ArratiaGordon1989> <|db-entry> L. > <\db-entry|+1zXxMKe2u4Uc2Nv|article|Okamoto1959> <|db-entry> > > <\db-entry|+qjMTxGW1JFgEXMZ|article|BHP01> <|db-entry> G. J. > <\db-entry|+tgCudvXqpzK6MZ|article|AKS04> <|db-entry> N. N. > <\db-entry|+1LCfUIVm228oQhYT|article|BGY80> <|db-entry> F. G. D. Y. Y. > Algorithms> <\db-entry|+8lDHqURSvmeX2n|inproceedings|Moe73> <|db-entry> > <\db-entry|+tgCudvXqpzK6Md|article|CZ81> <|db-entry> H. > <\db-entry|+1QNR1I4u1FLRwI5Z|inproceedings|vdH:rmodroots> <|db-entry> J. van der G. > <\db-entry|+2DPRky2cs1xL3QK|article|vdH:dmodroots> <|db-entry> J. van der G. > <\db-entry|+tgCudvXqpzK6Mw|inproceedings|vdH:ffft> <|db-entry> R. > <\db-entry|+tgCudvXqpzK6NB|article|BennettMartinOBryantRechnitzer2018> <|db-entry> G. K. A. > Math.> <\db-entry|+tgCudvXqpzK6NE|article|ErnvallHytonenPalojarvi2022> <|db-entry> N. > <\db-entry|+tgCudvXqpzK6Mh|techreport|GiorgiGrenetPerretduCrayRoche2022> <|db-entry> B. A. D. S. > > <\db-entry|+tgCudvXqpzK6MP|article|Sedunova2018> <|db-entry> > <\db-entry|+if2lDqI1qmMqAaZ|inproceedings|GrigorievKarpinski1987> <|db-entry> M. > <\db-entry|+2DPRky2cs1xL3Q8|inproceedings|Zip79> <|db-entry> > > <\db-entry|+1QNR1I4u1FLRwI5a|inproceedings|vdH:symtft> <|db-entry> R. É. > <\references> <\collection> > > > > > |math-font-series|||H>>|11>> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > |\>|29>> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > <\auxiliary> <\collection> <\associate|bib> vdH:ffsparse Pro1795 BenOrTiwari1988 CKL89 HuangRao1996 Kaltofen1985:stoc KaltofenLakshmanWiley1990 KaltofenTrager1990 KaltofenYagati1988 MF96 DiazKaltofen1998 FrImKaLa88 ArnoldGiesbrechtRoche2014 AGR16 ArnoldRoche2014 GargSchost2009 GR11 GiorgiGrenetPerretduCray2020 GiorgiGrenetPerretduCray2022 HuangRao1999 HuangGao2017 Huang19 HuangGao2020 Kaltofen2010:pasco HoevenLecerf2013 vdH:spinpol HM16 JaMo10 KaltofenLee2003 Roche2018 PerretduCray2023 PerretduCray2023 GiorgiGrenetPerretduCray2022 vdH:ffnlogn Sch77 vdH:ffmul HeintzSchnorr1980 RislerRonga1990 AGR16 Huang19 BS83 Huang19 GargSchost2009 GargSchost2009 AGR16 GR11 GR11 GR11 ArratiaGordon1989 ArratiaGordon1989 Okamoto1959 BHP01 AKS04 GR11 BGY80 Moe73 BoLeSc:2003:tellegen CKL89 CZ81 vdH:rmodroots vdH:dmodroots vdH:rmodroots vdH:rmodroots vdH:ffft BennettMartinOBryantRechnitzer2018 BennettMartinOBryantRechnitzer2018 ErnvallHytonenPalojarvi2022 GiorgiGrenetPerretduCrayRoche2022 Sedunova2018 GiorgiGrenetPerretduCray2022 vdH:spinpol BenOrTiwari1988 GrigorievKarpinski1987 HuangRao1996 Zip79 vdH:spinpol vdH:spinpol vdH:symtft vdH:symtft <\associate|table> |1>|> Heuristic complexity bounds for the three main approaches to sparse interpolation. |> |2>|> Prime factorizations of |2-1>, |3-1>, and |5-1> for various small smooth values of |s>. |> |3>|> Prime factorizations of |q-1> and |a> for |q=1299742> and various values of |s>. |> <\associate|toc> |math-font-series||font-shape||1.Introduction> |.>>>>|> |1.1.Complexity considerations |.>>>>|> > |1.2.Overview of the paper |.>>>>|> > |math-font-series||font-shape||2.Preliminaries on finite fields> |.>>>>|> |math-font-series||font-shape||3.General observations> |.>>>>|> |3.1.From approximate to full interpolation |.>>>>|> > |3.2.Super sparse interpolation in large characteristic |.>>>>|> > |3.3.Conclusion |.>>>>|> > |math-font-series||font-shape||4.Univariate interpolation using cyclic extensions> |.>>>>|> |4.1.Complexity analysis |.>>>>|> > |4.2.Survey of existing variants based on cyclic extensions |.>>>>|> > |4.2.1.Determining the exponents using Chinese remaindering |.>>>>|> > |4.2.2.Composite moduli |.>>>>|> > |4.2.3.Diversification |.>>>>|> > |4.3.An optimized probabilistic algorithm based on diversification |.>>>>|> > |4.4.Analysis of the expected number of correct terms |.>>>>|> > |4.5.Probabilistic complexity analysis |.>>>>|> > |4.6.Estimating the number of terms t |.>>>>|> > |4.7.Conclusion |.>>>>|> > |math-font-series||font-shape||5.Univariate interpolation using geometric progressions> |.>>>>|> |5.1.Root finding |.>>>>|> > |5.2.Discrete logarithms |.>>>>|> > |5.3.Field extensions |.>>>>|> > |5.4.Exploiting the Frobenius map |.>>>>|> > |5.5.Traces |.>>>>|> > |5.6.Combining interpolations for several moduli |r> |.>>>>|> > |5.7.Conclusion |.>>>>|> > |math-font-series||font-shape||6.Univariate sparse interpolation using FFTs> |.>>>>|> |6.1.Fast evaluation modulo |x-1> |.>>>>|> > |6.2.Recombination into approximate sparse interpolations |.>>>>|> > |6.3.Example over |\> |.>>>>|> > |6.4.Example over |\> |.>>>>|> > |6.5.Sparse interpolation over the rationals |.>>>>|> > |6.6.Further remarks |.>>>>|> > |6.7.Conclusion |.>>>>|> > |math-font-series||font-shape||7.Multivariate sparse interpolation> |.>>>>|> |7.1.Reduction to the univariate case using Kronecker segmentation |.>>>>|> > |7.2.Generalizing algorithms based on geometric progressions |.>>>>|> > |7.3.Multivariate interpolation by packets of coordinates |.>>>>|> > |7.4.Gathering information |.>>>>|> > |Zero test |.>>>>|> > |Constant test |.>>>>|> > |Linearity test |.>>>>|> > |Small total degree |.>>>>|> > |Small partial degrees |.>>>>|> > |Number of terms |.>>>>|> > |Active variables |.>>>>|> > |Linear cliques |.>>>>|> > |Low degree packets |.>>>>|> > |7.5.Packets of total degree one |.>>>>|> > |7.6.Multivariate FFTs |.>>>>|> > |7.7.Symmetries |.>>>>|> > |7.8.Conclusion |.>>>>|> > |math-font-series||font-shape||Bibliography> |.>>>>|>