Interpolation of Shifted-Lacunary Polynomials

INTERPOLA TION OF SHIFTED-LA CUNAR Y POL YNOMIALS Mark Giesbrecht and D aniel S. R oche Abstract. Giv en a “blac k box” function to ev aluate an unkno wn ra- tional p olynomial f ∈ Q [ x ] at p oin ts mo dulo a prime p , we exhibit algo- rithms to compute the represen tation of the p olynomial in the sparsest shifted p o wer basis. That is, we determine the sparsit y t ∈ Z > 0 , the shift α ∈ Q , the exponents 0 ≤ e 1 < e 2 < · · · < e t , and the co efficien ts c 1 , . . . , c t ∈ Q \ { 0 } such that f ( x ) = c 1 ( x − α ) e 1 + c 2 ( x − α ) e 2 + · · · + c t ( x − α ) e t . The computed sparsity t is absolutely minimal ov er an y shifted pow er basis. The no v elty of our algorithm is that the complexity is p olyno- mial in the (sparse) representation size, whic h may b e logarithmic in the degree of f . Our method combines previous celebrated results on sparse interpolation and computing sparsest shifts, and provides a wa y to handle polynomials with extremely high degree which are, in some sense, sparse in information. Keyw ords. Sparse in terp olation, Sparsest shift, Lacunary p olynomials Sub ject classification. Primary 68W30; Secondary 12Y05 1. Introduction In terp olating an unknown p olynomial from a set of ev aluations is a problem whic h has in terested mathematicians for h undreds of years, and which is now implemen ted as a standard function in most computer algebra systems. T o illustrate some differen t kinds of interpolation problems, consider the follo wing three representations for a p olynomial f of degree n : f ( x ) = a 0 + a 1 x + a 2 x 2 + · · · + a n x n , (1.1) f ( x ) = b 0 + b 1 x d 1 + b 2 x d 2 + · · · + b s x d s , (1.2) f ( x ) = c 0 + c 1 ( x − α ) e 1 + c 2 ( x − α ) e 2 + · · · + c t ( x − α ) e t . (1.3) 2 Giesbrec ht & Ro c he In (1.1) w e see the dense represen tation of the polynomial, where all co efficien ts (ev en zero es) are represented. Newton and W aring discov ered metho ds to in- terp olate f in time prop ortional to the size of this representation in the 18 th cen tury . The sparse or lacunary represen tation is shown in (1.2) , wherein only the terms with non-zero co efficien ts are written (with the p ossible exception of the constant coefficient b 0 ). Here w e sa y that f is s -sp arse b ecause it has exactly s non-zero and non-constan t terms; the constan t co efficien t requires sp ecial treatment in our algorithms regardless of whether or not it is zero, and so w e do not count it tow ards the total n umber of terms. Ben-Or & Tiwari ( 1988 ) disco vered a metho d to in terp olate in time p olynomial in the size of this represen tation. Kaltofen & Lee ( 2003 ) present and analyze very efficien t algo- rithms for this problem, in theory and practice. The Ben-Or & Tiwari metho d has also b een examined in the context of approximate (floating p oin t) p olyno- mials b y Giesbrec ht et al. ( 2006 ), where the similarit y to the 1795 method of de Pron y is also p oin ted out. Bl¨ aser et al. ( 2009 ) consider the more basic problem of iden tity testing of sparse p olynomials ov er Q , and present a deterministic p olynomial-time algorithm. In (1.3) , f is written in the shifted p o w er basis 1 , ( x − α ) , ( x − α ) 2 , . . . , and w e sa y that α is a t -sp arse shift of f b ecause the represen tation has exactly t non-zero and non-constan t terms in this basis. When α is chosen so that t is absolutely minimal in (1.3) , w e call this the sp arsest shift of f . W e presen t new algorithms to interpolate f ∈ Q [ x ], given a blac k b o x for ev aluation, in time prop ortional to the size of the shifted-lacunary represen tation corresp onding to (1.3) . It is easy to see that t could b e exp onen tially smaller than b oth n and s , for example when f = ( x + 1) n , demonstrating that our algorithms are pro viding a significant improv emen t in complexit y o v er those previously known, whose running times are p olynomial in n and s . The main applications of all these metho ds for p olynomial interpolation are signal pro cessing and reducing intermediate expression sw ell. Dense and sparse in terp olation hav e b een applied successfully to b oth these ends, and our new algorithms effectively extend the class of p olynomials for whic h suc h applications can b e made. The most significan t c hallenge here is computing the sparsest shift α ∈ Q . Computing this v alue from a set of ev aluation points was stated as an op en problem b y Boro din & Tiw ari ( 1991 ). An algorithm for a generalization of our problem in the dense representation w as given b y Grigoriev & Karpinski ( 1993 ), though its cost is exp onen tial in the size of the output; they admit that the dep endency on the degree of the p olynomial is probably not optimal. Our algorithm achiev es deterministic p olynomial-time complexit y for p olynomials In terp olation of Shifted-Lacunary Polynomials 3 o ver the rational n umbers. W e are alwa ys careful to count the bit c omplexity — the n umber of fixed-precision mac hine op erations — and hence account for an y co efficien t gro wth in the solution or in termediate expressions. The black b o x mo del we use is slightly mo dified from the traditional one: p ∈ N , θ ∈ Z p - - f ( θ ) mo d p f ( x ) ∈ Q [ x ] Giv en a prime p and an element θ in Z p , the blac k b o x computes the v alue of the unkno wn polynomial ev aluated at θ o v er the field Z p . (An error is produced exactly in those unfortunate circumstances that p divides the denominator of f ( θ ).) W e generally refer to this as a mo dular black b ox . T o accoun t for the reasonable p ossibilit y that the cost of black b o x calls dep ends on the size of p , w e define κ f to b e an upp er b ound on the n um b er of field op erations in Z p used in blac k b o x ev aluation, for a giv en p olynomial f ∈ Q [ x ]. Some kind of extension to the standard black b o x, such as the mo dular blac k b o x prop osed here, is in fact necessary , since the v alue of a p olynomial of degree n at an y p oin t other than 0 , ± 1 will t ypically hav e n bits or more. Th us, any algorithm whose complexit y is prop ortional to log n cannot p erform suc h an ev aluation ov er Q or Z . Other p ossibilities might include allowing for ev aluations on the unit circle in some represen tation of a subfield of C , or returning only a limited n umber of bits of precision for an ev aluation. T o be precise ab out our notion of size, first define size( q ) for q ∈ Q to be the num b er of bits needed to represent q . So if we write q = a b with a ∈ Z , b ∈ N , and gcd( a, b ) = 1, then size( q ) = d log 2 ( | a | + 1) e + d log 2 ( b + 1) e + 1. F or a rational p olynomial f as in (1.3) , define: (1.4) size( f ) = size( α ) + t X i =0 size( c i ) + t X i =1 size( e i ) . W e will often employ the following upp er b ound for simplicity: (1.5) size( f ) ≤ size( α ) + t ( H ( f ) + log 2 n ) , where H ( f ) is defined as max 0 ≤ i ≤ t size( c i ). Our algorithms will hav e p olynomial complexity in the smallest p ossible size( f ). F or the complexit y analysis, we use the normal notion of a “mul- tiplication time” function M ( n ), which is the n umber of field op erations re- quired to compute the pro duct of p olynomials with degrees less than n , or 4 Giesbrec ht & Ro c he in tegers with sizes at most n . W e alw ays assume that M ( n ) ∈ Ω( n ) and M ( n ) ∈ O ( n 2 ). Using the results from Cantor & Kaltofen ( 1991 ), w e can write M ( n ) ∈ O ( n log n log log n ). The remainder of the pap er is structured as follows. In Section 2 we sho w ho w to find the sparsest shift from ev aluation p oin ts in Z p , where p is a prime with some sp ecial prop erties pro vided b y some “oracle”. In Section 3 w e sho w ho w to p erform sparse interpolation given a mo dular black b o x for a poly- nomial. In Section 4 w e sho w how to generate primes suc h that a sufficient n umber satisfy the conditions of our oracle. Section 5 provides the complexity analysis of our algorithms. W e conclude in Section 6 , and introduce some op en questions. 2. Computing the Sparsest Shift F or a polynomial f ∈ Q [ x ], w e first fo cus on computing the sparsest shift α ∈ Q so that f ( x + α ) has a minimal n um b er of non-zero and non-constant terms. This information will later b e used to recov er a representation of the unkno wn p olynomial. 2.1. The p olynomial f ( p ) . Here, and for the remainder of this pap er, for a prime p and f ∈ Q [ x ], define f ( p ) ∈ Z p [ x ] to b e the unique p olynomial with degree less than p which is equiv alen t to f mo dulo x p − x and with all co efficien ts reduced mo dulo p . F rom F ermat’s Little Theorem, w e then see immediately that f ( p ) ( α ) ≡ f ( α ) mo d p for all α ∈ Z p . Hence f ( p ) can b e found by ev aluating f at eac h point 0 , 1 , . . . , p − 1 mo dulo p and using dense in terp olation ov er Z p [ x ]. Notice that, o ver Z p [ x ], ( x − α ) p ≡ x − α mo d x p − x , and therefore ( x − α ) e i ≡ ( x − α ) k for an y k 6 = 0 such that e i ≡ k mo d ( p − 1). The smallest such k is in the range { 1 , 2 , . . . , p } ; we now define this with some more notation. F or a ∈ Z and p ositiv e in teger m , define a rem 1 m to b e the unique in teger in the range { 1 , 2 , . . . , m } which is congruent to a mo dulo m . As usual, a rem m denotes the unique congruent in teger in the range { 0 , 1 , . . . , m − 1 } . If f is as in (1.3) , then b y reducing term-by-term w e can write (2.1) f ( p ) ( x ) = ( c 0 rem p ) + t X i =1 ( c i rem p )( x − α p ) e i rem 1 ( p − 1) , where α p is defined as α rem p . Hence, for some k ≤ t , α p is a k -sparse shift for f ( p ) . That is, the p olynomial f ( p ) ( x + α p ) ov er Z p [ x ] has at most t non-zero and non-constant terms. In terp olation of Shifted-Lacunary Polynomials 5 Computing f ( p ) from a mo dular blac k b o x for f is straigh tforward. First, use p blac k-b o x calls to determine f ( i ) rem p for i = 0 , 1 , . . . , p − 1. Recalling that κ f is the n um b er of field op erations in Z p for each black-box call, the cost of this step is O ( pκ f M (log p )) bit op erations. Second, w e use the w ell-known divide-and-conquer metho d to interpolate f ( p ) in to the dense representation (see, e.g., Boro din & Munro ( 1975 , Section 4.5)). Since deg f ( p ) < p , this step has bit complexit y O ( M ( p ) M (log p ) log p ). F urthermore, for any α ∈ Z p , the dense representation of f ( p ) ( x + α ) can b e computed in exactly the same wa y as the second step ab o ve, simply by shifting the indices of the already-ev aluated p oin ts by α . This immediately gives a na ¨ ıv e algorithm for computing the sparsest shift of f ( p ) : compute f ( p ) ( x + γ ) for γ = 0 , 1 , . . . , p − 1, and return the γ that minimizes the n um b er of non-zero, non- constan t terms. The bit complexity of this approac h is O ( p log p M ( p ) M (log p )), whic h for our applications will often b e less costly than the more sophisticated approac hes of, e.g., Lakshman & Saunders ( 1996 ) or Giesbrech t et al. ( 2003 ), precisely b ecause p will not b e v ery m uc h larger than deg f ( p ) . 2.2. Ov erview of Approach. W e will mak e rep eated use of the follo wing fundamen tal theorem from Lakshman & Saunders ( 1996 ): F act 2.2. Let F b e an arbitrary field and f ∈ F [ x ] , and supp ose α ∈ F is such that f ( x + α ) has t non-zero and non-constan t terms. If deg f ≥ 2 t + 1 then α is the unique sparsest shift of f . F rom this w e can see that, if α is the unique sparsest shift of f , then α p = α rem p is the unique sparsest shift of f ( p ) pr ovide d that deg f ( p ) ≥ 2 t + 1. This observ ation provides the basis for our algorithm. The input to the algorithms will b e a mo dular blac k b o x for ev aluating a rational p olynomial, as describ ed ab o v e, and bounds on the maximal size of the unkno wn p olynomial. Note that suc h b ounds are a necessity in an y type of blac k-b o x interpolation algorithm, since otherwise w e could nev er b e sure that the computed polynomial is really equal to the blac k-b o x function at every p oin t. Sp ecifically , we require B A , B T , B H , B N ∈ N such that size( α ) ≤ B A , t ≤ B T , size( c i ) ≤ B H , for 0 ≤ i ≤ t, log 2 n ≤ B N . 6 Giesbrec ht & Ro c he By considering the following p olynomial: c ( x − α ) n + ( x − α ) n − 1 + · · · + ( x − α ) n − t +1 , w e see that these b ounds are indep enden t — that is, none is polynomially- b ounded b y the others — and therefore are all necessary . W e are now ready to present the algorithm for computing the sparsest shift α almost in its entiret y . The only part of the algorithm left unsp ecified is an or acle which, based on the v alues of the b ounds, pro duces primes to use. W e wan t primes p suc h that deg f ( p ) ≥ 2 t + 1, which allo ws us to reco ver one mo dular image of the sparsest shift α . But since w e do not kno w the exact v alue of t or the degree n of f ov er Q [ x ], w e define some prime p to b e a go o d prime for sp arsest shift c omputation if and only if deg f ( p ) ≥ min { 2 B T + 1 , n } . F or the remainder of this section, “go od prime” means “goo d prime for sparsest shift computation.” Our oracle indicates when enough primes ha v e b een pro duced so that at least one of them is guaran teed to ha v e b een a go o d prime, which is necessary for the pro cedure to terminate. The details of how to construct such an oracle will b e considered in Section 4 . Algorithm 2.3. Computing the sparsest shift. Input: ◦ A mo dular black b o x for an unknown p olynomial f ∈ Q [ x ] ◦ Bounds B A , B T , B H , B N ∈ N as described ab o ve ◦ An oracle whic h pro duces primes and indicates when at least one go od prime must ha ve b een pro duced Output: A sparsest shift α of f . 1. P ← 1, G ← ∅ 2. While log 2 P < 2 B A + 1 do 3–14 3. p ← new prime from the oracle 4. Ev aluate f ( i ) rem p for i = 0 , 1 , . . . , p − 1 5. Use dense in terp olation to compute f ( p ) 6. If deg f ( p ) ≥ 2 B T + 1 then 7. Use dense in terp olation to compute f ( p ) ( x + γ ) for γ = 1 , 2 , . . . , p − 1 8. α p ← the unique sparsest shift of f ( p ) 9. P ← P · p , G ← G S { p } 10. Else if P = 1 and oracle indicates ≥ 1 go o d prime has b een pro duced then 11. q ← least prime suc h that log 2 q > 2 B T B A + B H (computed directly) 12. Ev aluate f ( i ) rem q for i = 0 , 1 , . . . , 2 B T In terp olation of Shifted-Lacunary Polynomials 7 13. Compute f ∈ Q [ x ] with deg f ≤ 2 B T b y dense in terp olation in Z q [ x ] follo wed b y rational reconstruction on the co efficien ts 14. Return A sparsest shift α computed b y a univ ariate algorithm from Giesbrec ht et al. ( 2003 ) on input f 15. Return The unique α = a/b ∈ Q such that | a | , b ≤ 2 B A and a ≡ bα p mo d p for each p ∈ G , using Chinese remaindering and rational reconstruction Theorem 2.4. With inputs as sp ecified, Algorithm 2.3 correctly returns a sparsest shift α of f . Pr oof. Let f , B A , B T , B H , B N b e the inputs to the algorithm, and suppose t, α are as sp ecified in (1.3) . First, consider the degenerate case where n ≤ 2 B T , i.e., the b ound on the sparsit y of the sparsest shift is at least half the actual degree of f . Then, since eac h f ( p ) can ha ve degree at most n (regardless of the c hoice of p ), the condition of Step 6 will never b e true. Hence Steps 10–14 will even tually b e executed. The size of coefficients o ver the standard p o wer basis is b ounded b y 2 B T B A + B H since deg f ≤ 2 B T , and therefore f will be correctly computed on Step 5 . In this case, F act 2.2 may not apply , i.e. the sparsest shift ma y not b e unique, but the algorithms from Giesbrech t et al. ( 2003 ) will still pro duce a sparsest shift of f . No w suppose instead that n ≥ 2 B T + 1. The oracle even tually pro duces a go od prime p , so that deg f ( p ) ≥ 2 B T + 1. Since t ≤ B T and f ( p ) has at most t non-zero and non-constant terms in the ( α rem p )-shifted p o w er basis, the v alue computed as α p on Step 8 is exactly α rem p , by F act 2.2 . The v alue of P will also be set to p > 1 here, and can only increase. So the condition of Step 10 is nev er true. Since the numerator and denominator of α are b oth b ounded ab o v e by 2 B A , w e can use rational reconstruction to compute α once we ha ve the image mo dulo P for P ≥ 2 2 B A +1 . Therefore, when we reac h Step 15 , w e ha ve enough images α p to recov er and return the correct v alue of α .  W e still need to specify which algorithm to use to compute the sparsest shift of a densely-represen ted f ∈ Q [ x ] on Step 14 . T o mak e Algorithm 2.3 completely deterministic, w e should use the univ ariate sym b olic algorithm from Giesbrec ht et al. ( 2003 , Section 3.1), although this will ha ve v ery high complex- it y . Using a probabilistic algorithm instead giv es the follo wing, whic h follo ws directly from the referenced w ork. 8 Giesbrec ht & Ro c he Theorem 2.5. If the “t w o pro jections” algorithm of Giesbrec ht et al. ( 2003 , Section 3.3) is used on Step 14 , then Steps 10–14 of Algorithm 2.3 can b e p er- formed with O ( B 2 T M ( B 4 T B A + B 3 T B H )) bit op erations, plus O ( κ f B T M ( B T B A + B H )) bit op erations for the blac k-b o x ev aluations. The precise complexit y analysis pro ving that the en tire Algorithm 2.3 has bit complexity polynomial in the b ounds given depends heavily on the size and n umber of primes p that are used, and so m ust b e p ostp oned until Section 5.1 , after our discussion on c ho osing primes. Example 2.6. Supp ose w e are given a mo dular blac k b o x for the following unkno wn p olynomial: f ( x ) = x 15 − 45 x 14 + 945 x 13 − 12285 x 12 + 110565 x 11 − 729729 x 10 + 3648645 x 9 − 14073345 x 8 + 42220035 x 7 − 98513415 x 6 + 177324145 x 5 − 241805625 x 4 + 241805475 x 3 − 167403375 x 2 + 71743725 x − 14348421 , along with the b ounds B A = 4, B T = 2, B H = 4, and B N = 4. One ma y easily confirm that f ( x ) = ( x − 3) 15 − 2( x − 3) 5 , and hence these b ounds are actually tigh t. No w suppose the oracle pro duces p = 7 in Step 3 . W e use the black b o x to find f (0) , f (1) , . . . , f (6) in Z 7 , and dense interpolation to compute f (7) ( x ) = 5 x 5 + 2 x 4 + 3 x 3 + 6 x 2 + x + 4 . Since deg f (7) = 5 ≥ 2 B T + 1, we mov e on to Step 8 and compute each f (7) ( x + γ ) with γ = 1 , 2 , . . . , 6. Examining these, we see that f (7) ( x + 3) = 5 x 5 + x 3 has the few est non-zero and non-constant terms, and so set α 7 to 3 on Step 8 . This means the sparsest shift must b e congruen t to 3 mo dulo 7. This provides a single modular image for use in Chinese remaindering and ra- tional reconstruction on Step 15 , after enough successful iterations for different primes p . ♦ 2.3. Conditions for Success. W e hav e seen that, provided deg f > 2 B T , a go od prime p is one such that deg f ( p ) > 2 B T . The follo wing theorem provides (quite lo ose) sufficient conditions on p to satisfy this requiremen t. Theorem 2.7. Let f ∈ Q [ x ] as in (1.3) and B T ∈ N such that t ≤ B T . Then, for some prime p , the degree of f ( p ) is greater than 2 B T whenev er the following hold: In terp olation of Shifted-Lacunary Polynomials 9 ◦ c t 6≡ 0 mod p ; ◦ ∀ i ∈ { 1 , 2 , . . . , t − 1 } , e t 6≡ e i mo d ( p − 1) ; ◦ ∀ i ∈ { 1 , . . . , 2 B T } , e t 6≡ i mod ( p − 1) . Pr oof. The first condition guarantees that the last term of f ( p ) ( x ) as in (2.1) do es not v anish. W e also kno w that there is no other term with the same degree from the second condition. Finally , the third condition tells us that the degree of the last term will be greater than 2 B T . Hence the degree of f ( p ) is greater than 2 B T .  F or purp oses of computation it will b e con v enient to simplify the ab o ve conditions to t wo non-divisibilit y requirements, on p and p − 1 resp ectiv ely: Cor ollar y 2.8. Let f , B T , B H , B N b e as in the input to Algorithm 2.3 with deg f > 2 B T . Then there exist C 1 , C 2 ∈ N with log 2 C 1 ≤ 2 B H and log 2 C 2 ≤ B N (3 B T − 1) suc h that deg f ( p ) > 2 B T whenev er p - C 1 and ( p − 1) - C 2 . Pr oof. W rite f as in (1 . 3) . W e will use the sufficient conditions giv en in Theorem 2.7 . W rite | c t | = a/b for a, b ∈ N relatively prime. In order for c t rem p to b e well-defined and not zero, neither a nor b can v anish modulo p . This is true whenever p - ab . Set C 1 = ab . Since a, b ≤ 2 B H , log 2 C 1 = log 2 ( ab ) ≤ 2 B H . No w write C 2 = t − 1 Y i =1 ( e t − e i ) · 2 B T Y i =1 ( e t − i ) . W e can see that the second and third conditions of Theorem 2.7 are satisfied whenev er ( p − 1) - C 2 . Now, since eac h integer e i is distinct and positive, and e t is the greatest of these, eac h ( e t − e i ) is a p ositiv e in teger less than e t . Similarly , since e t = deg f > 2 B T , eac h ( e t − i ) in the second pro duct is also a p ositiv e in teger less than e t . Therefore, using the fact that t ≤ B T , we see C 2 ≤ e 3 B T − 1 t . F urthermore, e t ≤ 2 B N , so w e kno w that log 2 C 2 ≤ B N (3 B T − 1).  A similar criteria for success is required in Bl¨ aser et al. ( 2009 ), and they emplo y Linnik’s theorem to obtain a polynomial-time algorithm for p olynomial iden tity testing. Linnik’s theorem was also employ ed in Giesbrech t & Ro c he ( 2007 ) to yield a m uch more exp ensiv e deterministic p olynomial-time algorithm for finding sparse shifts than the one presented here. 10 Giesbrec ht & Ro c he 3. Interpolation Once w e kno w the v alue of the sparsest shift α of f , we can trivially construct a modular black b o x for the t -sparse polynomial f ( x + α ) using the mo dular blac k b o x for f . Therefore, for the purp oses of in terp olation, we can assume α = 0, and fo cus only on interpolating a t -sparse p olynomial f ∈ Q [ x ] giv en a mo dular blac k b o x for its ev aluation. The basic techniques of this section are, for the most part, kno wn in the literature. Ho wev er, a unified presen tation in terms of bit complexity for our mo del of mo dular black b o xes will b e helpful. F or conv enience, we restate the notation for f and f ( p ) , given a prime p : f = c 0 + c 1 x e 1 + c 2 x e 2 + · · · + c t x e t , (3.1) f ( p ) = ( c 0 rem p ) + ( c 1 rem p ) x e 1 rem 1 ( p − 1) + · · · + ( c t rem p ) x e t rem 1 ( p − 1) . (3.2) Again, w e assume that w e are giv en b ounds B H , B T , and B N on max i size( c i ), t , and log 2 deg f , resp ectiv ely . W e also in tro duce the notation τ ( f ), whic h is defined to b e the n um b er of distinct non-zero, non-constan t terms in the uni- v ariate p olynomial f . This algorithm will again use the p olynomials f ( p ) for primes p , but now rather than a degree condition, we need f ( p ) to ha ve the maximal n umber of non-constan t terms. So w e define a prime p to be a go o d prime for interp olation if and only if τ ( f ( p ) ) = t . Again, the term “go o d prime” refers to this kind of prime for the remainder of this section. No w supp ose we hav e used mo dular ev aluation and dense interpolation (as in Algorithm 2.3 ) to recov er the p olynomials f ( p ) for k distinct go o d primes p 1 , . . . , p k . W e therefore ha ve k images of eac h exp onen t e i mo dulo ( p 1 − 1) , . . . , ( p k − 1). W rite eac h of these p olynomials as: (3.3) f ( p i ) = c ( i ) 0 + c ( i ) 1 x e ( i ) 1 + · · · + c ( i ) t x e ( i ) t . Note that it is not generally the case that e ( i ) j = e j rem 1 ( p i − 1). Because w e don’t kno w how to asso ciate the exp onen ts in each p olynomial f ( p i ) with their pre-image in Z , a simple Chinese remaindering on the exp onen ts will not work. P ossible approac hes are provided b y Kaltofen ( 1988 ), Kaltofen et al. ( 1990 ) or Av enda ˜ no et al. ( 2006 ). Ho wev er, the most suitable approac h for our purp oses is the clev er technique of Garg & Sc host ( 2009 ), based on ideas of Grigoriev & Karpinski ( 1987 ). W e in terp olate the p olynomial (3.4) g ( z ) = ( z − e 1 )( z − e 2 ) · · · ( z − e t ) , whose co efficien ts are symmetric functions in the e i ’s. Given f ( p i ) , we hav e all the v alues of e ( i ) j rem 1 ( p i − 1) for j = 1 , . . . , t ; we just don’t kno w the order. In terp olation of Shifted-Lacunary Polynomials 11 But since g is not dep enden t on the order, we can compute g mo d ( p i − 1) for i = 1 , . . . , k , and then find the ro ots of g ∈ Z [ x ] to determine the exp onen ts e 1 , . . . , e t . Once w e know the exp onen ts, w e reco ver the co efficien ts from their images modulo eac h prime. The correct coefficient in eac h f ( p ) can be iden tified b ecause the residues of the exp onen ts mo dulo p − 1 are unique, for each c hosen prime p . This approac h is made explicit in the following algorithm. Algorithm 3.5. Sparse Polynomial In terp olation ov er Q [ x ]. Input: ◦ A mo dular black b o x for unkno wn f ∈ Q [ x ] ◦ Bounds B H and B N as describ ed ab o ve ◦ An oracle whic h pro duces primes and indicates when at least one go od prime must ha ve b een returned Output: f ∈ Q [ x ] as in (3.1) 1. Q ← 1, P ← 1, k ← 1, t ← 0 2. While log 2 P < 2 B H + 1 or log 2 Q < B N or the oracle do es not guaran tee a go o d prime has b een pro duced do 3–8 3. p k ← new prime from the oracle 4. Compute f ( p k ) b y blac k b o x calls and dense interpolation 5. If τ ( f ( p k ) ) > t then 6. Q ← p k − 1, P ← p k , t ← τ ( f ( p k ) ), p 1 ← p k , f ( p 1 ) ← f ( p k ) , k ← 2 7. Else if τ ( f ( p k ) ) = t then 8. Q ← lcm( Q, p k − 1), P ← P · p k , k ← k + 1 9. F or i ∈ { 1 , . . . , k − 1 } do 10. g ( p i ) ← Q 1 ≤ j ≤ t ( z − e ( i ) j ) mo d p i − 1 11. Construct g = a 0 + a 1 z + a 2 z 2 + · · · + a t z t ∈ Z [ x ] suc h that g ≡ g ( p i ) mo d p i − 1 for 1 ≤ i < k , b y Chinese remaindering 12. F actor g as ( z − e 1 )( z − e 2 ) · · · ( z − e t ) to determine e 1 , . . . , e t ∈ Z 13. F or 1 ≤ i ≤ t do 14. F or 1 ≤ j ≤ k do 15. Find the exp onen t e ( j ) ` j of f ( p j ) suc h that e ( j ) ` j ≡ e i mo d p j − 1 16. Reconstruct c i ∈ Q by Chinese remaindering from residues c (1) ` 1 , . . . , c ( k ) ` k 17. Reconstruct c 0 ∈ Q by Chinese remaindering from residues c (1) 0 , . . . , c ( k ) 0 The following theorem follows from the ab o ve discussion. Theorem 3.6. Algorithm 3.5 w orks correctly as stated. 12 Giesbrec ht & Ro c he Again, this algorithm runs in polynomial time in the b ounds giv en, but w e p ostpone the detailed complexity analysis until Section 5.2 , after w e discuss ho w to choose primes from the “oracle”. Some small practical improv ements ma y b e gained if we use Algorithm 3.5 to in terp olate f ( x + α ) after running Algorithm 2.3 to determine the sparsest shift α , since in this case we will hav e a few previously-computed p olynomials f ( p ) . Ho wev er, w e do not explicitly consider this savings in our analysis, as there is not necessarily any asymptotic gain. No w we just need to analyze the conditions for primes p to b e go od. This is quite similar to the analysis of the sparsest shift algorithm ab o ve, so we omit man y of the details here. Theorem 3.7. Let f , B T , B H , B N b e as ab o ve. There exist C 1 , C 2 ∈ N with log 2 C 1 ≤ 2 B H B T and log 2 C 2 ≤ 1 2 B N B T ( B T − 1) such that τ ( f ( p ) ) is maximal whenev er p - C 1 and ( p − 1) - C 2 . Pr oof. Let f b e as in (3.1) , write | c i | = a i /b i in lo w est terms for i = 1 , . . . , t , and define C 1 = t Y i =1 a i b i , C 2 = t Y i =1 t Y j = i +1 ( e j − e i ) . No w supp ose p is a prime suc h that p - C 1 and ( p − 1) - C 2 . F rom the first condition, w e see that eac h c i mo d p is w ell-defined and nonzero, and so none of the terms of f ( p ) v anish. F urthermore, from the second condition, e i 6≡ e k mo d p − 1 for all i 6 = j , so that none of the terms of f ( p ) collide. Therefore f ( p ) con tains exactly t non-constant terms. The b ounds on C 1 and C 2 follo w from the facts that each a i , b i ≤ 2 B H and eac h difference of exp onen ts is at most 2 B N .  4. Generating primes W e no w turn our atten tion to the problem of generating primes for the sparsest shift and interpolation algorithms. In previous sections we assumed we had an “oracle” for this, but no w w e present an explicit and analyzed algorithm. The definition of a “goo d prime” is not the same for the algorithms in Section 2 and Section 3 . Ho wev er, Corollary 2.8 and Theorem 3.7 pro vide a unified presentation of sufficien t conditions for primes b eing “go o d”. Here w e call a prime which satisfies those sufficient conditions a useful prime . So ev ery useful prime is go o d (with the b ounds appropriately sp ecified for the relev an t algorithm), but some go o d primes might not b e useful. In terp olation of Shifted-Lacunary Polynomials 13 W e first describ e a set P of primes such that the num b er and densit y of useful primes within the set is sufficiently high. W e will assume that there exist n umbers C 1 , C 2 , and useful primes p are those such that p - C 1 and ( p − 1) - C 2 . The num b ers C 1 and C 2 will b e unknown, but we will assume we are given b ounds β 1 , β 2 suc h that log 2 C 1 ≤ β 1 and log 2 C 2 ≤ β 2 . Supp ose w e w ant to find ` useful primes. W e construct P explicitly , of a size guaranteed to con tain enough useful primes, then en umerate it. The follo wing fact is immediate from Mik aw a ( 2001 ), though it has b een somewhat simplified here, and the use of (unknown) constan ts is made more explicit. This will b e imp ortan t in our computational metho ds. F or q ∈ Z , let S ( q ) b e the smallest prime p such that q | ( p − 1). F act 4.1 ( Mik a w a 2001 ). There exists a constan t µ > 0 , such that for all n > µ , and for all in tegers q ∈ { n, . . . , 2 n } with fewer than µn/ log 2 n exceptions, w e ha ve S ( q ) < q 1 . 89 . Our algorithms for generating useful primes require explicit knowledge of the v alue of the constant µ in order to run correctly . So w e will assume that w e know µ in what follows. T o get around the fact that w e do not, w e simply start by assuming that µ = 1, and run any algorithm dep ending up on it. If the algorithm fails w e simply double our estimate for µ and rep eat. A t most a constan t n umber of doublings is required. W e mak e no claim this is particularly practical. F or conv enience we define Υ( x ) = 3 x 5 log x − µx log 2 x . Theorem 4.2. Let log 2 C 1 ≤ β 1 , log 2 C 2 ≤ β 2 and ` b e as ab o ve. Let n b e the smallest in teger suc h that n > 21 , n > µ and Υ( n ) > β 1 + β 2 + ` . Define Q = { q prime : n ≤ q < 2 n and S ( q ) < q 1 . 89 } , P = { S ( q ) : q ∈ Q} . Then the num b er of primes in P is at least β 1 + β 2 + ` , and the num b er of useful primes in P , such that p - C 1 and ( p − 1) - C 2 , is at least ` . F or all p ∈ P w e ha ve p ∈ O (( β 1 + β 2 + ` ) 1 . 89 · log 1 . 89 ( β 1 + β 2 + ` )) . Pr oof. By Rosser & Schoenfeld ( 1962 ), the num b er of primes b et w een n and 2 n is at least 3 n/ (5 log n ) for n ≥ 21. Applying F act 4.1 , we see # Q ≥ 3 n/ (5 log n ) − µn/ log 2 n when n ≥ max { µ, 21 } . Now suppose S ( q 1 ) = S ( q 2 ) 14 Giesbrec ht & Ro c he for q 1 , q 2 ∈ Q . If q 1 < q 2 , then S ( q 1 ) > q 2 1 , a con tradiction with the definition of Q . So w e m ust hav e q 1 = q 2 , and hence # P = # Q ≥ Υ( n ) > β 1 + β 2 + `. W e kno w that there are at most log 2 C 1 ≤ β 2 primes p ∈ P suc h that p | C 1 . W e also kno w that there are at most log 2 C 2 ≤ β 2 primes q ∈ Q suc h that q | C 2 , and hence at most log 2 C 2 primes p ∈ P suc h that p = S ( q ) and q | ( p − 1) | C 1 . Th us, b y construction P con tains at most β 1 + β 2 primes that are not useful out of β 1 + β 2 + ` total primes. T o analyze the size of the primes in P , we note that to mak e Υ( n ) > β 1 + β 2 + ` , w e hav e n ∈ Θ(( β 1 + β 2 + ` ) · log ( β 1 + β 2 + ` )) and eac h q ∈ Q satisfies q ∈ O ( n ). Elemen ts of P will b e of magnitude at most (2 n ) 1 . 89 and hence p ∈ O (( β 1 + β 2 + ` ) 1 . 89 log 1 . 89 ( β 1 + β 2 + ` )).  Giv en β 1 , β 2 and ` as ab o ve (where log 2 C 1 ≤ β 1 and log 2 C 2 ≤ β 2 for unkno wn C 1 and C 2 ), we generate the primes in P as follo ws. Start b y assuming that µ = 1, and compute n as the smallest integer such that Υ( n ) > β 1 + β 2 + ` , n ≥ µ and n ≥ 21. List all primes betw een n and 2 n using a Siev e of Eratosthenes. F or each prime q b et ween n and 2 n , determine S ( q ), if it is less than q 1 . 89 , b y simply c hecking if k q + 1 is prime for k = 1 , 2 , . . . , b q 0 . 89 c . If we find a prime p = S ( q ) < q 1 . 89 , add p to P . This is repeated un til P contains β 1 + β 2 + ` primes. If we are unable to find this n umber of primes, w e hav e underestimated µ (since Theorem 4.2 guarantees their existence), so we double µ and restart the process. Obviously in practice w e w ould not redo primality tests already p erformed for smaller µ , so really no w ork need b e wasted. Theorem 4.3. F or log 2 C 1 ≤ β 1 , log 2 C 2 ≤ β 2 , ` , and n as in Theorem 4.2 , we can generate β 1 + β 2 + ` elements of P with O (( β 1 + β 2 + ` ) 2 · log 7+ o (1) ( β 1 + β 2 + ` )) bit op erations. A t least ` of the primes in P will b e useful. Pr oof. The metho d and correctness follows from the ab o ve discussion. The Siev e of Eratosthenes can be run with O ( n log log log n ) bit operations (see Kn uth ( 1981 ), Section 4.5.4), and returns O ( n/ log n ) primes q betw een n and 2 n . Eac h primalit y test of k q + 1 can b e done with (log n ) 6+ o (1) bit op erations ( Lenstra & P omerance 2005 ), so the total cost is O ( n 2 (log n ) 5+ o (1) ) bit op er- ations. Since n ∈ O (( β 1 + β 2 + ` ) · log( β 1 + β 2 + ` )) the stated complexit y follo ws.  The analysis of our metho ds will be significan tly improv ed when more is disco vered ab out the b eha vior of the least prime congruent to one mo dulo In terp olation of Shifted-Lacunary Polynomials 15 a given prime, which w e ha ve denoted S ( q ). An asymptotic lo wer b ound of S ( q ) ∈ Ω( q log 2 q ) is conjectured in Granville & P omerance ( 1990 ), and we ha ve emplo yed the upper b ound from Mik a wa ( 2001 ) of S ( q ) ∈ O ( q 1 . 89 ) (with exceptions). F rom our o wn brief computational search w e ha ve evidence that the conjectured low er b ound may w ell be an upp er b ound: for all primes q ≤ 2 32 , S ( q ) < 2 q ln 2 q . If something similar could b e pro ven to hold asymptotically (ev en with some exceptions), the complexit y results of this and the next section w ould b e impro v ed significantly . In any case, the actual cost of the algorithms discussed will b e a reflection of the true b eha vior of S ( q ), even before it is completely understo o d b y us. Ev en more improv ements might b e p ossible if this rather complicated con- struction is abandoned altogether, as useful primes w ould naturally seem to b e relativ ely plentiful. In particular, one would exp ect that if we randomly c ho ose primes p directly from a set which has, say , 4( β 1 + β 2 + ` ) primes, we migh t exp ect that the probabilit y that p | C 1 or ( p − 1) | C 2 to less than, say , 1 / 4. Pro ving this directly app ears to b e difficult. Perhaps most germane results to this are low er b ounds on the Carmichael Lambda function (whic h for the pro duct of distinct primes p 1 , . . . , p m is the LCM of p 1 − 1, . . . , p m − 1), which are to o weak for our purp oses. See Erd¨ os et al. ( 1991 ). 5. Complexity analysis W e are no w ready to giv e a formal complexit y analysis for the algorithms presen ted in Section 2 and Section 3 . F or all algorithms, the complexity is p olynomial in the four b ounds B A , B T , B H , and B N defined in Section 2.2 , and since these are each b ounded ab o v e by size( f ), our algorithms will hav e p olynomial complexity in the size of the output if these b ounds are sufficien tly tigh t. 5.1. Complexit y of Sparsest Shift Computation. Algorithm 2.3 gives our algorithm to compute the sparsest shift α of an unkno wn polynomial f ∈ Q [ x ] given b ounds B A , B T , B H , and B N and an oracle for c ho osing primes. The details of this oracle are giv en in Section 4 . T o c ho ose primes, w e set ` = 2 B A + 1, and β 1 = 2 B H and β 2 = B N (3 B T − 1) (according to Corollary 2.8 ). F or the sake of notational brevity , define B Σ = B A + B H + B N B T so that β 1 + β 2 + ` ∈ O ( B Σ ). Theorem 5.1. Supp ose f ∈ Q [ x ] is an unkno wn p olynomial giv en by a blac k b o x, with b ounds B A , B T , B H , and B N giv en as ab o v e. If deg f > 2 B T , then 16 Giesbrec ht & Ro c he the sparsest shift α ∈ Q of f can be computed deterministically using O  B A B 1 . 89 Σ · log 2 . 89 B Σ · M ( B 1 . 89 Σ log 1 . 89 B Σ ) · M (log B Σ )  bit op erations, plus O ( κ f B 2 . 89 Σ log 1 . 89 B Σ M (log B Σ )) bit op erations for the black- b o x ev aluations. Pr oof. Algorithm 2.3 will alwa ys terminate (b y satisfying the conditions of Step 2 ) after 2 B A + 1 go od primes hav e b een pro duced by the oracle. Using the oracle to c ho ose primes, and b ecause β 1 + β 2 + ` ∈ O ( B Σ ), O ( B 2 Σ log 7+ o (1) B Σ ) bit operations are used to compute all the primes on Step 3 , b y Theorem 4.3 . And b y Theorem 4.2 , eac h chosen p is b ounded by O ( B 1 . 89 Σ log 1 . 89 B Σ ). All blac k-b o x ev aluations are p erformed on Step 4 ; there are p ev aluations at eac h iteration, and O ( B Σ ) iterations, for a total cost of O ( κ f B Σ p · M (log p )) bit op erations. The stated complexit y b ound follows from the size of eac h prime p . Steps 10–14 are nev er executed when deg f > 2 B T . Step 15 is only executed once and nev er dominates the complexit y . Dense p olynomial interpolation o v er Z p is p erformed at most O ( B Σ ) times on Step 5 and O ( p ) times at eac h of O ( B A ) iterations through Step 7 . Since p  B Σ , the latter step dominates. U sing asymptotically fast metho ds, eac h in terp olation of f ( p ) ( x + γ ) uses O ( M ( p ) log p ) field op erations in Z p , eac h of whic h costs O ( M (log p )) bit op erations. This gives a total cost ov er all iterations of O ( B A p · log p · M ( p ) · M (log p )) (a slight abuse of notation here since the v alue of p v aries). Again, using the fact that p ∈ O ( B 1 . 89 Σ log 1 . 89 B Σ ) gives the stated result.  T o simplify the discussion somewhat, consider the case that we hav e only a single b ound on the size of the output p olynomial, sa y B f ≥ size( f ). By setting eac h of B T , B H , and B N equal to B f , and b y using the m ultiplication algorithm from Can tor & Kaltofen ( 1991 ), w e obtain the follo wing comprehensiv e result: Cor ollar y 5.2. The sparsest shift α of an unkno wn p olynomial f ∈ Q [ x ] , whose shifted-lacunary size is b ounded b y B f , can b e computed using O  B 8 . 56 f · log 6 . 78 B f · (loglog B f ) 2 · logloglog B f  bit op erations, plus O  κ f B 5 . 78 f · log 2 . 89 B f · loglog B f · logloglog B f  bit op erations for the black-box ev aluations. In terp olation of Shifted-Lacunary Polynomials 17 Pr oof. The stated complexities follow directly from Theorem 5.1 abov e, using the fact that M ( n ) ∈ O ( n log n loglog n ) and B Σ ∈ O ( B 2 f ). Using the single b ound B f , we see that these costs alw ays dominate the cost of Steps 10– 14 given in Theorem 2.5 , and so we ha ve the stated general result.  In fact, if w e ha ve no bounds at all a priori , w e could start b y setting B f to some small v alue (p erhaps dep enden t on the size of the blac k b o x or κ f ), running Algorithm 2.3 , then doubling B f and running the algorithm again, and so forth until the same p olynomial f is computed in successive iterations. This can then b e tested on random ev aluations. Such an approac h yields an output-sensitiv e p olynomial-time algorithm which should b e correct on most input, though it could certainly b e fo oled into early termination. This is a significan t improv emen t ov er the algorithms from our original pa- p er ( Giesbrech t & Ro c he 2007 ), whic h had a dominating factor of B 78 f in the deterministic complexity . Also — and somewhat surprisingly — our algorithm is competitive even with the b est-known sparsest shift algorithms which require a (dense) f ∈ Q [ x ] to be giv en explicitly as input. By carefully constructing the mo dular black b o x from a given f ∈ Q [ x ], and b eing sure to set B T < (deg f ) / 2, w e can derive from Algorithm 2.3 a deterministic sparsest-shift algorithm with bit complexit y close to the fastest algorithms in Giesbrec ht et al. ( 2003 ); the d e- p endence on degree n and sparsit y t will b e somewhat less, but the dep endence on the size of the co efficien ts log k f k is greater. T o understand the limits of our computational techniques (as opp osed to our curren t understanding of the least prime in arithmetic progressions) we consider the cost of our algorithms under the optimistic assumption that S ( q ) ∈ O ( q ln 2 q ), p ossibly with a small n umber of exceptions. In this case the sparsest shift α of an unknown p olynomial f ∈ Q [ x ], whose shifted-lacunary size is b ounded b y B f , can b e computed using O  B 5 f · log 6 B f · (loglog B f ) 2 · logloglog B f  bit op erations. As noted in the previous section, we ha ve v erified computation- ally that S ( q ) ≤ 2 q ln 2 q for q < 2 32 . This w ould suggest the ab o ve complexit y for all sparsest-shift interpolation problems that we w ould exp ect to encounter. 5.2. Complexit y of Interpolation. The complexit y analysis of the sparse in terp olation algorithm given in Algorithm 3.5 will b e quite similar to that of the sparsest shift algorithm ab ov e. Here, w e need ` = max { 2 B H + 1 , B N } go od primes to satisfy the conditions of Step 2 , and from Theorem 3.7 , we set β 1 = 2 B H B T and β 2 = 1 2 B N B T ( B T − 1). Hence for this subsection w e set B S = B T ( B H + B N B T ) so that β 1 + β 2 + ` ∈ O ( B S ). 18 Giesbrec ht & Ro c he Theorem 5.3. Supp ose f ∈ Q [ x ] is an unknown p olynomial given b y a mo d- ular blac k b o x, with b ounds B T , B H , B N , and B S giv en as ab o ve. The sparse represen tation of f as in (1.3) can b e computed with O  B S log B S · M ( B 1 . 89 S log 1 . 89 B S ) · M (log B S ) + B 2 N · M  ( B N + log B T ) log( B N + log B T )   bit op erations, plus O ( κ f B 2 . 89 S log 1 . 89 B S M (log B S )) bit op erations for the black- b o x ev aluations. Pr oof. As in the sparsest-shift computation, the cost of choosing primes in Step 3 requires O ( B 2 S log 7+ o (1) B S ) bit op erations, and each chosen prime p k satisfies p k ∈ O ( B 1 . 89 S log 1 . 89 B S ). The total cost ov er all iterations of Step 4 is also similar to b efore, O ( B S · M ( p k ) log p k · M (log p k )) bit op erations, plus O ( κ f B S p M (log p k )) for the black-box ev aluations. W e can compute each g ( p i ) in Step 10 using O ( M ( t ) log t ) ring op erations mo dulo p i − 1. Note that k ∈ O ( ` ), whic h is O ( B H + B N ), so the total cost in bit op erations for all iterations of this step is O (( B H + B N ) log B T · M ( B T ) · M (log B S )). Step 11 p erforms t Chinese Remainderings each of k mo dular images, and the size of each resulting integer is b ounded b y 2 B N , for a cost of O ( B T log B N · M ( B N )) bit op erations. T o factor g in Step 12 , w e can use Algorithm 14.17 of von zur Gathen & Gerhard ( 2003 ), whic h has a total cost in bit op erations of O  B 2 T · M ( B N + log B T ) + B N ( B N + log B T ) · M (( B N + log B T ) log( B N + log B T ))  b ecause the degree of g is t , g has t distinct ro ots, and eac h co efficien t is b ounded b y 2 B N . In Step 15 , w e m ust first compute the mo dular image of e i mo d p j − 1 and then lo ok through all t exp onen ts of f ( p j ) to find a match. This is rep eated tk times. W e can use fast modular reduction to compute all the images of each e i using O ( M ( B N ) log B N ) bit op erations, so the total cost is O ( B T ( B H B T + B N B T + M ( B N ) log B N )) bit op erations. Finally , we p erform Chinese remaindering and rational reconstruction of t + 1 rational num b ers, each of whose size is b ounded by B H , for a total cost of O ( B T · M ( B H ) log B H ). Therefore w e see that the complexity is dominated either b y the dense in terp olation in Step 4 or the ro ot-finding algorithm in Step 12 , dep ending essen tially on whether B N dominates the other b ounds.  In terp olation of Shifted-Lacunary Polynomials 19 Once again, b y ha ving only a single b ound on the size of the output, the complexit y measures are greatly simplified. Cor ollar y 5.4. Given a mo dular blac k b o x for an unknown p olynomial f ∈ Q [ x ] and a b ound B f on the size of its lacunary represenation, that repre- sen tation can b e interpolated using O  B 8 . 67 f log 4 . 89 B f (loglog B f ) 2 logloglog B f  bit op erations, plus O  κ f B 8 . 67 f log 2 . 89 B f loglog B f logloglog B f  bit op erations for the black-box ev aluations. Similar improv ements to those discussed at the end of Section Section 5.1 can b e obtained under stronger (but unpro ven) num b er theoretic assumptions. 6. Conclusions and F uture W ork Here w e provide the first algorithm to in terp olate an unkno wn univ ariate ratio- nal p olynomial into the sparsest shifted p o wer basis in time p olynomial in the size of the output. The main to ol w e hav e in tro duced is mapping do wn mo dulo small primes where the sparse shift is also mapp ed nicely . This technique could b e useful for other problems in v olving lacunary polynomials as w ell, although it is not clear how it w ould apply in finite domains where there is no notion of “size”. There are man y further a ven ues to consider, the first of whic h might b e m ultiv ariate p olynomials with a shift in eac h v ariable (see, e.g., Grigoriev & Lakshman ( 2000 )). It would b e easy to adapt our algorithms to this case pro vided that the degree in e ach variable is more than twice the sparsit y (this is called a “very sparse” shift in Giesbrech t et al. ( 2003 )). Finding multiv ariate shifts in the general case seems more difficult. Even more c hallenging would b e allo wing m ultiple shifts, for one or more v ariables — for example, finding sparse g 1 , . . . , g k ∈ Q [ x ] and shifts α 1 , . . . , α k ∈ Q such that the unkno wn p olynomial f ( x ) equals g 1 ( x − α 1 ) + · · · + g k ( x − α k ). The most general problem of this t yp e, whic h w e are v ery far from solving, migh t b e to compute a minimal-length form ula or minimal-size algebraic circuit for an unknown function. W e hop e that the small step tak en here migh t pro vide some insigh t to wards this ultimate goal. 20 Giesbrec ht & Ro c he Ac kno wledgemen t The authors w ould like to thank Igor Shparlinski for p oin ting out the pap er of Mik a wa ( 2001 ), and for suggesting ho w to discard “exceptional” primes q . This a v oids the use of Linnik’s theorem, as employ ed in Giesbrec h t & Ro c he ( 2007 ), and impro ves the complexity considerably . The authors w ould also like to thank Erich Kaltofen for discussions and sharing of his early unpublished w ork on rational in terp olation, and ´ Eric Sc host for discussions and sharing a pre-print of Garg & Schost ( 2009 ). Finally , the authors would lik e to thank the anon ymous review ers for their careful readings and useful suggestions. An extended abstract of a preliminary version of this work app eared at the MA CIS 2007 conference ( Giesbrec ht & Roche 2007 ). References Mar tn A venda ˜ no, Teresa Kri ck, and Ariel P acetti , Newton-hensel interpo- lation lifting. F oundations of Computational Mathematics (2006), 81–120. Michael Ben-Or and Prasoon Tiw ari , A deterministic algorithm for sparse m ultiv ariate polynomial in terp olation. In STOC ’88: Pro ceedings of the tw entieth ann ual A CM symp osium on Theory of computing , New Y ork, NY, USA, 1988, ACM Press, 301–309. Markus Bl ¨ aser, Moritz Hardt, Richard J. Lipton, and Nisheeth K. Vish- noi , Deterministically testing sparse polynomial iden tities of unbounded degree. Inf. Pro cess. Lett. 109 (3) (2009), 187–192. Allan Borodin and Ian Munro , The computational complexity of algebraic and numeric problems . American Elsevier Publishing Co., Inc., New Y ork-London- Amsterdam, 1975. Elsevier Computer Science Library; Theory of Computation Se- ries, No. 1. Allan Borodin and Prasoon Tiw ari , On the decidabilit y of sparse univ ariate p olynomial in terp olation. Comput. Complexity 1 (1) (1991), 67–90. D a vid Cantor and Erich Kal tofen , F ast multiplication of p olynomials ov er arbitrary algebras. Acta Informatica 28 (1991), 693–701. P aul Erd ¨ os, Carl Pomerance, and Eric Schmutz , Carmichael’s lam b da func- tion. Acta Arithmetica 58 (1991), 363–385. Sanchit Gar g and ´ Eric Schost , In terp olation of p olynomials giv en b y straigh t- line programs. Theoretical Computer Science 410 (27–29) (2009), 2659–2662. In terp olation of Shifted-Lacunary Polynomials 21 Joa chim von zur Ga then and J ¨ urgen Gerhard , Modern computer algebra . Cam bridge Universit y Press, Cambridge, second edition, 2003. Mark Giesbrecht and Daniel R oche , In terp olation of shifted-lacunary p oly- nomials. In Pro c. Mathematical Asp ects of Computer and Information Sciences (MA CIS’07) , Paris, F rance, 2007. Mark Giesbrecht, Erich Kal tofen, and Wen-shin Lee , Algorithms for com- puting sparsest shifts of p olynomials in p o w er, Chebyshev and P o c hhammer bases. J. Sym b olic Comput. 36 (3-4) (2003), 401–424. International Symp osium on Symbolic and Algebraic Computation (ISSA C’2002) (Lille). Mark Giesbrecht, George Labahn, and Wen-shin Lee , Sym b olic-n umeric sparse in terp olation of multiv ariate p olynomials. In Pro c. ACM In ternational Sym- p osium on Sym b olic and Algebraic Computation (ISSA C) , 2006, 116–123. Andrew Granville and Carl Pomerance , On the Least Prime in Certain Arithmetic Progressions. J. London Math. Soc. s2-41 (2) (1990), 193–200. Dima Grigoriev and Marek Karpinski , The matching problem for bipartite graphs with p olynomially b ounded p ermanen ts is in NC. In F oundations of Computer Science (FOCS) , 1987, 166–172. Dima Grigoriev and Marek Karpinski , A zero-test and an interpolation algo- rithm for the shifted sparse polynomials. In Applied algebra, algebraic algorithms and error-correcting co des (San Juan, PR, 1993) , v ol. 673 of Lecture Notes in Com- put. Sci. , 162–169. Springer, Berlin, 1993. Dima Grigoriev and Y. Lakshman , Algorithms for computing sparse shifts for m ultiv ariate p olynomials. Appl. Algebra Engrg. Comm. Comput. 11 (1) (2000), 43– 67. Erich Kal tofen , Notes on polynomial and rational function in terp olation. (1988). Unpublished manuscript. Erich Kal tofen and Wen-shin Lee , Early termination in sparse in terp olation algorithms. J. Symbolic Comput. 36 (3-4) (2003), 365–400. In ternational Symposium on Symbolic and Algebraic Computation (ISSAC’2002) (Lille). E. Kal tofen, Y. N. Lakshman, and J.-M. Wiley , Mo dular rational sparse m ultiv ariate p olynomial interpolation. In ISSA C ’90: Proceedings of the in ternational symp osium on Symbolic and algebraic computation , New Y ork, NY, USA, 1990, A CM Press, 135–139. 22 Giesbrec ht & Ro c he Donald E. Knuth , The Art of Computer Programming, V ol.2, Seminumerical Algorithms . Addison-W esley , Reading MA, 2 edition, 1981. Y. N. Lakshman and B. Da vid Saunders , Sparse shifts for univ ariate p olynomi- als. Appl. Algebra Engrg. Comm. Comput. 7 (5) (1996), 351–364. Hendrik W. Lenstra, Jr. and Carl Pomerance , Primalit y testing with Gaus- sian p erio ds. Preprint, 2005. Hiroshi Mika w a , On primes in arithmetic progressions. Tsukuba journal of math- ematics 25 (1) (2001), 121–153. J. Barkley R osser and Lo well Schoenfeld , Appro ximate form ulas for some functions of prime n um b ers. Illinois J. Math. 6 (1962), 64–94. Man uscript received June 9, 2022 Mark Giesbrecht Sc ho ol of Computer Science Univ ersity of W aterlo o W aterlo o, ON, N2L 3G1, Canada mwg@cs.uwaterloo.ca http://www.cs.uwaterloo.ca/ ~ mwg D aniel S. Roche Sc ho ol of Computer Science Univ ersity of W aterlo o W aterlo o, ON, N2L 3G1, Canada droche@cs.uwaterloo.ca http://www.cs.uwaterloo.ca/ ~ droche