Diversification improves interpolation

Div ersification impro v es in terp olation Mark Giesbrec ht Daniel S. Ro c he Octob er 22, 2018 Abstract W e consider the problem of in terp olating an unkno wn multiv ariate p olynomial with co efficien ts taken from a finite field or as numerical appro ximations of complex num bers. Building on the recen t work of Garg and Schost, we impro ve on the b est-kno wn algorithm for inter- p olation ov er large finite fields b y presen ting a Las V egas randomized algorithm that uses fewer black b o x ev aluations. Using related tec h- niques, w e also address n umerical interpolation of sparse polynomials with complex co efficien ts, and provide the first prov ably stable algo- rithm (in the sense of relative error) for this problem, at the cost of mo destly more ev aluations. A key new technique is a randomization whic h mak es all co efficien ts of the unkno wn p olynomial distinguish- able, pro ducing what we call a diverse p olynomial. Another departure from most previous approaches is that our algorithms do not rely on ro ot finding as a subroutine. W e sho w how these impro vemen ts affect the practical p erformance with trial implementations. 1 In tro duction P olynomial in terp olation is a long-studied and imp ortan t problem in com- puter algebra and sym b olic computation. Given a wa y to ev aluate an un- kno wn p olynomial at an y chosen p oin t, and an upp er b ound on the degree, the interpolation problem is to determine a representation for the p oly- nomial. In sp arse in terpolation, we are also given an upp er b ound on the n umber of nonzero terms in the unkno wn polynomial, and the output is gen- erally returned in the sparse (also lacunary or sup ersp arse ) representation, wherein only the nonzero terms are explicitly stored. Applications of sparse interpolation include the manipulation and fac- torization of multiv ariate p olynomials and system solving (see, e.g., Canny et al. ( 1989 ); Kaltofen and T rager ( 1990 ); D ´ ıaz and Kaltofen ( 1995 , 1998 ); 1 Ja v adi and Monagan ( 2007 , 2009 ). With the adv ent of hybrid symbolic- n umeric algorithms for (systems of ) multiv ariate p olynomials with approx- imate co efficien ts, we find applications of approximate sparse interpolation, in particular for solving non-linear systems of equations (see, e.g., Sommese et al. ( 2001 , 2004 ); Stetter ( 2004 )) and factoring approximate multiv ariate p olynomials (see, e.g., Kaltofen et al. ( 2008 )). Sparse int erp olation is also a non-trivial generalization of the imp ortant problem of p olynomial identity testing : giv en a black b ox (esp ecially an algebraic circuit) computing an unkno wn p olynomial, determine whether the polynomial is zero. A relev an t result in our setting of sparse p olynomials is Bl¨ aser et al. ( 2009 ); for a more in-depth discussion, we recommend the recen t survey by Saxena ( 2009 ). Here w e examine the sparse in terp olation problem in t wo settings which ha ve received recent attention: when the co efficients are elemen ts of finite fields (particularly large finite fields, o ver which we hav e no choice) and when they are approximations to complex n umbers. W e give impro vemen ts o ver the state of the art in b oth of these cases, and demonstrate our new algorithms in practice with a full implemen tation in C++. 1.1 Problem definition Let F b e a field. A multiv ariate p olynomial f ∈ F [ x 1 , . . . , x n ] is said to b e t -sparse for some t ∈ N if f has at most t nonzero terms in the standard p o w er basis; that is, f can b e written f = t X i =1 c i x e i 1 1 x e i 2 2 · · · x e in n for co efficien ts c i ∈ F and exp onen t tuples ( e i 1 , . . . , e in ) ∈ N n for 1 ≤ i ≤ t . If each e i < d , then the size of this represen tation is O ( t ) field elemen ts plus O ( tn log d ) bits. W e seek algorithms whic h are p olynomial-time in the size of this represen tation. Let f ∈ F [ x 1 , . . . , x n ] ha ve degree less than d . A black b ox for f is a function which takes as input a vector ( a 1 , . . . , a n ) ∈ F n and pro duces f ( a 1 , . . . , a n ) ∈ F . The cost of the black b o x is the num b er of op erations in F required to ev aluate it at a giv en input. Clausen et al. ( 1991 ) sho wed that, if only ev aluations o ver the ground field F are allo wed, then for some instances at least Ω( n log t ) black b o x prob es are required. Hence if w e seek p olynomial-time algorithms, w e m ust extend the capabilities of the black b o x. T o this end, D ´ ıaz and Kaltofen ( 1998 ) 2 in tro duced the idea of an extende d domain black b ox which is capable of ev aluating f ( b 1 , . . . , b n ) ∈ E for any ( b 1 , . . . , b n ) ∈ E n where E is any exten- sion field of F . That is, w e c an change ev ery op eration in the black b ox to work ov er an extension field, usually paying an extra cost p er ev aluation prop ortional to the size of the extension. Motiv ated by the case of black b o xes that are division-free algebraic cir- cuits, w e will use the follo wing mo del whic h we b eliev e to b e fair and cov er all previous relev an t results. Here and for the remainder, M ( m ) is the n um b er of field op erations required to m ultiply tw o univ ariate p olynomials with degrees less than m , and O ˜ ( m ) represents any function b ounded by m (log m ) O (1) . Using Cantor and Kaltofen ( 1991 ), M ( m ) ∈ O ( m log m loglog m ), whic h is O ˜ ( m ). Definition 1.1. L et f ∈ F [ x 1 , . . . , x n ] and ` > 0 . A remainder blac k b o x for f with size ` is a pr o c e dur e which, given any monic squar e-fr e e p olynomial g ∈ F [ y ] with deg g = m , and any h 1 , . . . , h n ∈ F [ y ] with e ach deg h i < m , pr o duc es f ( h 1 , . . . , h n ) rem g using at most ` · M ( m ) op er ations in F . This definition is general enough to co ver the algorithms w e kno w of o v er finite fields, and w e submit that the cost mo del is fair to the standard black b o x, extended domain black b o x, and algebraic circuit settings. The mo del mak es sense ov er complex n umbers, as we will see. 1.2 In terp olation o v er finite fields W e first summarize previously known univ ariate in terp olation algorithms when F is a finite field with q elemen ts and identify our new contributions here. F or now, let f ∈ F q [ x ] hav e degree less than d and sparsit y t . W e will assume w e hav e a remainder blac k b o x for f with size ` . Since field elements can be represen ted with O (log q ) bits, a p olynomial-time algorithm will ha v e cost p olynomial in ` , t , log d , and log q . F or the dense output representation, one can use the classical metho d of Newton/W aring/Lagrange to interpolate in O ˜ ( `d ) time ( v on zur Gathen and Gerhard , 2003 , § 10.2). The algorithm of Ben-Or and Tiwari ( 1988 ) for sparse p olynomial in- terp olation, and in particular the version developed by Kaltofen and Y agati ( 1989 ), can b e adapted to arbitrary finite fields. Unfortunately , these algo- rithms require t discrete logarithm computations in F ∗ q , whose cost is small if the field size q is c hosen carefully (as in Kaltofen ( 2010 )), but not in general. F or arbitrary (and p oten tially large) q , we can tak e adv antage of the fact that each discrete logarithm that needs to b e computed falls in the range 3 Prob es Prob e degree Computation cost T otal cost Dense d 1 O ˜ ( d ) O ˜ ( `d ) Ben-Or & Tiw ari O ( t ) 1 O ( t 2 + t √ d ) O ˜ ( `t + t 2 + t √ d ) Garg & Sc host O ˜ ( t 2 log d ) O ˜ ( t 2 log d ) O ˜ ( t 4 log 2 d ) O ˜ ( `t 4 log 2 d ) Randomized G & S O ˜ ( t log d ) O ˜ ( t 2 log d ) O ˜ ( t 3 log 2 d ) O ˜ ( `t 3 log 2 d ) Ours O (log d ) O ˜ ( t 2 log d ) O ˜ ( t 2 log 2 d ) O ˜ ( `t 2 log 2 d ) T able 1: Sparse univ ariate in terp olation o ver large finite fields, with blac k b o x size ` , degree d , and t nonzero terms [0 , 1 , . . . , d − 1]. The “k angaro o metho d” of P ollard ( 1978 , 2000 ) can, with high probability , compute such a discrete log with O ( √ d ) field op erations. Using this algorithm mak es brings the total worst-case cost of Ben-Or and Tiw ari’s algorithm to O ( t` + t 2 + t √ d ). The curren t study builds most directly on the work of Garg and Sc host ( 2009 ), who ga ve the first p olynomial-time algorithm for sparse in terp olation o ver an arbitrary finite field. Their algorithm w orks roughly as follo ws. F or v ery small primes p , use the black b o x to compute f modulo x p − 1. A prime p is a “go o d prime” if and only if all the terms of f are still distinct mo dulo x p − 1. If we do this for all p in the range of roughly O ( t 2 log d ), then there will be sufficient go od primes to recov er the unique symmetric p olynomial ov er Z [ y ] whose ro ots are the exp onen ts of nonzero terms in f . W e then factor this p olynomial to find those exp onen ts, and correlate with any go o d prime image to determine the co efficien ts. The total cost is O ˜ ( `t 4 log 2 d ) field op erations. Using randomization, it is easy to reduce this to O ˜ ( `t 3 log 2 d ). Observ e that the co efficients of the symmetric integer p olynomial in Garg & Schost’s algorithm are b ounded by O ( d t ), which is muc h larger than the O ( d ) size of the exp onents ultimately recov ered. Our primary contribution o ver finite fields of size at least Ω( t 2 d ) is a new algorithm whic h av oids ev aluating the symmetric p olynomial and p erforming ro ot finding ov er Z [ y ]. As a result, w e reduce the total num b er of required ev aluations and develop a randomized algorithm with cost O ˜ ( `t 2 log 2 d ), which is roughly quadratic in the input and output sizes. Since this can b e deterministically v erified in the same time, our algorithm (as well as the randomized version of Garg & Sc host) is of the Las V egas type. The relev ant previous results men tioned ab o ve are summarized in T a- ble 1 , where we assume in all cases that the field size q is “large enough”. In the table, the “prob e degree” refers to the degree of g in eac h ev aluation 4 of the remainder blac k b o x as defined ab ov e. 1.3 Multiv ariate in terp olation An y of the univ ariate algorithms ab ov e can b e used to generate a multiv ari- ate p olynomial in terp olation algorithm in at least tw o differen t wa ys. F or what follo ws, write ρ ( d, t ) for the n um b er of remainder blac k b ox ev aluations required by some univ ariate interpolation algorithm, ∆( d, t ) for the degree of the remainder in eac h ev aluation, and ψ ( d, t ) for the n umber of other field op erations required b esides blac k b o x calls. Observe that these corresp ond to the first three columns in T able 1 . The first w ay to adapt a univ ariate interpolation algorithm to a mul- tiv ariate one is Kroneck er substitution: giv en a remainder blac k b o x for an unkno wn f ∈ F [ x 1 , . . . , x n ], with each partial degree less than d , we can easily construct a remainder black b o x for the univ ariate p olynomial ˆ f = f ( x, x d , x d 2 , . . . , x d n − 1 ) ∈ F [ x ], whose terms corresp ond one-to-one with terms of f . This is the approach taken for instance in Kaltofen ( 2010 , § 2) for the interpolation of m ultiv ariate p olynomials with rational co efficien ts. The cost is simply the cost of the chosen underlying univ ariate algorithm, with the degree increased to d n . The other metho d for constructing a multiv ariate interpolation algo- rithm is due to Zipp el ( 1990 ). The technique is inheren tly probabilistic and w orks v ariable-b y-v ariable, at eac h step solving a n umber of ˆ t × ˆ t transp osed V andermonde systems, for some ˆ t ≤ t . Sp ecifically , eac h system is of the form Ax = b , where A is a ˆ t × ˆ t matrix of scalars from the co efficien t field F q . The v ector v consists of the output of ˆ t remainder blac k b o x ev alua- tions, and so its elemen ts are in F q [ y ], and the system must b e solv ed modulo some g ∈ F q [ y ], as sp ecified by the underlying univ ariate algorithm. Observ e ho wev er that since A do es not con tain p olynomials, computing x = A − 1 b requires no mo dular p olynomial arithmetic. In fact, using the same tec h- niques as Kaltofen and Y agati ( 1989 , § 5), employing fast dense biv ariate p olynomial arithmetic, each system can b e solved using O  M  t · ∆( d, t )  · log  t · ∆( d, t )   field op erations. Eac h transp osed V andermonde system gives the remainder black b o x ev aluation of each of ˆ t univ ariate p olynomials that we are interpolating in that step. The num b er of such systems that must b e solv ed is therefore ρ ( d, t ), as determined by the underlying univ ariate algorithm. Finally , each 5 Kronec ker Zipp el Dense O ˜ ( `d n ) O ˜ ( `ntd ) Ben-Or & Tiw ari O ˜ ( `t + t 2 + td n/ 2 ) O ˜ ( nt 3 + nt 2 √ d + `nt 2 ) Garg & Sc host O ˜ ( `n 2 t 4 log 2 d ) O ˜ ( `nt 5 log 2 d ) Randomized G & S O ˜ ( `n 2 t 3 log 2 d ) O ˜ ( `nt 4 log 2 d ) Ours O ˜ ( `n 2 t 2 log 2 d ) O ˜ ( `nt 3 log 2 d ) T able 2: Sparse multiv ariate interpolation ov er large finite fields, with blac k b o x size ` , n v ariables, degree d , and t nonzero terms of the ˆ t univ ariate interpolations pro ceeds with the given ev aluations. The total cost, o ver all iterations, is O ˜  `nt · ∆( d, t ) · ρ ( d, t )  field op erations for the remainder black b ox ev aluations, plus O ˜  ntψ ( d, t ) + `nt · ∆( d, t )  field op erations for additional computation. Zipp el ( 1990 ) used the dense algorithm for univ ariate in terp olation; using Ben-Or and Tiwari’s algorithm instead w as studied by Kaltofen and Lee ( 2003 ). T able 2 summarizes the cost of the univ ariate algorithms mentioned ab o v e applied to sparse multiv ariate interpolation ov er a sufficiently large finite field, using Kronec ker’s and Zipp el’s metho ds. F or completeness, w e men tion a few more results on closely related prob- lems that do not hav e a direct b earing on the current study . Grigoriev et al. ( 1990 ) give a parallel algorithm with small depth but whic h is not com- p etitiv e in our mo del due to the large num b er of pro cessors required. A practical parallel version of Ben-Or and Tiwari’s algorithm has b een dev el- op ed by Ja v adi and Monagan ( 2010 ). Kaltofen et al. ( 1990 ) and Avenda˜ no et al. ( 2006 ) presen t mo dular algorithms for in terp olating p olynomials with rational and integer co efficien ts. How ev er, their metho ds do not seem to apply to finite fields. 1.4 Appro ximate P olynomial In terp olation In Section 4 w e consider the case of appro ximate sparse interpolation. Our goal is to provide b oth a numerically more robust practical algorithm, but also the first algorithm whic h is prov ably n umerically stable, with no heuris- tics or conjectures. W e define an “  -approximate black b o x” as one which 6 ev aluates an unknown t -sparse target polynomial f ∈ C [ x ], of degree d , with relativ e error at most  > 0. Our goal is to build a t -sparse p olynomial g suc h that k f − g k 2 ≤  k f k 2 . A b ound on the degree and sparsity of the tar- get p olynomial, as w ell as  , must also b e pro vided. In Section 4 w e formally define the ab o ve problem, and demonstrate that the problem of sparse in- terp olation is w ell-p osed. W e then adapt our v ariant of the Garg and Sc host ( 2009 ) algorithm for the approximate case, prov e it is numerically accurate in terms of the relative error of the output, and analyze its cost. W e also presen t a full implementation in Section 5 and v alidating exp erimen ts. Recen tly , a n umber of numerically-focussed sparse in terp olation algo- rithms ha ve been presen ted. The algorithm of Giesbrec ht et al. ( 2009 ) is a numerical adaptation of Ben-Or and Tiwari ( 1988 ), which samples f at O ( t ) randomly chosen ro ots of unit y ω ∈ C on the unit circle. In particular, ω is c hosen to ha ve (high) order at least the degree, and a randomization sc heme is used to av oid clustering of nodes which will cause dramatic ill- conditioning. A relatively w eak theoretical b ound is prov en there on the randomized conditioning scheme, though exp erimental and heuristic evi- dence suggests it is muc h b etter in practice. Cuyt and Lee ( 2008 ) adapt Rutishauser’s q d algorithm to alleviate the need for b ounds on the partial degrees and the sparsit y , but still ev aluate at high-order ro ots of unity . Ap- pro ximate sparse rational function in terp olation is considered by Kaltofen and Y ang ( 2007 ) and Kaltofen et al. ( 2007 ), using the Structured T otal Least Norm (STLN) metho d and, in the latter, randomization to impro ve conditioning. Approximate sparse interpolation is also considered for inte- ger p olynomials by Mansour ( 1995 ), where a p olynomial-time algorithm is presen ted in quite a differen t mo del from ours. In particular the ev aluation error is absolute (not relativ e) and the complexity is sensitiv e to the bit length of the in teger co efficien ts. Note that all these w orks ev aluate the polynomial only on the unit circle. This is necessary b ecause we allow and exp ect f to ha ve v ery large degree, whic h would cause a catastrophic loss of precision at data points of non-unit magnitude. Similarly , we assume that the complex argument of ev aluation p oin ts is exactly sp ecified, which is again necessary b ecause an y error in the argumen t would b e exp onen tially magnified b y the degree. The primary con tribution of the w ork in this pap er is to provide an algorithm with b oth rigorously prov able relative error and go o d practical p erformance. Our algorithm typically requires O ˜ ( t 2 log 2 d ) ev aluations at primitiv e ro ots of unity of order O ˜ ( t 2 log d ) (as opp osed to order d in pre- vious approac hes). W e guarantee that it finds a t -sparse p olynomial g such that k g − f k 2 ≤ 2  k f k 2 . An exp erimen tal demonstration of the numerical 7 robustness is giv en in Section 5 . 2 Sparse in terp olation for generic fields Here and for the remainder, we say a p olynomial f is t -sp arse if it can b e written as a sum of at most t nonzero coefficients times a monomial. W e assume the unknown p olynomial f is alwa ys univ ariate. This is without loss of generality , as w e can use the Kronec ker substitution as discussed ab o v e. The exp onen tial increase in the univ ariate degree only corresp onds to a factor of n increase in log deg f , and since our algorithms will ultimately ha ve cost p olynomial in log deg f , p olynomial time is preserved. Assume a fixed, unkno wn, t -sparse univ ariate p olynomial f ∈ F [ x ] with degree at most d . W e will use a remainder black b o x for f to ev aluate f rem ( x p − 1) for small primes p . W e say p is a “go o d prime” if the sparsit y of f rem ( x p − 1) is the same as that of f itself — that is, none of the exp onen ts are equiv alent mo dulo p . The follo wing lemma sho ws the size of primes required to randomly c ho ose go od primes with high probability . Lemma 2.1. L et f ∈ F [ x ] b e a t -sp arse p olynomial with de gr e e d , and let λ = max  21 ,  5 3 t ( t − 1) ln d  . A prime chosen at r andom in the r ange λ, . . . , 2 λ is a go o d prime for f with pr ob ability at le ast 1 / 2 . Pr o of. W rite e 1 , . . . , e t for the exp onen ts of nonzero terms in f . If p is a bad prime, then p divides ( e j − e i ) for some i < j . Eac h e j − e i ≤ d , so there can b e at most log λ d = ln d/ ln λ primes that divide each e j − e i . There are exactly  t 2  suc h pairs of exp onents, so the total num ber of bad primes is at most ( t ( t − 1) ln d ) / (2 ln λ ). F rom Rosser and Schoenfeld ( 1962 , Corollary 3 to Theorem 2), the total n umber of primes in the range λ, . . . , 2 λ is at least 3 λ/ (5 ln λ ) when λ ≥ 21, which is at least t ( t − 1) ln d/ ln λ , at least twice the n umber of bad primes. No w observe an easy case for the sparse interpolation problem. If a p olynomial f ∈ F [ x ], has all co efficients distinct; that is, f = P 1 ≤ i ≤ t c i x e i and c i = c j ⇒ i = j , then w e sa y f is diverse . T o interpolate a diverse p olynomial f ∈ F [ x ], we first follo w the metho d of Garg and Sc host ( 2009 ) b y computing f rem ( x p i − 1) for “go o d primes” p i suc h that the sparsity of f rem ( x p i − 1) is the same as that of f . Since f is div erse, f rem ( x p i − 1) is also div erse and in fact each mo dular image has the same set of co efficien ts. 8 Algorithm 1: Generic interpolation Input : µ ∈ R > 0 , T , D , q ∈ N , and a remainder blac k b o x for unknown T -sparse f ∈ F [ x ] with deg f < D Output : t ∈ N , e 1 , . . . , e t ∈ N , and c 1 , . . . , c t ∈ F suc h that f = P 1 ≤ i ≤ t c i x e i 1 t ← 0 2 λ ← max  21 ,  5 3 T ( T − 1) ln D  3 for d log 2 (3 /µ ) e primes p ∈ { λ, . . . , 2 λ } do 4 Use blac k b o x to compute f p = f ( x ) rem ( x p − 1) 5 if f p has mor e than t terms then 6 t ← sparsit y of f p 7 % ← p 8 α ← elemen t of F such that Pr[ f ( αx ) is not diverse] < µ/ 3 9 g % ← f ( α x ) rem ( x % − 1) 10 c 1 , . . . , c t ← nonzero co efficien ts of g % 11 e 1 , . . . , e t ← 0 12 for d 2 ln(3 /µ ) + 4(ln D ) / (ln λ ) e primes p ∈ { λ, . . . , 2 λ } do 13 Use blac k b o x to compute g p = f ( α x ) rem ( x p − 1) 14 if g p has exactly t nonzer o terms then 15 for i = 1 , . . . , t do Up date e i with exp onent of c i in g p mo dulo p via Chinese remaindering 16 for i = 1 , . . . , t do c i ← c i α − e i 17 return f ( x ) = P 1 ≤ i ≤ t c i x e i Using this fact, we av oid the need to construct and subsequently factor the symmetric p olynomial in the exp onen ts. Instead, w e correlate like terms based on the (unique) co efficien ts in each mo dular image, then use simple Chinese remaindering to construct each exp onen t e i from its image mo dulo eac h p i . This requires only O (log d ) remainder blac k box ev aluations at go od primes, gaining a factor of t impro vemen t ov er the randomized v ersion of Garg and Sc host ( 2009 ) for diverse p olynomials. In the following sections, w e will show how to choose an α ∈ F so that f ( α x ) — whic h we can easily construct a remainder black b o x for — is div erse. With such a pro cedure, Algorithm 1 gives a Monte Carlo algorithm for in terp olation o ver a general field. 9 Theorem 2.2. With inputs as sp e cifie d, Algorithm 1 c orr e ctly c omputes the unknown p olynomial f with pr ob ability at le ast 1 − µ . The total c ost in field op er ations (exc ept for step 8 ) is O  ` ·  log D log T + loglog D + log 1 µ  · M  T 2 log D   . Pr o of. The for lo op on line 3 searches for the true sparsity t and a single go od prime % . Since eac h prime p in the giv en range is go od with probabilit y at least 1 / 2 by Lemma 2.1 , the probability of failure at this stage is at most µ/ 3. The for lo op on line 12 searches for and uses sufficiently many go o d primes to recov er the exp onen ts of f . The pro duct of all the go o d primes m ust b e at least D , and since eac h prime is at least λ , at least (ln D ) / (ln λ ) go od primes are required. Let n = d 2 ln(3 /µ ) + 4(ln D ) / (ln λ ) e b e the num b er of primes sampled in this lo op, and k = d (ln D ) / (ln λ ) e the n umber of go o d primes required. W e can derive that ( n/ 2 − k ) 2 ≥ (ln(3 /µ ) + k ) 2 > ( n/ 2) ln(3 /µ ), and therefore exp( − 2( n 2 − k ) 2 /n ) < µ/ 3. Using Ho effding’s Inequalit y ( Ho effding , 1963 ), this means the probability of encountering fewer than k go od primes is less than µ/ 3. Therefore the total probability of failure is at most µ . F or the cost analysis, the dominating cost will b e the mo dular blac k b o x ev aluations in the last for lo op. The n umber of ev aluations in this lo op is O (log(1 /µ ) + (log D ) / (log λ )), and eac h ev aluation has cost O ( ` · M ( λ )). Since the size of eac h prime is Θ((log D ) / (log T + loglog D )), the complexity bound is correct as stated. In case the bound T on the num b er of nonzero terms is v ery bad, w e could c ho ose a smaller v alue of λ based on the true sparsit y t b efore line 8 , impro ving the cost of the remainder of the algorithm. In addition, as our b ound on p ossible num b er of “bad primes” seems to b e quite lo ose, a more efficient approach in practice w ould b e to replace the for lo op on line 12 with one that starts with a prime muc h smaller than λ and incrementally searches for the next larger primes until the pro duct of all go od primes is at least D . W e could c ho ose the low er b ound to start searc hing from based on low er b ounds on the birthda y problem. That is, assuming (falsely) that the exp onen ts are randomly distributed mo dulo p , start with the least p that will ha ve no exp onen ts collide modulo p with high probabilit y . This would yield an algorithm more sensitive to the true b ound on bad primes, but unfortunately giv es a worse formal cost analysis. 10 3 Sparse in terp olation o v er finite fields W e now examine the case that the ground field F is the finite field with q elemen ts, which w e denote F q . First w e show ho w to effectively div ersify the unknown p olynomial f in order to complete Algorithm 1 for the case of large finite fields. Then we show ho w to extend this to a Las V egas algorithm with the same complexit y . 3.1 Div ersification F or an unkno wn f ∈ F q [ x ] given by a remainder black b o x, we must find an α so that f ( αx ) is diverse. A surprisingly simple trick works: ev aluating f ( α x ) for a random nonzero α ∈ F q . Theorem 3.1. F or q ≥ T ( T − 1) T +1 and any T -sp arse p olynomial f ∈ F q [ x ] with deg f < D , if α is chosen uniformly at r andom fr om F ∗ q , the pr ob ability that f ( α x ) is diverse is at le ast 1 / 2 . Pr o of. Let t ≤ T b e the exact n umber of nonzero terms in f , and write f = P 1 ≤ i ≤ t c i x e i , with nonzero co efficien ts c i ∈ F ∗ q and e 1 < e 2 < · · · < e t . So the i th co efficien t of f ( αx ) is c i α e i . If f ( αx ) is not div erse, then we m ust hav e c i α e i = c j α e j for some i 6 = j . Therefore consider the p olynomial A ∈ F q [ y ] defined b y A = Y 1 ≤ i 0 and assume ther e exists an unknown t -sp arse f ∈ C [ x ] of de gr e e at most D . An  -appro ximate black b o x for f takes an input ξ ∈ C and pr o duc es a γ ∈ C such that | γ − f ( ξ ) | ≤  | f ( ξ ) | . That is, the relative error of any ev aluation is at most  . As noted in the in tro duction, we will sp ecify our input p oints exactly , at (relatively lo w order) ro ots of unity . The appr oximate sp arse univariate interp olation pr oblem is then as follows: given D, T ∈ N and δ ≥  > 0, and an  - appro ximate black b o x for an unknown T -sparse p olynomial f ∈ C [ x ] of degree at most D , find a T -sparse polynomial g ∈ C [ x ] such that k f − g k 2 ≤ δ k g k 2 . The following theorem shows that t -sparse p olynomials are well-defined b y go o d ev aluations on the unit circle. Theorem 4.3. L et  > 0 and f ∈ C [ x ] b e a t -sp arse p olynomial. Supp ose ther e exists a t -sp arse p olynomial g ∈ C [ x ] such that for a prime p which is go o d for b oth f and f − g , and p th primitive r o ot of unity ω ∈ C , we have | f ( ω i ) − g ( ω i ) | ≤  | f ( ω i ) | for 0 ≤ i < p. Then k f − g k 2 ≤  k f k 2 . Mor e over, if g 0 ∈ C [ x ] is forme d fr om g by deleting al l the terms not in the supp ort of f , then k f − g 0 k 2 ≤ 2  k f k 2 . Pr o of. Summing ov er p o w ers of ω w e hav e X 0 ≤ i

0 such that | c i | ≥ δ k f k 2 for 1 ≤ i ≤ t . The p olynomial f is said to b e  -diverse if and only if every p air of distinct c o efficients is at le ast  k f k 2 ap art. That is, for every 1 ≤ i < j ≤ t , | c i − c j | ≥  k f k 2 . In tuitively , if ( / 2) corresp onds to the mac hine precision, this means that an algorithm can reliably distinguish the co efficien ts of a  -div erse p olyno- mial. W e now sho w how to choose a random α to guarantee  -diversit y . Theorem 4.7. L et δ ≥  > 0 and f ∈ C [ x ] a t -sp arse p olynomial whose non- zer o c o efficients ar e of magnitude at le ast δ k f k 2 . If s is a prime satisfying s > 12 and t ( t − 1) ≤ s ≤ 3 . 1 δ  , then for ζ = e 2 π i /s an s -PRU and k ∈ N chosen uniformly at r andom fr om { 0 , 1 , . . . , s − 1 } , f ( ζ k x ) is  -diverse with pr ob ability at le ast 1 2 . Pr o of. F or each 1 ≤ i ≤ t , write the co efficien t c i in p olar notation to base ζ as c i = r i ζ θ i , where each r i and θ i are nonnegative real n umbers and r i ≥ δ k f k 2 . Supp ose f ( ζ k x ) is not  -diverse. Then there exist indices 1 ≤ i < j ≤ t suc h that    r i ζ θ i ζ ke i − r j ζ θ j ζ ke j    ≤  k f k 2 . Because min( r i , r j ) ≥ δ k f k 2 , the v alue of the left hand side is at least δ k f k 2 ·   ζ θ i + ke i − ζ θ j + ke j   . Dividing out ζ θ j + ke i , w e get    ζ θ i − θ j − ζ k ( e j − e i )    ≤  δ . 17 By wa y of contradiction, assume there exist distinct choices of k that satisfy the ab o ve inequalit y , sa y k 1 , k 2 ∈ { 0 , . . . , s − 1 } . Since ζ θ i − θ j and ζ e j − e i are a fixed p o w ers of ζ not dep ending on the c hoice of k , this means    ζ k 1 ( e j − e i ) − ζ k 2 ( e j − e i )    ≤ 2  δ . Because s is prime, e i 6 = e j , and we assumed k 1 6 = k 2 , the left hand side is at least | ζ − 1 | . Observe that 2 π /s , the distance on the unit circle from 1 to ζ , is a go o d appro ximation for this Euclidean distance when s is large. In particular, since s > 12, | ζ − 1 | 2 π /s > √ 2  √ 3 − 1  / 2 π / 6 , and therefore | ζ − 1 | > 6 √ 2( √ 3 − 1) /s > 6 . 2 /s , whic h from the statement of the theorem is at least 2 /δ . This is a contradiction, and therefore the assumption w as false; namely , there is at most one c hoice of k such that the i ’th and j ’th co efficients collide. Then, since there are exactly  t 2  distinct pairs of co efficien ts, and s ≥ t ( t − 1) = 2  t 2  , f ( ζ k x ) is div erse for at least half of the choices for k . W e note that the diversification which maps f ( x ) to f ( ζ k x ) and back is n umerically stable since ζ is on the unit circle. In practice, the previous theorem will b e far to o p essimistic. W e there- fore prop ose the metho d of Algorithm 5 to adaptively choose s , δ , and ζ k sim ultaneously , given a go od prime p . Supp ose there exists a threshold S ∈ N such that for all primes s > S , a random s th primitiv e ro ot of unity ζ k mak es f ( ζ k x )  -div erse with high probabilit y . Then Algorithm 5 will return a ro ot of unity whose order is within a constan t factor of S , with high probability . F rom the previous theorem, if such an S exists it m ust b e O ( t 2 ), and hence the num b er of iterations required is O (log t ). Otherwise, if no such S exists, then we cannot diversify the p olynomial. Roughly sp eaking, this corresp onds to the situation that f has to o many co efficien ts with absolute v alue close to the machine precision. In this case, w e can simply use the algorithm of Garg and Sc host ( 2009 ) numerically , ac hieving the same stability but using a greater num b er of ev aluations and bit op erations. It is p ossible to establish an adaptive h ybrid b et ween our algorithm and that of Garg and Schost ( 2009 ) b y making f as  -diverse as p ossible given our precision. The non-zero co efficien ts of f are clustered in to groups which are not  -diverse (i.e., are within  k f k 2 of each other). W e 18 Algorithm 5: Adaptive diversification Input :  -appro ximate black b o x for f , known go o d prime p , known sparsit y t Output : ζ , k such that f ( ζ k x ) is  -div erse, or FAIL 1 s ← 1, δ ← ∞ , f p ← 0 2 while s ≤ t 2 and # { c o effs c of f s s.t. | c | ≥ δ } < t do 3 s ← least prime ≥ 2 s 4 ζ ← exp(2 π i /s ) 5 k ← random integer in { 0 , 1 , . . . , s − 1 } 6 Compute f s = f ( ζ k x ) rem x p − 1 7 δ ← least num b er s.t. all co efficien ts of f s at least δ in absolute v alue are pairwise  -distinct 8 if δ > 2  then return FAIL 9 else return ζ k can use the symmetric p olynomial reconstruction of Garg and Schost ( 2009 ) to extract the exp onen ts within eac h group. 4.4 Appro ximate in terp olation algorithm W e now plug our  -appro ximate remainder black b o x, and metho d for mak- ing f  -diverse, in to our generic Algorithm 1 to complete our algorithm for appro ximate interpolation. Theorem 4.8. L et δ > 0 , f ∈ C [ x ] with de gr e e at most D and sp arsity at most T , and supp ose al l nonzer o c o efficients c of f satisfy | c | > δ k f k 2 . Supp ose also that  < 1 . 5 δ / ( T ( T − 1)) , and we ar e given an  -appr oximate black b ox for f . Then, for any µ < 1 / 2 we have an algorithm to pr o duc e a g ∈ C [ x ] satisfying the c onditions of The or em 4.3 . The algorithm suc c e e ds with pr ob ability at le ast 1 − µ and uses O ˜ ( T 2 · log (1 /µ ) · log 2 D ) black b ox evaluations and flo ating p oint op er ations. Pr o of. Construct an appro ximate remainder black b o x for f using Algo- rithm 4 . Then run Algorithm 1 using this blac k b o x as input. On step 8 of Algorithm 1 , run Algorithm 5 , iterating steps 5 – 7 d log 2 (3 /µ ) e times on eac h iteration through the while lo op to c ho ose a diversifying α = ζ k with probabilit y at least 1 − µ/ 3. The cost comes from Theorems 2.2 and 4.5 along with the previous discussion and Theorem 4.7 . 19 Observ e that the resulting algorithm is Monte Carlo, but could b e made Las V egas b y combining the finite fields zero testing algorithm discussed in Section 3.2 with the guaran tees of Theorem 4.3 . 5 Implemen tation results W e implemen ted our algorithms in C++ using the GNU Multiple Precision Arithmetic Library (GMP , http://gmplib.org/ ) and Victor Shoup’s Num- b er Theory Library (NTL, http://www.shoup.net/ntl/ ) for the exp onen t arithmetic. F or comparison with the algorithm of Garg and Schost ( 2009 ), w e also used NTL’s squarefree p olynomial factorization routines. W e note that, in our exp erimen ts, the cost of integer p olynomial factorization (for Garg & Sc host) and Chinese remaindering were alwa ys negligible. In our timing results, “Determ” refers to the deterministic algorithm as stated in Garg and Sc host ( 2009 ) and “Alg 1” is the algorithm w e ha ve presen ted here ov er finite fields, without the v erification step. W e also de- v elop ed and implemen ted a more adaptive, Monte Carlo version of these algorithms, as briefly describ ed at the end of Section 2 . The basic idea is to sample mo dulo x p − 1 for just one prime p ∈ Θ( t 2 log d ) that is go od with high probabilit y , then to search for m uch smaller go o d primes. This go o d prime search starts at a lo wer b ound of order Θ( t 2 ) based on the birthda y problem, and finds consecutiv ely larger primes un til enough primes hav e b een found to reco ver the symmetric p olynomial in the exp onen ts (for Garg & Schost) or just the exp onen ts (for our metho d). The corresp onding im- pro ved algorithms are referred to as “G&S MC” and “Alg 1++” b elo w, resp ectiv ely . T able 3 summarizes some timings for these four algorithms ov er the finite field Z / 65521 Z . This mo dulus was c hosen for con venience of implemen ta- tion, although other metho ds such as the Ben-Or and Tiw ari algorithm migh t b e more efficien t in this particlar field since discrete logarithms could b e computed quickly . How ever, observ e that our algorithms (and those from Garg and Schost ) hav e only p oly-logarithmic dep endence on the field size, and so will ev entually dominate. The timings are giv en in seconds of CPU time on a 64-bit AMD Phenom I I 3.2GHz pro cessor with 512K/2M/6M cache, compiled using GCC 4.4.3 with the -O3 flag. Note that the num b ers listed reflect the b ase-2 lo garithm of the degree b ound and the sparsity bound for the randomly-generated test cases. The timings are mostly as exp ected based on our complexit y estimates, 20 log 2 D T Determ G&S MC Alg 1 Alg 1++ 12 10 3.77 0.03 0.03 0.01 16 10 46.82 0.11 0.11 0.08 20 10 — 0.38 0.52 0.33 24 10 — 0.68 0.85 0.38 28 10 — 1.12 2.35 0.53 32 10 — 1.58 2.11 0.66 12 20 37.32 0.15 0.02 0.02 16 20 — 0.91 0.52 0.28 20 20 — 3.5 3.37 1.94 24 20 — 6.59 5.94 2.99 28 20 — 10.91 10.22 3.71 32 20 — 14.83 16.22 4.24 12 30 — 0.31 0.01 0.01 16 30 — 3.66 1.06 0.65 20 30 — 10.95 6.7 3.56 24 30 — 25.04 12.42 9.32 28 30 — 38.86 19.36 13.8 32 30 — 62.53 68.1 14.66 12 40 — 0.58 0.01 0.02 16 40 — 8.98 3.7 1.54 20 40 — 30.1 12.9 8.42 24 40 — 67.97 38.34 16.57 28 40 — — 73.69 36.24 32 40 — — — 40.79 T able 3: Finite Fields Algorithm Timings 21 Noise Mean Error Median Error Max Error 0 4.440 e − 16 4.402 e − 16 8.003 e − 16 ± 10 − 12 1.113 e − 14 1.119 e − 14 1.179 e − 14 ± 10 − 9 1.149 e − 11 1.191 e − 11 1.248 e − 11 ± 10 − 6 1.145 e − 8 1.149 e − 8 1.281 e − 8 T able 4: Approximate Algorithm Stability and also confirm our suspicion that primes of size O ( t 2 ) are sufficient to a void exp onen t collisions. It is satisfying but not particularly surprising to see that our “Alg 1++” is the fastest on all inputs, as all the algorithms ha ve a similar basic structure. Had we compared to the Ben-Or and Tiwari or Zipp el’s metho d, they w ould probably be more efficient for small sizes, but w ould be easily b eaten for large degree and arbitrary finite fields as their costs are sup er-polynomial. The implemen tation of the appro ximate algorithm uses mac hine double precision (IEEE), the built-in C++ complex type, and the p opu- lar F astest F ourier T ransform in the W est (FFTW, http://www.fftw.org/ ) pac k age for computing FFTs. Our stabilit y results are summarized in T a- ble 4 . Each test case was randomly generated with degree at most 2 20 and at most 50 nonzero terms. W e v aried the precision as sp ecified in the table and ran 10 tests in each range. Observe that the error in our results w as often less than the  error on the ev aluations themselves. Both implemen tations are released under an MIT-style licence and are a v ailable from the second author’s w ebsite at http://www.cs.uwaterloo.ca/ ~ droche/diverse/ . 6 Conclusions W e ha ve sho wn how to use the idea of diversification to improv e the complex- it y of sparse interpolation ov er large finite fields by a factor of t , the num ber of nonzero terms. W e achiev e a similar complexity for approximate sparse in terp olation, and provide the first prov ably n umerically stable algorithm for this purp ose. Our exp erimen ts confirm these theoretical results. Numerous op en problems remain. A primary shortcoming of our algo- rithms is the quadratic dep endence on t , as opp osed to linear in the case of dense in terp olation or even sparse interpolation in smaller or chosen finite fields using the Ben-Or and Tiwari algorithm. It seems that reducing this 22 quadratic dep endency will not b e p ossible without a different approach, b e- cause of the birthday problem embedded in the div ersification step. In the appro ximate case, a pro v ably n umerically stable algorithm for sparse inter- p olation with only O ( t ) prob es is still an op en question. And, while general bac kward error stability is not p ossible in the high degree case, it would b e in teresting in the case of low degree and many v ariables. Ac kno wledgemen ts W e thank Reinhold Burger and ´ Eric Schost for p ointing out an error in an earlier draft. The commen ts and suggestions of the anonymous referees were also v ery helpful, in particular regarding connections to previous results and the pro of of Theorem 3.1 . References Mart ´ ın Av enda ˜ no, T eresa Krick, and Ariel P acetti. Newton-Hensel interpo- lation lifting. F ound. Comput. Math. , 6(1):81–120, 2006. Eric Bach and Jeffrey Shallit. A lgorithmic numb er the ory. Vol. 1 . F ounda- tions of Computing Series. MIT Press, Cam bridge, MA, 1996. Efficient algorithms. Mic hael Ben-Or and Praso on Tiw ari. A deterministic algorithm for sparse m ultiv ariate p olynomial interpolation. In Pr o c. STOC’88 , pages 301–309, New Y ork, NY, USA, 1988. A CM. Markus Bl¨ aser, Moritz Hardt, Richard J. Lipton, and Nisheeth K. Vish- noi. Deterministically te sting sparse p olynomial identities of un b ounded degree. Information Pr o c essing L etters , 109(3):187 – 192, 2009. John Cann y , Erich Kaltofen, and Lakshman Y agati. Solving systems of nonlinear p olynomial equations faster. In Pr o c. ISSAC’89 , pages 121– 128, 1989. Da vid G. Cantor and Eric h Kaltofen. On fast multiplication of p olynomials o ver arbitrary algebras. A cta Informatic a , 28:693–701, 1991. Mic hael Clausen, Andreas Dress, Johannes Grabmeier, and Marek Karpin- ski. On zero-testing and interpolation of k-sparse m ultiv ariate p olynomials o ver finite fields. The or etic al Computer Scienc e , 84(2):151 – 164, 1991. 23 Annie Cuyt and W en-shin Lee. A new algorithm for sparse interpolation of m ultiv ariate p olynomials. The or etic al Computer Scienc e , 409(2):180–185, 2008. Angel D ´ ıaz and Eric h Kaltofen. On computing greatest common divisors with polynomials given by black b o xes for their ev aluations. In Pr o c. ISSA C’95 , pages 232–239, 1995. Angel D ´ ıaz and Eric h Kaltofen. F OXBO X: a system for manipulating sym- b olic ob jects in black b o x represen tation. In Pr o c. ISSAC’98 , pages 30–37, 1998. Sanc hit Garg and ´ Eric Schost. In terp olation of p olynomials given b y straigh t-line programs. The or etic al Computer Scienc e , 410(27-29):2659 – 2662, 2009. J. von zur Gathen and J. Gerhard. Mo dern Computer Algebr a . Cam bridge Univ ersity Press, Cambridge, second edition, 2003. Mark Giesbrech t and Daniel S. Ro che. Detecting lacunary p erfect p o w ers and computing their ro ots. Journal of Symb olic Computation , 2010. T o app ear; preprint at . Mark Giesbrec ht, George Labahn, and W en-shin Lee. Sym b olic-n umeric sparse in terp olation of multiv ariate p olynomials. Journal of Symb olic Computation , 44(8):943 – 959, 2009. Dima Y u. Grigoriev, Marek Karpinski, and Michael F. Singer. F ast paral- lel algorithms for sparse multiv ariate polynomial interpolation ov er finite fields. SIAM J. on Computing , 19(6):1059–1063, 1990. W assily Ho effding. Probability inequalities for sums of b ounded random v ariables. J. Amer. Statist. Asso c. , 58:13–30, 1963. Sa yed Mohammad Mahdi Ja v adi and Mic hael Monagan. On factorization of multiv ariate p olynomials ov er algebraic num b er and function fields. In Pr o c. ISSAC’09 , pages 199–206, 2009. Sey ed Jav adi and Mic hael Monagan. P arallel sparse polynomial interpolation o ver finite fields. In Pr o c. Intl. Wkshp. Par al lel and Symb olic Computation (P ASCO) , pages 160–168, 2010. Sey ed Mohammad Mahdi Jav adi and Mic hael Monagan. A sparse mo dular GCD algorithm for p olynomials o ver algebraic function fields. In Pr o c. ISSA C’07 , pages 187–194, 2007. 24 E. Kaltofen and Z. Y ang. On exact and appro ximate interpolation of sparse rational functions. In Pr o c. ISSAC’07 , pages 11–18, 2007. doi: 10.1145/ 1277548.1277577. E. Kaltofen, Y. N. Lakshman, and J.-M. Wiley . Mo dular rational sparse m ultiv ariate polynomial in terp olation. In Pr o c. ISSAC’90 , pages 135–139, New Y ork, NY, USA, 1990. A CM. E. Kaltofen, Z. Y ang, and L. Zhi. On probabilistic analysis of randomization in hybrid symbolic-numeric algorithms. In Pr o c. Workshop on Symb olic- Numeric Computation (SNC 2007) , pages 203–210, 2007. doi: 10.1145/ 1277500.1277503. Eric h Kaltofen and W en-shin Lee. Early termination in sparse interpolation algorithms. J. Symb olic Comput. , 36(3-4):365–400, 2003. Eric h Kaltofen and Barry M. T rager. Computing with p olynomials given b y blac k b o xes for their ev aluations: Greatest common divisors, factor- ization, separation of numerators and denominators. Journal of Symb olic Computation , 9:301–320, 1990. Eric h Kaltofen and Lakshman Y agati. Improv ed sparse m ultiv ariate p olyno- mial in terp olation algorithms. In Pr o c. ISSAC’88 , pages 467–474, 1989. Eric h Kaltofen, John P . Ma y , Zhenfeng Y ang, and Lihong Zhi. Approximate factorization of m ultiv ariate p olynomials using singular v alue decomp osi- tion. Journal of Symb olic Computation , 2008. Eric h L. Kaltofen. Fifteen years after DSC and WLSS2: What parallel com- putations I do to day . In Pr o c. Intl. Wkshp. Par al lel and Symb olic Com- putation (P ASCO) , pages 10–17, 2010. doi: 10.1145/1837210.1837213. Y. Mansour. Randomized approximation and interpolation of sparse p oly- nomials. SIAM Journal on Computing , 24(2):357–368, 1995. J. M. Pollard. Mon te Carlo metho ds for index computation (mo d p ). Math. Comp. , 32(143):918–924, 1978. doi: 10.1090/S0025- 5718- 1978- 0491431- 9. J. M. Pollard. Kangaro os, monop oly and discrete logarithms. Journal of Cryptolo gy , 13:437–447, 2000. doi: 10.1007/s001450010010. J. B. Rosser and L. Schoenfeld. Appro ximate formulas for some functions of prime n umbers. Il l. J. Math. , 6:64–94, 1962. 25 Nitin Saxena. Progress on p olynomial identit y testing. Bul l. EA TCS , 99: 49–79, 2009. Andrew J. Sommese, Jan V ersc helde, and Charles W. W ampler. Numer- ical decomp osition of the solution sets of p olynomial systems in to irre- ducible comp onen ts. SIAM Journal on Numeric al Analysis , 38(6):2022– 2046, 2001. Andrew J. Sommese, Jan V ersc helde, and Charles W. W ampler. Numerical factorization of m ultiv ariate complex p olynomials. The or etic al Computer Scienc e , pages 651–669, 2004. Hans J. Stetter. Numeric al Polynomial Algebr a . SIAM, 2004. Ric hard Zippel. In terp olating p olynomials from their v alues. Journal of Symb olic Computation , 9(3):375 – 403, 1990. 26