Faster Interpolation Algorithms for Sparse Multivariate Polynomials Given by Straight-Line Programs\
FFaster Interpolation Algorithms for Sparse MultivariatePolynomials Given by Straight-Line Programs ∗ Qiao-Long Huang , and Xiao-Shan Gao , KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences University of Chinese Academy of Sciences
Abstract
In this paper, we propose new deterministic and Monte Carlo interpolation algorithms for sparsemultivariate polynomials represented by straight-line programs. Let f be an n -variate polyno-mial given by a straight-line program, which has a degree bound D and a term bound T . Ourdeterministic algorithm is quadratic in n, T and cubic in log D in the Soft-Oh sense, whichhas better complexities than existing deterministic interpolation algorithms in most cases. OurMonte Carlo interpolation algorithms have better complexities than existing Monte Carlo in-terpolation algorithms and are the first algorithms whose complexities are linear in nT in theSoft-Oh sense. Since nT is a factor of the size of f , our Monte Carlo algorithms are optimal in n and T in the Soft-Oh sense. The sparse interpolation for multivariate polynomials has received considerable interest. There aretwo basic models for this problem: the polynomial is either given as a straight-line program (SLP)[7, 11, 12, 14] or a more general black-box [1, 8, 13, 16, 21]. In this paper, we consider the problemof interpolation for a sparse multivariate polynomial given by an SLP.
Let f ∈ R [ x , . . . , x n ] be an SLP polynomial of length L , with a degree bound D and a term bound T , where R is a computable ring. In this paper, all complexity analysis relies on the “Soft-Oh”notation O ∼ ( φ ) = O ( φ · polylog ( φ )), where polylog means log c for some fixed c >
0. When wesay “linear”, “optimal”, etc., we mean linear and optimal in the sense of “Soft-Oh” complexity.We propose a new deterministic interpolation algorithm for f given by an SLP, whose complexityis O ∼ ( Ln T log D + LnT log D ) R arithmetic operations and a similar number of extra bit opera-tions. We also propose two new Monte Carlo interpolation algorithms for SLP multivariate polyno-mials. For a given µ ∈ (0 , O ∼ ( LnT (log D + log D log µ )) R arithmetic operations, and with probability at least 1 − µ , it returns the correct polynomial. Inour second algorithm, R is a finite field F q and we can evaluate f in a proper extension field of F q . The bit complexity of our second algorithm is O ∼ ( LnT log D (log D + log q ) log µ ), and withprobability at least 1 − µ , it returns the correct polynomial. These algorithms are the first oneswhose complexity is linear and optimal in n and T . ∗ Partially supported by a NSFC grant No.11688101. a r X i v : . [ c s . S C ] J u l lgorithms Total Cost TypeDense LD n DeterministicGarg & Schost [11] Ln T log D DeterministicRandomized G & S [12] Ln T log D Las VegasArnold, Giesbrecht & Roche [4] Ln T log D Monte CarloThis paper (Thm 5.9) Ln T log D + LnT log D DeterministicThis paper (Thm 6.8)
LnT log D Monte Carlo
Table 1: A “soft-Oh” comparison for SLP polynomials over an arbitrary ring R In Table 1, we list the complexities for the SLP interpolation algorithms which are similar tothe methods proposed in this paper. In the table, L is the size of the SLP and “total cost” meansthe number of arithmetic operations in R . For deterministic algorithms, the ratio of the cost ofour deterministic method with that of the algorithm given in [11] is nT +log DnT . So our method hasbetter complexities unless f is super sparse, more precisely, unless T < (cid:113) log Dn in the Soft-Ohsense. Noting that a dense polynomial with degree d has (cid:18) dd + n (cid:19) terms, our algorithm worksbetter in most cases. For probabilistic algorithms, our method is the only one whose complexity islinear in nT and is better than that of the Monte Carlo method given in [4].Kaltofen gave an interpolation algorithm for SLP polynomials [14], whose complexity is poly-nomial in D . Avenda˜no-Krick-Pacetti [7] gave an algorithm for interpolating an SLP f ∈ Z [ x ],with bit complexity polynomial in L, log D, h and h (cid:48) , where h is an upper bound on the height of f and h (cid:48) an upper bound on the height of the values of f (and some of its derivatives) at somesample points. This method does not seem to extend to arbitrary rings. Mansour [19] and Alon-Mansour [1] gave deterministic algorithms for polynomials in Z [ x ], with bit complexity polynomialin n, log D, T, H , where H is an upper bound on the bit-length of the output coefficients. Also,this method is hard to extend to arbitrary rings. Note that interpolation algorithms for black-boxpolynomials [8, 16, 21, 17] can also be used to SLP polynomials, but the complexities of thesealgorithms are polynomial in D instead of log D . Bit Field AlgorithmComplexity Extension typeGarg & Schost [11] Ln T log D log q not DeterministicRandomized Garg-Schost [12] Ln T log D log q not Las VegasGiesbrecht & Roche [12] Ln T log D ( n log D + log q ) yes Las VegasArnold, Giesbrecht & Roche [4] Ln T log D log q not Monte CarloArnold, Giesbrecht & Roche [5] LnT log D (log D + log q ) + n ω T yes Monte CarloArnold, Giesbrecht & Roche [6] Ln log D ( T log D + n )(log D + log q ) yes Monte Carlo+ n ω − T log D + n ω log D This paper (Thm 5.9) Ln T log D log q + LnT log D log q not DeterministicThis paper (Thm 6.8) LnT log D log q not Monte CarloThis paper (Thm 6.11) LnT log D (log q + log D ) yes Monte Carlo Table 2: A “soft-Oh” comparison for SLP polynomials over finite field F q In Table 2, we list the bit complexities for SLP interpolation algorithms over the finite field F q . In the table, “field extension” means that in the probe of the SLP, whether elements in anextension field of F q are needed. It is easy to see that our probabilistic algorithms are the only oneswhich are linear in n and T . Also, our second probabilistic algorithm has better complexity thanall algorithms except randomized Garg-Schost [12] which is Las Vegas. The ratio of the cost of our2rst probabilistic algorithm with that of the randomized Garg-Schost [12] is log DnT , so our methodis faster in most cases. Our methods build on the work by Garg-Schost [11], Arnold-Giesbrecht-Roche [12], Giesbrecht-Roche [13], and Klivans-Spielman [17], where the basic idea is to reduce multivariate interpolation tounivariate interpolation. Three new techniques are introduced in this paper: a criterion for checkingwhether a term belongs to a polynomial (see Section 2 for details), a deterministic method to findan “ok” prime (see Section 3 for details), a new Kronecker type substitution to reduce multivariatepolynomial interpolation to univariate polynomial interpolation (see Section 5.1 for details). Ourmethods have three major steps: • First, we find an “ok” prime p such that at least half of the terms of f do not collide or mergewith other terms of f in the univariate polynomial f mod ( D,p ) = f ( x, x D , . . . , x D n − ) mod ( x p − • Second, we obtain a set S of terms containing those non-colliding terms of f . In the univariatecase, these terms are found by the Chinese Remaindering Theorem and in the multivariatecase, these terms are found by a new Kronecker substitution. • Finally, we use our criterion for checking whether a term belongs to a polynomial to find atleast half of the terms of f from S .Repeating these three steps for at most log T times, we obtain f . In the rest of this section, wegive detailed comparison with related work.Grag and Schost [11] gave a deterministic interpolation algorithm for a univariate SLP polyno-mial f by recovering f from f mod ( x p −
1) for O ( T log D ) different primes p . The randomizedLas Vegas version of this method needs O ( T log D ) probes. The multivariate interpolation comesdirectly from the Kronecker substitution [18]. Our univariate interpolation algorithm has two ma-jor differences from that given in [11]. First, we compute f mod ( x p −
1) for O ( T log D ) differentprimes p , and second we introduce a criterion to check whether a term really belongs to f . Ourmultivariate interpolation method is similar to our univariate interpolation algorithm, where a newKronecker type substitution is introduced to recover the exponents.Giesbrecht-Roche [12] introduced the idea of diversification and a probabilistic method to choose“good” primes. It improves Grag and Schost’s algorithm by a factor O ( T ), but becomes a LasVegas algorithm. In our Monte Carlo algorithms, we use a new Kronecker substitution insteadof the method of diversification to find the same term in different remainders of f , where onlyadditions are used for the coefficients. Hence, our algorithm can work for more general rings andhas better complexity.In Arnold, Giesbrecht, and Roche [4], the concept of “ok” prime is introduced and a MonteCarlo univariate algorithm is given, which has complexity linear in T but cubic in n . The “ok”prime in [4] is probabilistic. In our deterministic method, we give a method to find an exact “ok”prime.In Arnold, Giesbrecht, and Roche [5], their univariate interpolation algorithm is extended tofinite fields. By combining the idea of diversification, the complexity becomes better. This algo-rithm will be used in our second probabilistic algorithm. The bit complexity of their multivariateinterpolation algorithm is linear in n ω , where ω is the constant of matrix multiplication, while our3lgorithm is linear in n . The reason is that, our method uses a new Kronecker substitution to findthe exponents and does not need to solve linear systems.In Arnold, Giesbrecht, and Roche [6], they further improved their interpolation algorithm forfinite fields. By combining the random Kronecker substitution and diversification, the complexitybecomes better, but still linear in n ω .Finally, the new Kroneceker type substitution introduced in this paper is inspired by the works ofKlivans-Spielman [17], Arnold [2] and Arnold-Roche [3]. In [17], the substitution f ( q x, q x mod ( D,p ) ,. . . , q n x mod ( D n − ,p ) ) was used, where q , . . . , q n are primes. In this paper, we introduced the substi-tution f ( x, x mod ( D,p ) , . . . , x p + mod ( D k − ,p ) , . . . , x mod ( D n − ,p ) ) (see section 5.1 for exact definition).Our substitution has the following advantages: (1) For the complex filed, the size of data is notchanged after our substitution, while the size of data for the substitution in [17] is increased by afactor of D . (2) Only arithmetic operations for the coefficients are used in our algorithm and thusthe algorithm works for general computable rings, while the substitution in [17] needs factoriza-tion and R should be a UFD at least. In [2], the substitution f ( x s , x s , x s k + p . . . , x s n ) was used,where s i are random integers. Comparing to the randomized Kronecker substitution in [2, 3], oursubstitution is deterministic. In this section, we give a criterion to check whether a term belongs to a polynomial.Throughout this paper, let f = c m + c m + · · · + c t m t ∈ R [ X ] be a multivariate polynomialwith terms c i m i , where R is a computable ring, X = { x , x , . . . , x n } are n indeterminates, and m i , i = 1 , , . . . , t are distinct monomials. Denote f = t to be the number of terms of f and M f = { c m , c m , . . . , c t m t } to be the set of terms of f . Let D, T ∈ N such that D > deg( f ) and T ≥ f . For p ∈ N > , let f mod ( D,p ) = f ( x, x D , . . . , x D n − ) mod ( x p − ∈ R [ x ] . (1)We have the following key concept. Definition 2.1
A term cm ∈ M f is called a collision in f mod ( D,p ) if there exists an aw ∈ M f \{ cm } such that m mod ( D,p ) = w mod ( D,p ) . The following fact is obvious.
Lemma 2.2
Let f ∈ R [ x ] , deg( f ) < D and cm ∈ M f . If cm is not a collision in f mod ( D,p ) , then forany prime q , cm is also not a collision in f mod ( D,pq ) . Lemma 2.3
Let f = (cid:80) ti =1 c i m i , T ≥ f, D > deg f , N = max { , (cid:100) n ( T −
1) log D (cid:101)} . For each cm ∈ M f , there exist at most N − primes p , . . . , p N − such that cm is a collision in f mod ( D,p i ) forall i = 1 , . . . , N − .Proof. If T = 1, then N = 1. The lemma is obvious. Now we assume T ≥
2, then N = (cid:100) n ( T −
1) log D (cid:101) . It suffices to show that for any N different primes p , p , . . . , p N , there exists atleast one p j , such that cm is not a collision in f mod ( D,p j ) .4ssume m i = x e i, x e i, · · · x e i,n n , i = 1 , , . . . , t . We prove it by contradiction. It suffices toconsider the case of c m . We assume by contradiction that for every p j , j = 1 , , . . . , N , c m isa collision in f mod ( D,p j ) . Let B = t (cid:89) s =2 ( n (cid:88) i =1 e ,i D i − − n (cid:88) i =1 e s,i D i − ) . First, we show that if c m is a collision in f mod ( D,p j ) , then mod ( B, p j ) = 0. Since c m is a collisionin f mod ( D,p j ) , without loss of generality, assume m mod D,p j ) = m mod D,p j ) . Then 0 = deg( m mod D,p j ) ) − deg( m mod D,p j ) ) = mod ( (cid:80) ni =1 e ,i D i − , p j ) − mod ( (cid:80) ni =1 e ,i D i − , p j ). So we have mod ( (cid:80) ni =1 e ,i D i − − (cid:80) ni =1 e ,i D i − , p j ) = 0. So mod ( B, p j ) = 0.Since p , p , . . . , p N are different primes, (cid:81) N j =1 p j divides B . Note that | (cid:80) ni =1 e ,i D i − − (cid:80) ni =1 e s,i D i − | ≤ ( D − (cid:80) ni =1 D i − ) = D n −
1. So | B | = | (cid:81) ts =2 ( (cid:80) ni =1 e ,i D i − − (cid:80) ni =1 e s,i D i − ) | ≤ ( D n − t − . Thus (cid:81) N j =1 p j ≥ N ≥ n ( T −
1) log D = D n ( T − > ( D n − T − ≥ | B | , which contra-dicts the fact that (cid:81) N j =1 p j divides B . The lemma is proved.Now we give a criterion for testing whether a term cm is in M f . Theorem 2.4
Let f = (cid:80) ti =1 c i m i , T ≥ f, D > deg f , N = max { , (cid:100) n ( T −
1) log D (cid:101)} , N = (cid:100) nT log D (cid:101) , and P = { p , p , . . . , p N + N − } be N + N − different primes. For a term cm satisfying deg( m ) < D , cm ∈ M f if and only if there exist at least N integers j ∈ [1 , N + N − such that f − cm ) mod ( D,p j ) < f mod ( D,p j ) .Proof. If T = 1, then N = 1, the proof is obvious. So we assume T ≥
2, then N = (cid:100) n ( T −
1) log D (cid:101) .Let cm ∈ M f . If p j is a prime such that cm is not a collision in f mod ( D,p j ) , then f − cm ) mod ( D,p j ) = f mod ( D,p j ) −
1. So f − cm ) mod ( D,p j ) < f mod ( D,p j ) . By Lemma 2.3, there exist at most N − q j such that cm is a collision in f mod ( D,q j ) . In P , as N + N − − ( N −
1) = N , there exist at least N primes such that f − cm ) mod ( D,p j ) < f mod ( D,p j ) .For the other direction, assume cm / ∈ M f . We show there exist at most N − j ∈ [1 , N + N −
1] such that f − cm ) mod ( D,p j ) < f mod ( D,p j ) . Consider two cases: Case 1: m is not amonomial in f . Case 2: m is a monomial in f , but cm is not a term in f .Case 1. Since m is not a monomial in f , − cm is a term in M f − cm and f − cm ) ≤ T + 1. ByLemma 2.3, there exist at most N − P such that − cm is a collision in f − cm . For allother primes p j in P , f − cm ) − ( − cm )) mod ( D,p j ) = f − cm ) mod ( D,p j ) −
1, that is, f − cm ) mod ( D,p j ) = f mod ( D,p j ) + 1. So there exist at most N − p j in P such that f − cm ) mod ( D,p j ) < f mod ( D,p j ) .Case 2. Since m is a monomial in f and cm / ∈ M f , f − cm has the same number of terms as f .Assume the term of f with monomial m is c m . Then ( c − c ) m ∈ M f − cm . By Lemma 2.3, for atmost N − ≤ N −
1) primes p j in P , ( c − c ) m is a collision in ( f − cm ) mod ( D,p j ) . For all other primes p k in P , ( c − c ) m is not a collision in ( f − cm ) mod ( D,p k ) , or equivalently, f − cm ) mod ( D,p k ) = f − cm − ( c − c ) m ) mod ( D,p k ) + 1 = f − c m ) mod ( D,p k ) + 1. But we always have f − c m ) mod ( D,p k ) + 1 ≥ f mod ( D,p k ) ,and so f − cm ) mod ( D,p k ) ≥ f ) mod ( D,p k ) . So there exist at most N − p j in P such that f − cm ) mod ( D,p j ) < f mod ( D,p j ) . The theorem is proved.As a corollary, we can deterministically recover f from f mod ( D,p j ) , j = 1 , , . . . , N + N − orollary 2.5 Use the notations in Theorem 2.4. We can uniquely recover f from f mod ( D,p j ) , j =1 , , . . . , N + N − .Proof. Let C = { c , c , . . . , c k } be all the different coefficients in f mod ( D,p j ) , j = 1 , , . . . , N + N − N + N − ≥ N , all the coefficients of f are in C . Let M = { m , m , . . . , m s } be the set of all the monomials with degrees less than D . So all the terms of f are in { c i m j | i =1 , , . . . , k, j = 1 , , . . . , s } . By Theorem 2.4, we can check if c i m j is in M f . So we can find all theterms of f .The above result can be changed into a deterministic algorithm for interpolating f . But thealgorithm is not efficient due to the reason that s is linear in D n . In the following, we will showhow to find a smaller alternative set M and give an efficient interpolation algorithm. A prime p is called an “ok” prime if at least half of the terms in f are not collisions in f mod ( D,p ) , where D > deg f . In this section, we give a deterministic method to find an “ok” prime for f .Denote C f ( D,p ) to be the number of collision terms of f in f mod ( D,p ) . We need the following lemmafrom [4]. Lemma 3.1 [4] Let f ∈ R [ x ] . If f mod ( D,q ) ≤ f mod ( D,p ) , then C f ( D,p ) ≤ C f ( D,q ) . It is easy to modify the above lemma into multivariate case.
Corollary 3.2
Let f ∈ R [ X ] , D > deg( f ) . If f mod ( D,q ) ≤ f mod ( D,p ) , then C f ( D,p ) ≤ C f ( D,q ) .Proof. Let F = f ( x, x D , . . . , x D n − ), then f mod ( D,p ) = F mod ( x p − C f ( D,p ) = C F ( D,p ) . Thecorollary follows from Lemma 3.1. Lemma 3.3
Let f = (cid:80) ti =1 c i m i ∈ R [ X ] , m i = x e i, x e i, · · · x e i,n n , D > deg( f ) , A = (cid:81) i
3, then
K > s . So K ≥ (cid:100) s (cid:101) . We haveproved the lemma. Theorem 3.4
Let f = (cid:80) ti =1 c i m i ∈ R [ X ] , T ≥ f, D > deg f , N = max { , (cid:100) n ( T −
1) log D (cid:101)} ,and p , p , . . . , p N be N different primes. Let j be an integer in [1 , N ] such that f mod ( D,p j ) ≥ f mod ( D,p j ) for all j . Then at least (cid:100) t (cid:101) of the terms of f are not collisions in f mod ( D,p j ) . roof. If T = 1, then N = 1, the proof is obvious. So we assume T ≥
2, then N = (cid:100) n ( T −
1) log D (cid:101) .We first claim that there exists at least one p j in p , p , . . . , p N such that C f ( D,p j ) < t . We proveit by contradiction. Assume for j = 1 , , . . . , N , C f ( D,p j ) ≥ t . Then by Lemma 3.3, for all j in[1 , N ], p (cid:100) C f ( D,pj ) (cid:101) j divides A , where A is defined in Lemma 3.3. Since p j , j = 1 , , . . . , N aredifferent primes, then (cid:81) N j =1 p (cid:100) C f ( D,pj ) (cid:101) j divides A . Now (cid:81) N j =1 p (cid:100) C f ( D,pj ) (cid:101) j ≥ (cid:81) N j =1 p (cid:100) · t (cid:101) j ≥ N t ≥ tN ≥ tn ( T −
1) log D = D nt ( T − , which contradicts to the inequality A ≤ ( D n − t ( t − . Weproved the claim.By Corollary 3.2, we have C f ( D,p j ) ≤ C f ( D,p j ) , j = 1 , , . . . , N . So C f ( D,p j ) < · t = t . So thenumber of no collision terms of f in f mod ( D,p j ) is > t − t = t . We proved the theorem. In this section, we consider the interpolation of a univariate polynomial f with deg f < D . Thealgorithm works as follows. First, we use Theorem 3.4 to find an “ok” prime p such that at leasthalf of the terms of f are not collisions in f mod ( D,p ) . Second, we use f mod ( D,pp k ) , k = 1 , , . . . , K D to finda set S containing these non-collision terms of f by the Chinese Remaindering Theorem, where p k is the k -th prime and K D is the smallest number such that p p . . . p K D ≥ D . Finally, we useTheorem 2.4 to pick up the terms of f from S . x p − In this section, let f be a univariate polynomial in R [ x ]. We will give an algorithm to recover thoseterms of f from f mod ( D,p ) , which are not collisions in f mod ( D,p ) .Let f be a univariate polynomial, D > deg f , and p ∈ N >
0. In this case, f mod ( D,p ) = f ( x ) mod ( x p − f mod ( D,p ) = a x d + a x d + · · · + a r x d r (2) f mod ( D,pp k ) = f k, + f k, + · · · + f k,r + g k , k = 1 , , . . . , K D where d < d < · · · < d r , p k is the k -th prime, K D is the smallest number such that p p . . . p K D ≥ D , and f k,i mod ( x p −
1) = a i x d i and g k mod ( x p −
1) = 0. f mod ( D,pp k ) can be written as the aboveform, because f mod ( D,pp k ) mod ( x p −
1) = f mod ( D,p ) . We now introduce the following key notation U fD,p = { a i x e i | such that i ∈ [1 , r ] , a i is from (2), e i ∈ [0 , D − , andU1 : f k,i = a i x b k,i , k = 1 , , . . . , K D . (3)U2 : mod ( e i , p k ) = mod ( b k,i , p k ) , k = 1 , , . . . , K D . } The following lemma gives the geometric meaning of U fD,p . Lemma 4.1
Let f ∈ R [ x ] , deg( f ) < D and cm ∈ M f . If cm is not a collision in f mod ( D,p ) , then cm ∈ U fD,p . roof. It suffices to show that cm satisfies the conditions of the definition of U fD,p . Assume m = x e . Since cm is not a collision in f mod ( D,p ) , without loss of generality, assume cm mod ( D,p ) = a x d and d = mod ( e, p ), where a x d is defined in (2). By Lemma 2.2, cm is also not a collisionin f mod ( D,pp k ) and f k, = a x b k, for b k, = mod ( e, pp k ), k = 1 , , . . . , K D . Since mod ( e, p k ) = mod ( mod ( e, pp k ) , p k ) = mod ( b k, , p k ), conditions U1 and U2 are satisfied and the lemma isproved.Note that U fD,p may contain terms not in M f . The following algorithm computes the set U fD,p . Algorithm 4.2 (UTerms)Input:
Univariate polynomials f mod ( D,p ) , f mod ( D,pp k ) , k = 1 , , . . . , K D , a prime p , a degree bound D > deg( f ). Output: U fD,p . Step 1:
Write f mod ( D,p ) and f mod ( D,pp k ) in as form (2): f mod ( D,p ) = a x d + a x d + · · · + a r x d r f mod ( D,pp k ) = f k, + f k, + · · · + f k,r + g k where k = 1 , . . . , K D . Step 2:
Let U = {} . For i = 1 , , . . . , r a: If for k = 1 , , . . . , K D , one of f k,i ) (cid:54) = 1 or one of the coefficient of f k,i is not a i , then break. b: for k = 1 , , . . . , K D , assume f k,i = a i x b k,i and let e k = mod ( b k,i , p k ). c: Let β = Chinese Remainder ([ e , e , . . . , e K D ] , [ p , p , . . . , p K D ]). d: If β < D then let U = U (cid:83) { a i x β } . Step 3:
Return U . Remark 4.3 In c of Step 2, Chinese Remainder ([ e , e , . . . , e K D ] , [ p , p , . . . , p K D ]) means tofind an integer ≤ w < D such that mod ( w, p i ) = e i , i = 1 , , . . . , K D . Since p p · · · p K D ≥ D ,the integer w is unique. This can be done by the Chinese remainder algorithm. Lemma 4.4
Algorithm 4.2 needs O ( T log D ) arithmetic operations in R and O ∼ ( T log D log p + T log D ) bit operations.Proof. In Step 1, we need to do a traversal for the terms of f mod ( D,pp k ) . Since K D is the smallestnumber such that p p · · · p K D ≥ D , then 2 K D − ≤ p p · · · p K D − < D . So K D ≤ log D + 1,which is O (log D ). In order to write f mod ( D,pp k ) as the form (2), we needs to perform the modularoperation mod p on every degree of f mod ( D,pp k ) , then use the quick sorting method to write theirterms in ascending order according to the degree. Since f mod ( D,pp k ) has no more than T terms and k = 1 , , . . . , K D , it needs O ∼ ( T log D ) arithmetic operations. Since the height of the data is O (log( pp K D )) and the prime p K D is O ∼ (log D ), it needs O ∼ ( T log D log( pp K D )) bit operations,which is O ∼ ( T log D log p + T log D ) bit operations.8n a of Step 2, since f mod ( D,pp k ) ≤ T , it totally needs O ∼ ( T log D ) bit operations to determinewhether f k,i ) = 1. To compare the coefficients of f k,i , it needs O ( T log D ) arithmetic operationsin R . In b of Step 2, since the height of the data is O (log( pp K D )), it needs O ( T log D + T log D log p )bit operations. In c , we need to call at most T times Chinese remaindering. By [20, p.290], thecost of the Chinese remaindering algorithm is O ∼ (log D ) arithmetic operations in Z . Since theheight of the data is O (log D ), it needs O ∼ ( T log D ) bit operations. So the complexity of Step 2is O ∼ ( T log D + T log D log p ) bit operations and O ( T log D ) arithmetic operations. We first give a precise definition for SLP polynomials.
Definition 4.5
An SLP over a ring R is a branchless sequence of arithmetic instructions thatrepresents a polynomial function. It takes as input a vector ( a , . . . , a n ) and outputs a vector ( b , . . . , b L ) by way of a series of instructions Γ i : 1 ≤ i ≤ L of the form Γ i : b i ← α (cid:63) i β , where (cid:63) i is an operation (cid:48) + (cid:48) , (cid:48) − (cid:48) or (cid:48) × (cid:48) , and α, β ∈ R (cid:83) { a , . . . , a n } (cid:83) { b , . . . , b i − } . The inputs and outputsmay belong to R or a ring extension of R . We say that an SLP computes a multivariate polynomial f ∈ R [ x , . . . , x n ] if it sets b L to be f ( a , . . . , a n ) . We now give the interpolation algorithm for univariate polynomials, where f ∗ in the input isintroduced because the multivariate interpolation algorithm in Section 5 will use it. Algorithm 4.6 (UIPoly)Input:
An SLP S f that computes f ( x ) ∈ R [ x ], f ∗ ∈ R [ x ], T ≥ max { f, f ∗ } , T ≥ f − f ∗ ),( T ≥ T ), D > max { deg f, deg f ∗ } . Output:
The exact form of f − f ∗ . Step 1:
Let N = max { , (cid:100) ( T −
1) log D (cid:101)} , N = (cid:100) T log D (cid:101) , N = max { N , N + N − } . Step 2:
Find the first N primes p , . . . , p N . Step 3:
Compute the smallest K D such that p · · · p K D ≥ D . Step 4:
For j = 1 , , . . . , N , probe f mod ( D,p j ) . Let f j = f mod ( D,p j ) − f ∗ mod ( D,p j ) and h ∗ = 0. Step 5:
Loop
Let α = max { f j | j = 1 , , . . . , N } and j the smallest number such that f j = α . If α = 0, then return h ∗ . For k = 1 , , . . . , K D , probe f mod ( D,p j p k ) and let g k = f mod ( D,p j p k ) − f ∗ mod ( D,p j p k ) . Let U f − f ∗ − h ∗ D,p j := UTerms ( f j , g , g , . . . , g K D , p j , D ). Let h = 0. For each u ∈ U f − f ∗ − h ∗ D,p j , if { j | f j − u mod ( D,p j ) ) < f j ) , j = 1 , . . . , N + N − } ≥ N then h := h + u . 9 .6: Let h ∗ = h ∗ + h , T = T − h , N = max { , (cid:100) ( T −
1) log D (cid:101)} , N = (cid:100) T log D (cid:101) , N =max { N , N + N − } . For j = 1 , , . . . , N , let f j = f j − h mod ( D,p j ) . Theorem 4.7
Algorithm 4.6 returns f − f ∗ using O ∼ ( LT log D + LT log D ) ring operations in R and similarly many bit operations, where L is the size of the SLP representation for f . Specially,when f ∗ = 0 , the algorithm returns f .Proof. We first prove the correctness of the theorem. We claim that each loop of Step 5 will obtainat least half of the terms of f − f ∗ − h ∗ . Then, the algorithm will return the correct f − f ∗ byrunning at most log T times of the loop in Step 5. In Step 5.1, Theorem 3.4 is used to find an okayprime p j . In Step 5.4, by Lemma 4.1 and Theorem 3.4, at least half of the terms of f − f ∗ − h ∗ arein U f − f ∗ − h ∗ D,p j . In Step 5.5, Theorem 2.4 is used to select the elements of M f − f ∗ − h ∗ from U f − f ∗ − h ∗ D,p j .In summary, at Step 5.6, h contains at least half of the terms of f − f ∗ − h ∗ and the claim is proved.Then, the correctness of the algorithm is proved.We now analyse the complexity of the algorithm, which comes from Step 4 and Step 5. Thecomplexity of other steps are lower than these two steps.In Step 3, since the bit complexity of finding the first N primes is O ( N log N log log N ) by [20,p.500,Thm.18.10] and N is O ∼ ( T log D ), the bit complexity of Step 3 is O ∼ ( T log D ).In Step 4, we probe N univariate polynomials f mod ( D,p j ) and probing f mod ( D,p j ) costs O ∼ ( Lp j ), since S f is of length L and the univariate polynomials in the procedure have degrees < p j . Since p i is of O ∼ ( T log D ) and N is O ∼ ( T log D ), the cost of probing f mod ( D,p j ) is O ∼ ( LT log D ) ring andbit operations. It needs O ∼ ( T log D ) ring operations and O ∼ ( T log D ) bit operations to obtain f ∗ mod ( D,p j ) . Then the total complexity of Step 4 is O ∼ ( LT log D ) ring and bit operations.We now consider Step 5. We first consider the complexity of each loop of the this step. InStep 5.1, since N is of O ( T log D ) and the terms of ( f − f ∗ − h ∗ ) mod ( D,p j ) is no more than T , it needs O ∼ ( T log D ) bit operations.In Step 5.3, since K D is of O (log D ) and p j p k is of O ∼ ( T log D ), we need O ∼ ( LT log D )arithmetic operations in R and similarly many bit operations to obtain f mod ( D,p j p k ) . We need O ( K D f ∗ ) ring operations and O ∼ ( K D f ∗ max { log deg f ∗ , log( p j p k ) } ) bit operations to obtain f ∗ mod ( D,p j p k ) . Since f ∗ ≤ T , deg f ∗ ≤ D , p k , K D is of O (log D ), and p j is of O ( T log D ), the costis O ( T log D ) ring operations and O ∼ ( T log D max { log D, log T + log log D } ) = O ∼ ( T log D ) bitoperations. Then, the total complexity of this step is O ∼ ( LT log D ) arithmetic operations in R and similarly many bit operations.In Step 5.4, by Lemma 4.4, the complexity is O ( T log D ) arithmetic operations in R and O ∼ ( T log D ) bit operations.In Step 5 .
5, in order to determine whether f j − u mod ( D,p j ) ) < f j ), we just need to determinewhether u mod ( D,p j ) is a term of f j . We sort the terms of f j such that they are in ascending orderaccording to their degrees, which costs O ∼ (( N + N ) T ) = O ∼ ( T log D ) bit operations, since N + N is O ∼ ( T log D ). To find whether f j has a term with degree deg( u mod ( D,p j ) ), we need O (log T )comparisons. Since the height of the degree is O (log D ), it needs O (log T log D ) bit operations.To compare the coefficient, it needs one arithmetic operation. So it totally needs O (log T log D )bit operations and O (1) arithmetic operation to compare f j − u mod ( D,p j ) ) with f j ). Hence, thetotal complexity of Step 5.5 is O ∼ ( U f − f ∗ D,p j ( N + N ) log T log D + ( N + N ) T ) bit operations10nd O ( U f − f ∗ D,p j ( N + N )) arithmetic operations in R . Since U f − f ∗ D,p j ≤ T and N + N − O ( T log D ), the total complexity is O ∼ ( T log D ) arithmetic operations and O ∼ ( T log D ) bitoperations.So the total complexity of each loop of Step 5 is O ∼ ( LT log D + T log D ) arithmetic operationsand O ∼ ( LT log D + T log D ) bit operations, which comes from Step 5.3 and Step 5.5, respectively.Since each loop of Step 5 will obtain at least half of f − f ∗ − h ∗ , Step 5 has at most O (log T ) loops.So the total complexity of Step 5 is O ∼ ( LT log D + T log D ) ring operations and O ∼ ( LT log D + T log D ) bit operations.Combing with the complexity of Step 4, the total complexity of the algorithm is O ∼ ( LT log D + LT log D + T log D ) = O ∼ ( LT log D + LT log D ) ring and bit operations, which comes fromStep 4 and Step 5.3. The theorem is proved.For an n -variate polynomial f of degree < D , we can use the Kronecker substitution [18]to reduce the interpolation of f to that of a univariate polynomial of degree D n , which can becomputed with Algorithm UIPoly ( S f , , T, T, D ). By Theorem 4.7, we have Corollary 4.8
For an SLP f ∈ R [ X ] with T ≥ f and D < deg( f ) , we can find f using O ∼ ( Ln T log D + Ln T log D ) ring operations in R and a similar number of bit operations. In the next section, we will give a multivariate polynomial interpolation algorithm which hasbetter complexity.
In this section, we will give a new multivariate interpolation algorithm which is quadratic in n ,while the algorithm given in Corollary 4.8 is cubic in n . The algorithm is quite similar to Algorithm4.6 and works as follows. First, we use Theorem 3.4 to find an “ok” prime p for f . Second, weuse a modified Kronecker substitution to obtain a set S of terms, which contains at leats half ofthe terms of f . Finally, we use Theorem 2.4 to identify the terms of f from S . The multivariateinterpolation algorithm will call Algorithm 4.6. x p − Let f be a multivariate polynomial, M f the set of terms in f , t = f , D > deg( f ), T ≥ f , and p ∈ N > . Consider the modified Kronecker substitutions: f ( D,p ) = f ( x, x mod ( D,p ) , . . . , x mod ( D n − ,p ) ) (4) f ( D,p,k ) = f ( x, x mod ( D,p ) , . . . , x p + mod ( D k − ,p ) , . . . , x mod ( D n − ,p ) ) (5)where k ∈ { , , . . . , n } , f ( D,p ) comes from the substitutions x i = x mod ( D i − ,p ) , i = 1 , , . . . , n , and f ( D,p,k ) comes from the substitutions x i = x mod ( D i − ,p ) , i = 1 , , . . . , n, i (cid:54) = k, x k = x p + mod ( D k − ,p ) .Note that when n = 1, f ( D,p ) = f ( D,p,k ) = f ( x ). Substitution (4) was introduced in [17] andsubstitution (5) is introduced in this paper. We havedeg f ( D,p ) ≤ Dp and deg f ( D,p,k ) ≤ Dp. (6)Similar to Definition 2.1, a term cm is said to be a collision in f ( D,p ) or in f ( D,p,k ) , if there existsan aw ∈ M f \{ cm } such that m ( D,p ) = w ( D,p ) or m ( D,p,k ) = w ( D,p,k ) .11e now show how to compute S f ( D,p ) and S f ( D,p,k ) . Lemma 5.1
Let S f be an SLP procedure to compute f , which has length L . Then we can designa procedure S f ( D,p ) ( S f ( D,p,k ) ) for f ( D,p ) ( f ( D,p,k ) ), which has length L and costs extra O (log D + n log p + L log p ) bit operations. Probing f ( D,p ) mod ( x q − from S f ( D,p ) costs O ∼ ( Lq + L log p ) arithmetic operations and similarly many bit operations.Proof. Define a procedure S f ( D,p ) for f ( D,p ) as follows. Suppose we want to compute f ( D,p ) ( a ) forsome a in R or an extension of R . Assume S f consists of the operations Γ i : b i ← α i (cid:63) i β i , i =1 , , . . . , L with input { a , a , . . . , a n } . Now we define the i -th instruction Γ i in S f ( D,p ) .Γ i := b i ← α i (cid:63) i β i if both α i and β i are not in { a , a , . . . , a n } b i ← a mod ( D j − ,p ) (cid:63) i β i if α i is a j , but β i is not in { a , a , . . . , a n } b i ← α i (cid:63) i a mod ( D j − ,p ) if α i is not in { a , a , . . . , a n } , but β i is a j b i ← a mod ( D j − ,p ) (cid:63) i a mod ( D j − ,p ) if α i is a j and β i is a j (7)Now we analyse the complexity of the procedure. In order to obtain mod ( D j − , p ) , j =1 , , . . . , n , it needs O (log D + n log p ) bit operations. To obtain all Γ i , it needs O ( L ) arithmeticoperations. Since the height of the data is log p , it needs O ( L log p ) bit operations. So it totallyneeds O (log D + n log p + L log p ) bit operations.The univariate polynomial f ( D,p ) mod ( x q −
1) can be computed from S f ( D,p ) as follow: firstwe replace a i by x . During the computing, we always use the mod ( x q −
1) to reduce the degree.So the degree of x is less than q . If the length of S f is L , then probing f ( D,p ) mod ( x q −
1) from S f ( D,p ) costs O ∼ ( LM ( q ) + L log p ) arithmetic operations in R plus similar bit operations, where LM ( q ) is the complexity of multiplying two univariate polynomials with degrees < q . By [10], wemay assume M ( q ) is O ( q log q log log q ). So it costs O ∼ ( Lq + L log p ) ring operations and similarlymany bit operations. The definition of S f ( D,p,k ) is the same as S f ( D,p ) . The only difference is thatwhen j = k , then replace a k by a mod ( D k − ,p )+ p . Remark 5.2
From the above Lemma, although S f ( D,p ) , S f ( D,p,k ) are not SLP procedures, we stillcan probe f ( D,p ) mod ( x q − from them. Since in the following algorithms, p is O ∼ ( nT log D ) and q is O ∼ ( T log( nD )) , the complexity of the probing is O ∼ ( Lq ) . So we can still regard S f ( D,p ) , S f ( D,p,k ) as SLP procedures of length L . Let f mod ( D,p ) = a x d + a x d + · · · + a r x d r ( d < d < · · · < d r ) (8)Since f ( D,p ) mod ( x p −
1) = f ( D,p,k ) mod ( x p −
1) = f mod ( D,p ) , for k = 1 , , . . . , n , we can write f ( D,p ) = f + f + · · · + f r + g (9) f ( D,p,k ) = f k, + f k, + · · · + f k,r + g k where f i mod ( x p −
1) = f k,i mod ( x p −
1) = a i x d i , g mod ( x p −
1) = g k mod ( x p −
1) = 0.Similar to (3), we define the following key notation M fD,p = { a i x e i, · · · x e i,n n | a i is from (8) for some i ∈ [1 , r ] , and121 : f i = a i x u i , f k,i = a i x b k,i , k = 1 , , . . . , n. (10)M2 : e i,k = b k,i − u i p ∈ N , k = 1 , , . . . , n. M3 : u i = e i, + e i, mod ( D, p ) + · · · + e i,n mod ( D n − , p ) . M4 : n (cid:88) j =1 e i,j < D. } Lemma 5.3
Let f = (cid:80) ti =1 c i m i ∈ R [ X ] , D > deg f . If c i m i is not a collision in f mod ( D,p ) , then c i m i ∈ M fD,p .Proof. It suffices to show that c i m i satisfies the conditions of the definition of M fD,p . Assume m i = x e x e · · · x e n n . Since c i m i is not a collision in f mod ( D,p ) , without loss of generality, assume ( c i m i ) mod ( D,p ) = a x d and d = mod ( (cid:80) nj =1 e j D j − , p ), where a x d is defined in (8). It is easy to see that c i m i is also not a collision in f ( D,p ) and in f ( D,p,k ) . Hence, f = a x u for u = (cid:80) ni =1 e i mod ( D i − , p ); b k, = u + pe k . Clearly, M1, M2 and M3 are correct. Since deg( m i ) = (cid:80) nj =1 e i,j < D , M4 iscorrect.Now we give the following algorithm to compute M fD,p , whose correctness is obvious. Algorithm 5.4 (MTerms)Input:
Univariate polynomials f mod ( D,p ) , f ( D,p ) , f ( D,p,k ) , where k = 1 , , . . . , n , a prime p , D > deg( f ). Output: M fD,p . Step 1:
Write f mod ( D,p ) , f ( D,p ) , and f ( D,p,k ) as forms (8) and (9) f mod ( D,p ) = a x d + a x d + · · · + a r x d r f ( D,p ) = f + f + · · · + f r + gf ( D,p,k ) = f k, + f k, + · · · + f k,r + g k , ( k = 1 , , . . . , n ) . Step 2:
Let S = {} . For i = 1 , , . . . , r do a: If one of f i , f ,i , . . . , f n,i is not of the following form: f i = a i x u i , f k,i = a i x b k,i , then break. b: Let e i,k = b k,i − u i p for k = 1 , , . . . , n . If e i,k / ∈ N , then break. c: If u i (cid:54) = e i, + e i, mod ( D, p ) + · · · + e i,n mod ( D n − , p ), then break; d: If (cid:80) nj =1 e i,j ≥ D , then break; e: Let S = S (cid:83) { a i x e i, · · · x e i,n n } . Step 3
Return S . Lemma 5.5
Algorithm 5.4 needs O ( nT ) arithmetic operations in R and O ∼ ( nT log( pD )) bit op-erations.Proof. In Step 1, we need to write f ( D,p ) and f ( D,p,k ) as the desired form. This can be done in threesteps. First, we perform the modular operation mod p on every degree of f ( D,p ) , f ( D,p,k ) , whichcosts O ( nT log( pD )) bit operations, since each of f ( D,p ) , f ( D,p,k ) has no more than T terms and the13eight of the degree is O (log( pD )). Second, we sort the terms of f mod ( D,p ) , f ( D,p ) , f ( D,p,k ) into ascendingorder according to the new degree module p , which costs O ∼ ( nT log( p )) bit operations, since thedegrees are < p . In order to check whether f i mod ( x p −
1) = f k,i mod ( x p −
1) = a i x d i , we need O ( T n ) operations over R . Finally, f i , f k,i can be obtained with O ( T n ) comparisons of the degrees,which costs O ∼ ( nT log( p )) bit operations and O ( T n ) R -operations. So, the total complexity ofStep 1 is O ∼ ( nT log( pD )) bit operations.For Step 2, we first consider the complexity of one loop. Since the height of the degrees of f mod ( D,p ) , f ( D,p ) , f ( D,p,k ) are O (log( pD )), Steps a , b , c , and d costs O ( n log( pD )) bit operations. Sincewe have at most T loops, the total complexity is is O ( nT log( pD )) bit operations.In a of Step 2, since f ( D,p ) , f ( D,p,k ) ≤ T , it totally needs O ∼ ( nT ) bit operations to deter-mine whether f i , f k,i ) = 1. To compare the coefficients of f i , f k,i , it needs O ( nT ) arithmeticoperations in R . We prove the lemma In this section, we give the interpolation algorithm for multivariate polynomial. We first give asub-algorithm, which computes f ( D,p ) , and f ( D,p,k ) efficiently. Algorithm 5.6 (Substitution)Input:
A polynomial f ∈ R [ X ], a prime p , a number D ∈ N with deg f < D . Output:
The univariate polynomials f ( D,p ) , and f ( D,p,k ) , k = 1 , , . . . , n . Step 1:
Assume f = c m + c m + · · · + c t m t , where m i = x e i, x e i, · · · x e i,n n , i = 1 , , . . . , t . Step 2:
Let u = 1;For i = 2 , , . . . , n do u i = mod ( u i − D, p ). Step 3:
Let h = 0. For i = 1 , , . . . , n , let h i = 0; Step 4:
For i = 1 , , . . . , t do a: Let d = 0. b: For k = 1 , , . . . , n , let d = d + e i,k u k . c: h = h + c i x d ; d: For k = 1 , , . . . , n , let h k := h k + c i x d + e i,k p . Step 5:
Return h , h i , i = 1 , , . . . , n ; Lemma 5.7
Algorithm 5.6 is correct. The complexity is O ∼ ( nt log p + nt log( D )) bit operationsand O ( nt ) arithmetic operations in R .Proof. In Step 2, u i = mod ( D i − , p ). In b of Step 4, d is the degree of m i ( D,p ) , so h is f ( D,p ) afterfinishing Step 4. In d , since deg( m i ( D,p,k ) ) = deg( m i ( D,p ) ) + pe i,k , h k is f ( D,p,k ) after finishing Step4. So the correctness is proved.Now we analyse the complexity. In Step 2, it needs O ( n max { log D, log p } ) bit operations. In b of Step 4, it needs O ( nt ) arithmetic operations in Z . Since deg( m i ( D,p ) ) is O ( p · deg f ) ≤ O ( pD ),14he bit operations is O ( nt log( pD )). In c and d , it needs O ∼ ( nt log( pD )) bit operations and O ( nt )arithmetic operations in R .We now give the interpolation algorithm. Algorithm 5.8 (MPolySI)Input:
An SLP S f that computes f ∈ R [ X ], T ≥ f , D > deg f . Output:
The exact form of f . Step 1:
Let N = max { , (cid:100) n ( T −
1) log D (cid:101)} , N = (cid:100) nT log D (cid:101) , N = max { N , N + N − } . Step 2:
Find the first N different primes p , p , . . . , p N . Step 3:
For j = 1 , , . . . , N , probes f mod ( D,p j ) . Let f mod j = f mod ( D,p j ) . Step 4:
Let h = 0 and T = T . Step 5:
Loop
Let α = max { f mod j | j = 1 , , . . . , N } and j the smallest number such that f mod j = α . If α = 0 return h . { f ∗ , f ∗ , . . . , f ∗ n } = Substitution ( h, p j , D ). Let f j = UIPoly ( S f ( D,pj , f ∗ , T, T , Dp j ). For k = 1 , , . . . , n , let g k = UIPoly ( S f ( D,pj ,k ) , f ∗ k , T, T , Dp j ). Let M f − hD,p j := MTerms ( f mod j , f j , g , g , . . . , g n , p j , D ). Let s = 0. For each u ∈ M f − hD,p j , if { j | f mod j − u mod ( D,p j ) ) < f mod j ) , j = 1 , . . . , N + N − } ≥ N , then s := s + u . Let h = h + s , T = T − s , N = max { , (cid:100) n ( T −
1) log D (cid:101)} , N = (cid:100) nT log D (cid:101) , N =max { N , N + N − } . For j = 1 , , . . . , N , let f mod j = f mod j − s mod ( D,p j ) . Theorem 5.9
Algorithm 5.8 finds f using O ∼ ( Ln T log D + LnT log D ) ring operations in R and similar bit operations.Proof. The algorithm is quite similar to the univariate interpolation Algorithm 4.6. So, we will onlygive the sketch of the proof and give detailed proof only for those steps and are essentially differentfrom that in Algorithm 4.6. By Theorem 3.4 and Lemma 5.3, at least half of terms of f − h arein M f − hD,p j obtained in Step 5.6. In Step 5.7, Theorem 2.4 is used to select the elements of M f − h from M f − hD,p j . So at least half of the terms of f − h will be found in each loop of Step 5. Then thecorrectness of the algorithm is proved.We now analyse the complexity of the algorithm, which comes from that of Steps 3 and 5.In Step 2, since the bit complexity of finding the first N primes is O ( N log N log log N ) by [20,p.500,Thm.18.10] and N is O ∼ ( nT log D ), the bit complexity of Step 2 is O ∼ ( nT log D ).15n Step 3, we probe S f for O ( nT log D ) times. Since p i is O ∼ ( nT log D ), the cost of probes is O ∼ ( Ln T log D ) ring and bit operations.We now consider the complexity of one loop for Step 5. In Step 5.1, since N is of O ( nT log D )and f mod j ≤ T , the bit complexity is O ∼ ( nT log D ).In Steps 5.3, by Lemma 5.7, since p i is O ∼ ( nT log D ), the complexity is O ∼ ( n log D + nT log p j + nT log D ) = O ∼ ( nT log D ) bit operations and O ( nT ) arithmetic operations in R .In Steps 5.4 and 5.5, by Remark 5.1, we can regard S f ( D,pj and S f ( D,pj ,k ) as SLP proceduresof length L . Since the numbers of terms and degrees of ( f − h ) ( D,p j ) , ( f − h ) ( D,p j ,k ) are boundedby O ( T ) and O ∼ ( nT D ), by Theorem 4.7, the complexity is O ∼ ( LnT log D + LnT log D ) ringand bit operations.In Step 5.6, by Theorem 5.5, the complexity is O ∼ ( nT log D ) bit operations and O ( nT ) ringoperations.In Step 5.7, to all the u mod ( D,p j ) , we need O ∼ ( n ( N + N ) M f − hD,p j log( Dp ) bit operations. The prooffor rest of this step is similar to that of Step 5.5 of Algorithm 4.6. The complexity is O ∼ ( n ( N + N ) M f − hD,p j log( Dp j )+ M f − hD,p j log T ( N + N ) log( Dp j )) bit operations and O ( M f − hD,p j ( N + N ))ring operations. Since M f − hD,p j ≤ T , p j = O ∼ ( nT log D ), and N + N = O ( nT log D ), it needs O ∼ ( n T log D ) bit operations and O ∼ ( nT log D ) ring operations.In Step 5.9, we need n s operations in Z to obtain s mod ( D,p j ) . Subtract s mod ( D,p j ) from f j needs s log T operations in Z and s arithmetic operation in R . Since the height of the data is O (log( nT D )) and we need update N polynomials, the complexity is O (( n s log( nT D ) + s log T log( nT D )) N ) bit operations and O ( sN ) ring operations. Since the sum of s is t , it total costs O ∼ ( nT log D ) ring operations and O ∼ ( n T log D ) bit operations.Then the total complexity of one loop of Step 5 is O ∼ ( LnT log D + LnT log D ) ring operationsand O ∼ ( LnT log D + LnT log D + n T log D ) bit operations, which come from Steps 5.5, 5.7,and 5.9. Since every loop of Step 5 finds at least half of the terms in f − h , the loop runs at most O (log T ) times. So, the total complexity of Step 5 is O ∼ ( LnT log D + LnT log D ) ring operationsand O ∼ ( LnT log D + LnT log D + n T log D ) bit operations. Plus the complexity of Step 3,the complexity of the algorithm is O ∼ ( Ln T log D + LnT log D ) ring and bit operations, whichare from Step 3 and Step 5.5. In this section, we give Monte Carlo interpolation algorithms for multivariate polynomials, whichcould be considered as probabilistic versions of Algorithm 5.8.The following theorem shows how to use a probabilistic method to obtain a p such that thenumber of collision terms of f in f mod ( D,p ) is very small, which is a probabilistic version of Theorem3.4. Theorem 6.1
Let f = (cid:80) ti =1 c i m i ∈ R [ X ] , T ≥ f, D > deg f , N = max { , (cid:98) n ( T −
1) log D (cid:99)} ,and p , p , . . . , p N be N different primes, < ε < , k = (cid:100) log ( ε ) (cid:101) . If j , j , . . . , j k are randomlychosen from [1 , N ] and j ∈ { j , . . . , j k } is the integer such that f mod ( D,p j ) = max { f mod ( D,p j ) , f mod ( D,p j ) ,. . . , f mod ( D,p jk ) } . Then with probability ≥ − ε , C f ( D,p j ) ≤ (cid:98) T (cid:99) .Proof. If T = 1, then N = 1, the proof is obvious. So we assume T ≥
2, then N = (cid:98) n ( T −
16) log D (cid:99) . First we claim that if randomly choose an integer j in [1 , N ], then with probabilityat least , C f ( D,p j ) < T . It suffices to show that there exist at least N integers j in [1 , N ] suchthat C f ( D,p j ) < T . We prove this by contradiction. Assume α , α , . . . , α N +1 are N + 1 differentintegers in [1 , N ] such that C f ( D,p αi ) ≥ T . By Lemma 3.3, p (cid:100)C f ( D,pαi ) / (cid:101) α i divides A , where A isdefined in Lemma 3.3. Since p α i are different and C f ( D,p αi ) ≥ T , ( p α p α · · · p α N +1 ) (cid:100) T (cid:101) divides A .Now we have ( p α p α · · · p α N +1 ) (cid:100) T (cid:101) ≥ ( N +1) (cid:100) T (cid:101) ≥ (cid:100) n ( T −
1) log D (cid:101)(cid:100) T (cid:101) ≥ nT ( T −
1) log D = D nT ( T − , which contradicts to A ≤ ( D n − t ( t − . So in { p , p , · · · , p N } , there are at least N primes p i such that C f ( D,p i ) < T . We have proved the claim.If there exists at least one j i in { j , j , . . . , j k } such that C f ( D,p ji ) < T , then by Corollary3.2, C f ( D,p j ) ≤ C f ( D,p ji ) < T . So C f ( D,p j ) ≥ T only when all C f ( D,p ji ) ≥ T, i = 1 , , . . . , k , bythe claim just proved, the probability for this to happen is at most ( ) k ≤ ( ) log ε = ε . Since C f ( D,p j ) < T implies that C f ( D,p j ) ≤ (cid:98) T (cid:99) , the probability of C f ( D,p j ) ≤ (cid:98) T (cid:99) is at least 1 − ε .The theorem is proved. Remark 6.2
Note that the result C f ( D,p j ) ≤ (cid:98) T (cid:99) of Theorem 6.1 is different with that of Theorem3.4. To find a p such that at least (cid:100) t (cid:101) of the terms of f are not collisions in f mod ( D,p j ) is not enoughfor our probabilistic algorithm. Lemma 6.3 If ε ∈ (0 , and a ≥ , then (1 − ε ) a ≥ − aε .Proof. By Taylor expansion, we have (1 + x ) a = 1 + ax + a ( a − (1+ θx ) a − , where θ ∈ (0 , x = − ε , then (1 − ε ) a = 1 − aε + a ( a − (1 − θε ) a − . Since θ ∈ (0 , ε ∈ (0 ,
1) and a ≥
1, wehave a ( a − (1 − θε ) a − ≥
0. So we have (1 − ε ) a ≥ − aε .We first consider interpolation over an arbitrary computable ring. For the univariate interpola-tion algorithm, we use the following algorithm given in [4], which is the fastest known probabilisticalgorithm over arbitrary rings. Theorem 6.4 [4] Let f ∈ R [ x ] , where R is any ring. Given any SLP of length L that computes f ,and bounds T and D for the sparsity and degree of f , one can find all coefficients and exponents of f using O ∼ ( LT log D + LT log D log ν ) ring operations in R , plus a similar number of bit operations.The algorithm is probabilistic of the Monte Carlo type: it can generate random bits at unit cost andon any invocation returns the correct answer with probability greater than − ν , for a user-suppliedtolerance < ν < . We use
PUniPoly ( S f , f ∗ , T , D, ν ) to denote the algorithm in Theorem 6.4, where f ∗ is acurrent approximation to f and f − f ∗ ) ≤ T , f ∗ ≤ T, max { deg f ∗ , deg f } < D .Now we give an algorithm which interpolates at least half of the terms. Algorithm 6.5 (HalfPoly)Input:
An SLP S f that computes f , f ∗ ∈ R [ X ], T ≥ max { f, f ∗ } , T ≥ f − f ∗ ), D > max { deg f, deg f ∗ } , a tolerance ν such that 0 < ν < Output:
With probability ≥ − ν , return a polynomial f ∗∗ such that f − f ∗ − f ∗∗ ) ≤ (cid:98) T (cid:99) .17 tep 1: Let N = max { , (cid:100) n ( T −
1) log D (cid:101)} , ε = νn +1 , and k = (cid:100) log ε (cid:101) . Find the first 2 N primes p , p , . . . , p N . Step 2:
Let j , j , . . . , j k be randomly chosen from [1 , N ]. Delete the repeated numbers, we stilldenote these integers as j , j , . . . , j k . Step 3:
For i = 1 , , . . . , k , probe f mod ( D,p ji ) . Let g i = f mod ( D,p ji ) − f ∗ mod ( D,p ji ) . Step 4:
Let α = max { g i | i = 1 , , . . . , k } and j satisfying g j = α . If α ≥ T , then returnfailure. Step 5:
Let { q, q , . . . , q n } = Substitution ( f ∗ , p j , D ). Step 6:
Let η = PUniPoly ( S f ( D,pj , q, T , p j D, ε ), η i = PUniPoly ( S f ( D,pj ,i ) , q i , T , p j D, ε ) ,i = 1 , . . . , n . If η or η i is failure, then return failure. Step 7:
Let M = MTerms ( g j , η, η , η , . . . , η n , D, p j ). Step 8:
Return f ∗∗ = (cid:80) s ∈ M s . Lemma 6.6
Algorithm 6.5 computes f ∗∗ such that f − f ∗ − f ∗∗ ) ≤ (cid:98) T (cid:99) with probability at least − ν . The algorithm costs O ∼ ( LnT log D + LnT log D log ν ) ring operations in R and a similarmany bit operations.Proof. We first show that Algorithm 6.5 returns the polynomial f ∗∗ such that f − f ∗ − f ∗∗ ) ≤ (cid:98) T (cid:99) with probability 1 − ν . In Step 4, by Theorem 6.1, with probability 1 − ε , C f − f ∗ ( D,p j ) ≤ (cid:98) T (cid:99) . If j satisfies C f − f ∗ ( D,p j ) ≤ (cid:98) T (cid:99) and η = ( f − f ∗ ) ( D,p j ) , η i = ( f − f ∗ ) ( D,p j ,i ) , i = 1 , , . . . , n , thenby Lemma 5.3, f − f ∗ contains at most (cid:98) T (cid:99) terms which are not in f ∗∗ . Since the terms of f ∗∗ which are not in M f − f ∗ come from at least two terms in f − f ∗ , there exist at most (cid:98) T (cid:99) terms of f ∗∗ not in f − f ∗ . So f − f ∗ − f ∗∗ ) ≤ (cid:98) T (cid:99) + (cid:98) T (cid:99) ≤ T < T and we have f − f ∗ − f ∗∗ ) ≤ (cid:98) T (cid:99) . In Step6, by Theorem 6.4, the probability of obtaining the correct( f − f ∗ ) ( D,p j ) (( f − f ∗ ) ( D,p j ,i ) ), or η = ( f − f ∗ ) ( D,p j ) ( η i = ( f − f ∗ ) ( D,p j ,i ) ), is ≥ − ε , so theprobability of obtaining correct polynomials ( f − f ∗ ) ( D,p j ) and (( f − f ∗ ) ( D,p j ,i ) ) , i = 1 , . . . , n is ≥ (1 − ε ) n +2 . By Lemma 6.3, (1 − ε ) n +1 ≥ − ( n + 1) ε . Since ε = νn +1 , (1 − ε ) n +1 ≥ − ν . Hence,with probability ≥ − ν , we obtain an f ∗∗ satisfying f − f ∗ − f ∗∗ ) ≤ (cid:98) T (cid:99) . The correctness ofthe lemma is proved.Now we analyse the complexity. Since the bit complexity of finding the first 2 N primes is O ( N log N log log N ) by [20, p.500,Thm.18.10] and N is O ∼ ( nT log D ), the bit complexity of Step1 is O ∼ ( nT log D ).In Step 2, since probabilistic machines flip coins to decide binary digits, each of these randomchoices can be simulated with a machine with complexity O (log 2 N ). So the complexity of Step 2is O (log ε log( nT log D )) bit operations. Since O (log ε ) is O (log n + log ν ), the bit complexity ofStep 2 is O (log n log( nT log D ) + log( nT log D ) log ν ).In Step 3, we probe k = (cid:100) log ε (cid:101) times for f . Since p j i is of O ∼ ( nT log D ), the cost of the probesis O ∼ ( LnT log D log ε ) ring operations and a similar many bit operations. So the complexity ofstep 3 is O ∼ ( LnT log D log ν ) arithmetic operations in R and a similar number of bit operations.In Step 4, we find the integer j . Since f − f ∗ ) mod ( D,p ji ) ≤ T , it needs at most O ∼ ( T log ε ) bitoperations to compute all f − f ∗ ) mod ( D,p ji ) , i = 1 , , . . . , k . Find j needs O ∼ ( T ) bit operations. Sothe bit complexity of Step 4 is O ∼ ( T log n + T log ν ).18n Step 5, by Lemma 5.7, it needs O ∼ ( nT log D ) bit operations and O ( nT ) ring operations.In Step 6, we call n + 1 times Algorithm PUniPoly . Since the terms and degrees of ( f − f ∗ ) ( D,p j ) , ( f − f ∗ ) ( D,p ji ,k ) are respectively bounded by T and 2 p j D , by Theorem 6.4, the complexityof Step 6 is O ∼ ( LnT log ( p j D ) + LnT log( p j D ) log ε )) arithmetic operations in R and plus asimilar number of bit operations. Since ε = νn +1 , p j is O ∼ ( nT log D ), and 2 p j D is O ∼ ( nT D ),the complexity of step 6 is O ∼ ( LnT log D + LnT log D log ν )) ring operations and similar bitoperations.In Step 7, by Theorem 5.5, the complexity is O ( nT ) arithmetic operation in ring R and O ∼ ( nT log pD ) bit operations.It is easy to see that the complexity is dominated by Step 3 and Step 6. The theorem is proved.We now give the complete interpolation algorithm. Algorithm 6.7 (MCMulPoly)Input:
An SLP S f that computes f ∈ R [ X ], T > f , D > deg f , and µ ∈ (0 , Output:
With probability at least 1 − µ , return f . Step 1:
Let h = 0 , T = T, ν = µ (cid:100) log T (cid:101) +1 . Step 2:
While T > b: Let g = HalfPoly ( S f , h, T, T , D, ν ). If g = f ailure , then return failure. c: Let h = h + g , T = (cid:98) T (cid:99) . Step 3:
Return h . Theorem 6.8
Algorithm 6.7 computes f with probability ≥ − µ . The algorithm costs O ∼ ( LnT log D + LnT log D log µ ) ring operations in R and a similar number of bit operations.Proof. In Step 2, by Lemma 6.6, f − h − g ) ≤ (cid:98) T (cid:99) with probability ≥ − ν . By Lemma 6.3, Step 2will run at most k = (cid:100) log T (cid:101) + 1 times and return the correct f with probability ≥ (1 − ν ) k ≥ − µ .The correctness of the theorem is proved.Now we analyse the complexity. In Step 2, we call at most O (log T ) times Algorithm HalfPoly .Since the terms and degrees of f − h, f are respectively bounded by T and D , by Theorem 6.6,the complexity of Step 2 is O ∼ ( LnT log D + LnT log D log ν )) ring and bit operations. Since ν = µ (cid:100) log T (cid:101) +1 , the complexity of Step 2 is O ∼ ( LnT log D + LnT log D log µ )) ring and bit operations.The theorem is proved. Remark 6.9
Set µ = 1 / . Then Algorithm 6.7 computes f with probability at lest . The algorithmcosts O ∼ ( LnT log D ) ring operations in R and a similar number of bit operations. We now consider interpolation over finite fields. In [5], Arnold, Giesbrecht & Roche gave a newunivariate interpolation algorithm for the finite field with better complexities.
Theorem 6.10 [5] Let f ∈ F q [ x ] with at most T non-zero terms and degree at most D , and let <η ≤ . Suppose we are given an SLP S f of length L that computes f . Then there exists an algorithmthat interpolates f , with probability at least − η , with a cost of O ∼ ( LT log D (log D + log q ) log η ) bit operations.
19e use
PUniPoly ( f, f ∗ , T , D, η ) to denote the algorithm in Theorem 6.10, where f ∗ is acurrent approximation to f and f − f ∗ ) ≤ T , f ∗ ≤ T, deg f ∗ < D, deg f < D .Replacing Algorithm PUniPoly with Algorithm PUniPoly in Step 6 of Algorithm HalfPoly ,we obtain a multivariate interpolation algorithm for finite fields. We assume S f can evaluate in anextension field of F q and have the following result. Theorem 6.11
Let f ∈ F q [ X ] be given as an SLP, with at most T non-zero terms and degree atmost D , and let ≤ µ < / . Then we can interpolate f , with probability at least − µ , with acost of O ∼ ( LnT log D (log D + log q ) log µ ) bit operations.Proof. The proof is the same as that of Theorem 6.8. The only difference is to use Theorem 6.11instead of Theorem 6.4.
In this section, practical performances of the new algorithms implemented in Maple will be pre-sented. The data are collected on a desktop with Windows system, 3.60GHz Core i7-4790 CPU,and 8GB RAM memory. The Maple codes can be found in
We randomly construct five polynomials, then regard them as SLP polynomials and reconstructthem with the Algorithms 4.6 , 5.8 and 6.7. We do not collect the time of probes. We just test thetime of recovering f from the univariate polynomial f mod ( D,p ) , f ( D,p ) and f k ( D,p ) . The average times forthe five examples are collected.For Algorithm 4.6, the relations between the running times and T, T , T are respectively givenin Figures 1, 2, and 3, where the parameter d is fixed. The relations between the running times and D, log D are respectively given in Figures 4 and 5, where the parameter t is fixed. Figures 4 and 5show that the complexity of Algorithm 4.6 is sensitive to D . Overall, these figures are basically inaccordance with the theoretical complexity bound O ∼ ( LT log D + LT log D ) for Algorithm 4.6.For Algorithm 5.8, the relations between the running times and T, T , T are respectively givenin Figures 6, 7, 8, where the parameters n, d are fixed. Similarly, the relations between the runningtimes and d, log d, n, n are respectively given in Figures Figures 9, 10, 11, 12. From these, wecan see that the practical performances are basically in accordance with the theoretical complexitybound O ∼ ( Ln T log D + LnT log D ) of Algorithm 5.8.For Algorithm 6.7, the relations between the running times and T, T are respectively givenin Figures 13, 14, where the parameters n, d, µ = 1 / D, log D, n, n are respectively given in Figures 15, 16, 17, 18.These figures show that the practical performance is worse than the theoretical complexity bound O ∼ ( LnT log D ) for Algorithm 6.7, because the logarithm factors omitted in the soft-Oh analysishave significant impact on the running time. In this paper, we give a new deterministic interpolation algorithm and two Monte Carlo interpola-tion algorithms for SLP sparse multivariate polynomials. Our deterministic algorithm has better20igure 1:
Horizontal-axis: T Figure 2:
Horizontal-axis: T Figure 3:
Horizontal-axis: T Figure 4: Horizontal-axis: d Figure 5: Horizontal-axis: log d Figure 6:
Horizontal-axis: T Figure 7:
Horizontal-axis: T Figure 8:
Horizontal-axis: T complexities than existing deterministic interpolation algorithms in most cases. Our Monte Carlointerpolation algorithms are the first algorithms whose complexities are linear in nT and polyno-mial in log D . The algorithms are based on several new ideas. In order to have a deterministicalgorithm, we give a criterion for checking whether a term belongs to a polynomial. We also give adeterministic method to find a “good” prime p in the sense that at least half of the terms in f arenot collisions in f mod ( D,p ) . Finally, a new Kronecker type substitution is given to reduce multivariatepolynomial interpolations to univariate polynomial interpolations.Figure 9: Horizontal-axis: d Figure 10: Horizontal-axis: ln d n Figure 12: Horizontal-axis: n Figure 13: Horizontal-axis: T Figure 14: Horizontal-axis: T Figure 15: Horizontal-axis: d Figure 16: Horizontal-axis: ln d Figure 17: Horizontal-axis: n Figure 18: Horizontal-axis: n eferences [1] N. Alon and Y. Mansour, “Epsilon-discrepancy sets and their application for interpolation ofsparse polynomials”, Inform. Process. Lett. 54(6), 337-342, 1995.[2] A. Arnold, “Sparse polynomial interpolation and testing,” PhD Thesis, Waterloo Unversity,2016.[3] A. Arnold and D.S. Roche “Multivariate sparse interpolation using randomized Kroneckersubstitutions,” Proc. ISSAC’14, ACM Press, 35-42, 2014.[4] A. Arnold, M. Giesbrecht, D.S. Roche, “Faster sparse interpolation of straight-line programs,”CASC 2013, LNCS 8136, 61-74, 2013.[5] A. Arnold, M. Giesbrecht, D.S. Roche, “Sparse interpolation over finite fields via low-orderroots of unity,” Proc. ISSAC’14, ACM Press, 2014.[6] A. Arnold, M. Giesbrecht, D.S. Roche, “Faster sparse multivariate polynomial interpolationof straight-line programs,” Journal of Symbolic Computation, 75, 4-24, 2016.[7] M. Avenda˜no, T. Krick, A. Pacetti, “Newton-Hensel Interpolation Lifting[J]”. Foundations ofComputational Mathematics, 2006, 6(1):82-120.[8] M. Ben-Or and P. Tiwari, “A deterministic algorithm for sparse multivariate polynomialinterpolation,” Proc. STOC’88 , 301-309, ACM Press, 1988.[9] M. Bl¨ a ser, M. Hardt, R.J. Lipton, N.K. Vishnoi, “Deterministically testing sparse polynomialidentities of unbounded degree,” Information Processing Letters, 109, 187-192, 2009.[10] D.G. Cantor and E. Kaltofen, “On fast multiplication of polynomials over arbitrary algebras.”Acta Informatica 28.7(1991):693-701.[11] S. Garg and E. Schost, “Interpolation of polynomials given by straight-line programs,” The-oretical Computer Science, 410, 2659-2662, 2009.[12] M. Giesbrecht and D.S. Roche, “Diversification improves interpolation,” ISSAC’11, ACMPress, 123-130, 2011.[13] Q.L. Huang, X.S. Gao, “Sparse interpolation of black-box multivariate polynomials usingkronecker type substitutions,” arXiv 1710.01301, 2017.[14] E. Kaltofen, “Computing with polynomials given by straight-line programs I: greatest commondivisors.” Seventeenth ACM Symposium on Theory of Computing ACM, 1985:131-142.[15] E. Kaltofen, “Greatest common divisors of polynomials given by straight-line programs.”Journal of the Acm 35.1(1988):231-264.[16] E. Kaltofen and Y.N. Lakshman, “Improved sparse multivariate polynomial interpolationalgorithms,” Proc. ISSAC’88, 467-474, ACM Press, 1988.[17] A.R. Klivans and D. Spielman, “Randomness efficient identity testing of multivariate polyno-mials,” Proc. STOC01, 216-223, ACM Press, 2001.2318] L. Kronecker, “Grundz¨ u ge einer arithmetischen theorie der algebraischen Gr¨ o ssen,” Journalf¨ uu