Faster interpolation algorithms for sparse multivariate polynomials given by straight-line programs

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.1 (1-20)

Journal of Symbolic Computation ••• (••••) •••–•••

Contents lists available at ScienceDirect

Journal of Symbolic Computation www.elsevier.com/locate/jsc

Faster interpolation algorithms for sparse multivariate polynomials given by straight-line programs Qiao-Long Huang a,b , Xiao-Shan Gao a,b,∗ a b

KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, China University of Chinese Academy of Sciences, China

a r t i c l e

i n f o

Article history: Received 5 June 2018 Accepted 26 September 2019 Available online xxxx MSC: 00-01 99-00 Keywords: Black box Interpolation Kronecker substitution

a b s t r a c t In this paper, we propose new deterministic and Monte Carlo interpolation algorithms for sparse multivariate polynomials represented by straight-line programs. Let f be an n-variate polynomial given by a straight-line program, which has a total degree bound D and a term bound T . Our deterministic algorithm is quadratic in n, T and cubic in log D in the Soft-Oh sense, which has better complexities than existing deterministic interpolation algorithms in most cases. Our Monte Carlo interpolation algorithms have better complexities than existing Monte Carlo interpolation algorithms and are the first algorithms whose complexities are linear in nT and polynomial in log D in the Soft-Oh sense. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The sparse interpolation for multivariate polynomials has received considerable interest. There are two basic models for this problem: the polynomial is either given as a straight-line program (SLP) (Avendaño et al., 2006; Garg and Schost, 2009; Giesbrecht and Roche, 2011; Kaltofen, 1985) or a more general black-box (Alon and Mansour, 1995; Ben-Or and Tiwari, 1988; Huang and Gao, 2019;

*

Corresponding author. E-mail address: [email protected] (X.-S. Gao). URL: http://www.mmrc.iss.ac.cn/~xgao/ (X.-S. Gao).

https://doi.org/10.1016/j.jsc.2019.10.005 0747-7171/© 2019 Elsevier Ltd. All rights reserved.

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.2 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

2

Table 1 A “soft-Oh” comparison for SLP polynomials over an arbitrary ring R. Algorithms

Total Cost

Type

Dense Garg and Schost (2009) Randomized (Giesbrecht and Roche, 2011) Arnold et al. (2013)

L Dn Ln2 T 4 log2 D Ln2 T 3 log2 D Ln3 T log3 D

Deterministic Deterministic Las Vegas Monte Carlo

This paper (Theorem 5.8) This paper (Theorem 6.8)

Ln2 T 2 log2 D + LnT log3 D LnT log3 D

Deterministic Monte Carlo

Kaltofen and Lakshman, 1988; Zippel, 1990). In this paper, we consider the problem of interpolation for a sparse multivariate polynomial given by an SLP. 1.1. Main results Let f ∈ R[x1 , . . . , xn ] be an SLP polynomial of length L, with a total 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 logc for some fixed c > 0. When we say “linear”, “cubic”, etc., we mean linear and cubic in the sense of “Soft-Oh” complexity. We propose a new deterministic interpolation algorithm for f given by an SLP, whose complexity is O ∼ ( Ln2 T 2 log2 D + LnT log3 D ) operations in R and a similar number of extra bit operations. We also propose two new Monte Carlo interpolation algorithms for SLP multivariate polynomials. For a 1 given μ ∈ (0, 1), the complexity of our first algorithm is O ∼ ( LnT (log3 D + log D log μ )) operations in R, and with probability at least 1 − μ, it returns the correct polynomial. In our second algorithm, R is a finite field Fq and we can evaluate f in a proper extension field of Fq . The bit complexity of 1 our second algorithm is O ∼ ( LnT log2 D (log D + log q) log μ ), and with probability at least 1 − μ, it returns the correct polynomial. These algorithms are the first ones whose complexity is linear in nT and polynomial in log D in the Soft-Oh sense. In Table 1, we list the complexities for the SLP interpolation algorithms which are similar to the methods proposed in this paper. In the table, L is the size of the SLP and “total cost” means the number of operations in R. For deterministic algorithms, the ratio of the cost of our deterministic nT +log D method with that of the algorithm given in Garg and Schost (2009) is . So our method has nT 3



better complexity unless f is super sparse, more precisely, unless T <





3

log D n

in the Soft-Oh sense.

d+n terms, our algorithm works better in d most cases. For probabilistic algorithms, our method is the only one whose complexity is linear in nT and polynomial in log D and is better than that of the Monte Carlo method given in Arnold et al. (2013). Kaltofen gave an interpolation algorithm for SLP polynomials (Kaltofen, 1985), whose complexity is polynomial in D. Avendaño et al. (2006) gave an algorithm for interpolating an SLP f ∈ Z[x], with bit complexity polynomial in L , log D , h and h , where h is an upper bound on the height of f and h an upper bound on the height of the values of f (and some of its derivatives) at some sample points. This method does not seem to extend to arbitrary rings. Mansour (1995) and Alon and Mansour (1995) gave deterministic algorithms for polynomials in Z[x], with bit complexity polynomial in 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-box polynomials (Ben-Or and Tiwari, 1988; Kaltofen and Lakshman, 1988; Zippel, 1990; Klivans and Spielman, 2001) can also be used to SLP polynomials, but the complexities of these algorithms are polynomial in D instead of log D. In Table 2, we list the bit complexities for SLP interpolation algorithms over finite field Fq . From Table 2, we can see that our probabilistic algorithms are the only ones which are linear in n, T and Noting that a dense polynomial with degree d has

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.3 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

3

Table 2 A “soft-Oh” comparison for SLP polynomials over finite field Fq . Algorithms

Bit Complexity

Algorithm type

Garg and Schost (2009) Randomized Garg-Schost (Giesbrecht and Roche, 2011) Giesbrecht and Roche (2011) Arnold et al. (2013) Arnold et al. (2014) Arnold et al. (2016)

Ln2 T 4 log2 D log q Ln2 T 3 log2 D log q Ln2 T 2 log2 D (n log D + log q) Ln3 T log3 D log q LnT log2 D (log D + log q) + nω T Ln log D ( T log D + n)(log D + log q) +nω−1 T log D + nω log D

Deterministic Las Vegas Las Vegas Monte Carlo Monte Carlo Monte Carlo

Ln2 T 2 log2 D log q + LnT log3 D log q

Deterministic Monte Carlo Monte Carlo

This paper (Theorem 5.8) This paper (Theorem 6.8) This paper (Theorem 6.11)

LnT log3 D log q LnT log2 D (log q + log D )

polynomial in log D. Also, our second probabilistic algorithm has better complexity than all algorithms except randomized Garg-Schost (Giesbrecht and Roche, 2011) which is Las Vegas. The ratio of the cost of our first probabilistic algorithm with that of the randomized Garg-Schost (Giesbrecht and Roche, log D 2011) is nT 2 , so our method is faster in most cases. 1.2. Main idea and relation with existing work Our methods build on the work by Garg and Schost (2009), Arnold et al. (2013), Giesbrecht and Roche (2011), and Klivans and Spielman (2001), where the basic idea is to reduce multivariate interpolation to univariate interpolation. Three new techniques are introduced in this paper: a criterion for checking whether a term belongs to a polynomial (see Section 2 for details), a deterministic method to find an “ok” prime (see Section 3 for details), a new Kronecker type substitution to reduce multivariate polynomial interpolation to univariate polynomial interpolation (see Section 5.1 for details). Our methods 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 merge D D n−1 with other terms of f in the univariate polynomial f (mod ) mod (x p − 1). D , p ) = f (x, x , . . . , x • Second, we obtain a set S of terms containing those non-colliding terms of f . In the univariate case, these terms are found by the Chinese Remaindering Theorem and in the multivariate case, 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 at least half of the terms of f from S. Repeating these three steps for at most log2 T times, we obtain f . In the rest of this section, we give detailed comparison with related work. Garg and Schost (2009) gave a deterministic interpolation algorithm for a univariate SLP polynomial f by recovering f from f mod (x p − 1) for O ( T 2 log D ) different primes p. The randomized Las Vegas version of this method needs O ( T log D ) probes. The multivariate interpolation comes directly from the Kronecker substitution (Kronecker, 1882). Our univariate interpolation algorithm has two major differences from that given in Garg and Schost (2009). First, we compute f mod (x p − 1) for O( T log D ) different primes p, and second we introduce a criterion to check whether a term really belongs to f . Our multivariate interpolation method is similar to our univariate interpolation algorithm, where a new Kronecker type substitution is introduced to recover the exponents. Giesbrecht and Roche (2011) introduced the idea of diversification and a probabilistic method to choose “good” primes. It improves Garg and Schost’s algorithm by a factor O ( T 2 ), but becomes a Las Vegas algorithm. In our Monte Carlo algorithms, we use a new Kronecker substitution instead of the method of diversification to find the same term in different remainders of f , where only additions are used for the coefficients. Hence, our algorithm can work for more general rings and has better complexity.

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.4 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

4

In Arnold et al. (2013), the concept of “ok” prime is introduced and a Monte Carlo univariate algorithm is given, which has complexity linear in T but cubic in n. The “ok” prime in Arnold et al. (2013) is probabilistic. In our deterministic method, we give a method to find a guaranteed “ok” prime. In Arnold et al. (2014), their univariate interpolation algorithm is extended to finite fields. By combining the idea of diversification, the complexity becomes better. This algorithm will be used in our second probabilistic algorithm. The bit complexity of their multivariate interpolation algorithm is linear in nω , where ω is the constant of matrix multiplication, while our algorithm is linear in n. The reason is that, our method uses a new Kronecker substitution to find the exponents and does not need to solve linear systems. In Arnold et al. (2016), they further improved their interpolation algorithm for finite fields. By combining the random Kronecker substitution and diversification, the complexity becomes better, but still linear in nω . Finally, the new Kroneceker type substitution introduced in this paper is inspired by the work of Klivans and Spielman (2001), Arnold (2016) and Arnold and Roche (2014). In Klivans and Spielman n −1 (2001), the substitution f (q1 x, q2 xmod( D , p ) , . . . , qn xmod( D , p ) ) was used, where q1 , . . . , qn are primes. k −1

n−1

In this paper, we introduce the substitution f (x, xmod( D , p ) , . . . , x p +mod( D , p ) , . . ., xmod( D , p ) ) (see section 5.1 for exact definition). Our substitution has the following advantages: (1) For the field of complex numbers, the size of data is not changed after our substitution, while the size of data for the substitution in Klivans and Spielman (2001) is increased by a factor of D. (2) Only arithmetic operations for the coefficients are used in our algorithm and thus the algorithm works for general computable rings, while the substitution in Klivans and Spielman (2001) needs factorization and R should be a UFD at least. In Arnold (2016), the substitution f (xs1 , . . . , xsk + p . . . , xsn ) was used, where si are random integers. Comparing to the randomized Kronecker substitution in Arnold, Arnold and Roche (2016, 2014), our substitution is deterministic. 2. A criterion for term testing

In this section, we give a criterion to check whether a term belongs to a polynomial. Throughout this paper, let f = c 1 m1 + c 2 m2 + · · · + ct mt ∈ R[X] be a multivariate polynomial with terms c i mi , where R is a computable ring, X = {x1 , x2 , . . . , xn } are n indeterminates, and mi , i = 1, 2, . . . , t are distinct monomials. Denote # f = t to be the number of terms of f , deg( f ) to be the total degree of f , and M f = {c 1 m1 , c 2 m2 , . . . , ct mt } to be the set of terms of f . Let D , T ∈ N such that D > deg( f ) and T ≥ # f . For p ∈ N>0 , let D D f (mod D , p ) = f (x, x , . . . , x

n −1

) mod (x p − 1) ∈ 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 mod that mmod ( D , p) = w ( 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 for any prime q, mod cm is also not a collision in f ( D , pq) .

t

Lemma 2.3. Let f = i =1 c i mi , T ≥ # f , D > deg( f ), N 1 = max{1, n( T − 1) log2 D }. For each cm ∈ M f , there exist at most N 1 − 1 primes p 1 , . . . , p N 1 −1 such that cm is a collision in f (mod D , p ) for all i = 1, . . . , N 1 − 1. i

Proof. If T = 1, then N 1 = 1. Since f has only one term, f has no collision in f (mod D , p ) for any prime p and the lemma is correct. Now we assume T ≥ 2, then N 1 = n( T − 1) log2 D . It suffices to show

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.5 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

5

that for any N 1 different primes p 1 , p 2 , . . . , p N 1 , there exists at least one p j , such that cm is not a collision in f (mod D,p ) . j

e i ,1 e i ,2

e i ,n

Assume mi = x1 x2 · · · xn , i = 1, 2, . . . , t. We prove it by contradiction. It suffices to consider the case of c 1 m1 . We assume by contradiction that for every p j , j = 1, 2, . . . , N 1 , c 1 m1 is a collision in f (mod D , p ) . Let j

B=

t  n n   ( e 1 , i D i −1 − e s,i D i −1 ). s =2 i =1

i =1

First, we show that if c 1 m1 is a collision in f (mod D , p j ) , then mod( B , p j ) = 0. Since c 1 m1 is a colmod mod , without loss of generality, assume mmod lision in f (mod 1( D , p ) = m2( D , p ) . Then 0 = deg(m1( D , p ) ) − D,p ) j

n

deg(mmod 2( D , p ) ) = mod(

−

n

j

i =1 e 2,i D

i −1

n

i −1 , p j ) − mod( i =1 e 1,i D

j

n

j

i −1 , p j ). So we have mod( i =1 e 2,i D

j

i −1 i =1 e 1,i D

, p j ) = 0. So mod( B , p j ) = 0.

N1 n p divides B. Note that | i =1 e 1,i D i −1 − j =1 j     i −1 | ≤ ( D − 1)( ni=1 D i −1 ) = D n − 1. So | B | = | ts=2 ( ni=1 e 1,i D i −1 − ni=1 e s,i D i −1 )| ≤ i =1 e s,i D  n t −1 ( D − 1) . Thus Nj =1 1 p j ≥ 2N 1 ≥ 2n(T −1) log2 D = D n(T −1) > ( D n − 1)T −1 ≥ | B |, which contradicts the N1 fact that p divides B. The lemma is proved. 2 j =1 j Since p 1 , p 2 , . . . , p N 1 are different primes,

n

Now we give a criterion for testing whether a term cm is in M f .

t

Theorem 2.4. Let f = i =1 c i mi , T ≥ # f , D > deg( f ), N 1 = max{1, n( T − 1) log2 D }, N 2 = max{1, nT log2 D }, and P = { p 1 , p 2 , . . . , p N 1 + N 2 −1 } be N 1 + N 2 − 1 different primes. For a term cm satisfying deg(m) < D, cm ∈ M f if and only if there exist at least N 2 integers j ∈ [1, N 1 + N 2 − 1] such that mod #( f − cm)mod (D,p ) < # f (D,p ) . j

j

Proof. If T = 1, then N 1 = 1, N 2 ≥ 1 and N 1 + N 2 − 1 = N 2 . In this case, cm ∈ M f if and only if mod #( f − cm)mod ( D , p ) = 0 and # f ( D , p ) = 1 for any prime p j , so the theorem is correct. Now we assume j

j

T ≥ 2, then N 1 = n( T − 1) log2 D . mod 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)( D , p j ) = mod mod # f (mod D , p ) − 1. So #( f − cm)( D , p ) < # f ( D , p ) . By Lemma 2.3, there exist at most N 1 − 1 primes q j such j

j

j

that cm is a collision in f (mod D ,q j ) . In P , as N 1 + N 2 − 1 − ( N 1 − 1) = N 2 , there exist at least N 2 primes mod such that #( f − cm)( D , p ) < # f (mod D,p ). j

j

/ M f . We show there exist at most N 2 − 1 integers j ∈ [1, N 1 + For the other direction, assume cm ∈ N 2 − 1] such that #( f − cm)mod < # f (mod (D,p ) D , p ) . Consider two cases: Case 1: m is not a monomial in f . j

j

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. By Lemma 2.3, there exist at most N 2 − 1 primes in P such that −cm is a collision in f − cm. For all other mod mod mod primes p j in P , #(( f − cm) −(−cm))mod ( D , p ) = #( f − cm)( D , p ) − 1, that is, #( f − cm)( D , p ) = # f ( D , p ) + 1. j

j

j

j

mod So there exist at most N 2 − 1 primes p j in P such that #( f − cm)mod (D,p j ) < # f (D,p j ) . / M f , f − cm has the same number of terms as Case 2. Since m is a monomial in f and cm ∈ f . Assume the term of f with monomial m is c 1 m. Then (c 1 − c )m ∈ M f −cm . By Lemma 2.3, for at most N 1 − 1(≤ N 2 − 1) primes p j in P , (c 1 − c )m is a collision in ( f − cm)mod ( D , p ) . For all other j

mod primes pk in P , (c 1 − c )m is not a collision in ( f − cm)mod ( D , pk ) , or equivalently, #( f − cm)( D , pk ) = mod mod mod #( f − cm − (c 1 − c )m)( D , p ) + 1 = #( f − c 1 m)( D , p ) + 1. But we always have #( f − c 1 m)( D , p ) + 1 ≥ k

k

k

mod mod # f (mod D , pk ) , and so #( f − cm)( D , pk ) ≥ #( f )( D , pk ) . So there exist at most N 2 − 1 primes p j in P such that mod mod #( f − cm)( D , p ) < # f ( D , p ) . The theorem is proved. 2 j

j

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.6 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

6

As a corollary, we can deterministically recover f from f (mod D , p j ) , j = 1, 2, . . . , N 1 + N 2 − 1. Corollary 2.5. Use the notations in Theorem 2.4. We can uniquely recover f from f (mod D , p j ) , j = 1, 2, . . . , N 1 + N 2 − 1. Proof. Let C = {c 1 , c 2 , . . . , ck } be all the different coefficients in f (mod D , p ) , j = 1, 2, . . . , N 1 + N 2 − 1. By j

Lemma 2.3, since N 1 + N 2 − 1 ≥ N 1 , all the coefficients of f are in C . Let M = {m1 , m2 , . . . , ms } 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, 2, . . . , k, j = 1, 2, . . . , s}. By Theorem 2.4, we can check if c i m j is in M f . So we can find all the terms of f . 2 The above result can be changed into a deterministic algorithm for interpolating f . But the algorithm is not efficient due to the reason that s is linear in D n . In the following, we will show how to find a smaller alternative set M and give an efficient interpolation algorithm. 3. Find an “ok” prime 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 . f

Denote C( D , p ) to be the number of collision terms of f in f (mod D , p ) . We need the following lemma from Arnold et al. (2013). f

f

mod Lemma 3.1. (Arnold et al., 2013) Let f ∈ R[x]. If # f (mod D ,q) ≤ # f ( D , p ) , then C( D , p ) ≤ 2C( D ,q) .

It is easy to modify the above lemma into multivariate case. f

f

mod Corollary 3.2. Let f ∈ R[X], D > deg( f ). If # f (mod D ,q) ≤ # f ( D , p ) , then C( D , p ) ≤ 2C( D ,q) .

Proof. Let F = f (x, x D , . . . , x D follows from Lemma 3.1. 2

e j ,k ) D

, and p a

f

F p ), then f (mod D , p ) = F mod (x − 1). So C( D , p ) = C( D , p ) . The corollary

t

e i ,1 e i ,2 e i ,n i =1 c i mi ∈ R[X], mi = x1 x2 · · · xn , s f prime. If C( D , p ) = s, then p  2  divides A.

Lemma 3.3. Let f = k−1

n −1

D > deg( f ), A =

i < j

i , j ∈{1,...,t }

n

k=1 (e i ,k

−

Proof. We divide the terms of f into groups, called collision blocks, such that two terms of f collide in f (mod D , p ) if and only if they are in the same group. Let ni be the number of collision blocks containing i terms. Assume c j 1 m j 1 + c j 2 m j 2 + · · · + c j i m j i is the sum of all the terms in a collision block with mod i terms. For any u , v ∈ { j 1 , j 2 , . . . , j i }, we have deg(mmod u ( D , p ) ) = deg(m v ( D , p ) ). So (e u ,1 + e u ,2 D + · · · +

e u ,n D n−1 ) mod p = (e v ,1 + e v ,2 D + · · · + e v ,n D n−1 ) mod p, which implies that p divides (e u ,1 − e v ,1 ) + (e u ,2 − e v ,2 ) D + . . . + (e u ,n − e v ,n ) D n−1 .

t

i (i −1)

i (i −1)

There are pairs such u , v, so p 2 is a factor of A. Let K = i =1 21 (i 2 − i )ni . Since there 2 exist ni  such collisionblocks, p K is a factor of A. Now  we give a lower t bound ofKt . First we see t t t t that t = i =1 ini , s = i =2 ini . K = i =1 21 (i 2 − i )ni = 12 i =1 i 2 ni − 12 i =1 ini = 12 i =1 i 2 ni − 12 t = 1 n 2 1

+

1 2

t

i =2 i

2

ni − 12 t ≥ 12 n1 + t − n1 − 12 t = 12 t − 12 n1 =

there is at least one ni > 0, i > 3, then K >

1 s. 2

So K ≥

1 s. If 2 1  2 s. We

=  12 s. If have proved the lemma. 2 ni = 0, i ≥ 3, then K =

1 s 2

t

Theorem 3.4. Let f = i =1 c i mi ∈ R[X], T ≥ # f , D > deg( f ), N 1 = max{1, n( T − 1) log2 D }, and mod p 1 , p 2 , . . . , p 4N 1 be 4N 1 different primes. Let j 0 be an integer in [1, 4N 1 ] such that # f (mod D,p ) ≥ # f (D,p ) for all j. Then at least  2t  of the terms of f are not collisions in f (mod D , p j0 ) .

j0

j

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.7 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

7

Proof. If T = 1, then f has only one term, and f has no collision in f (mod D , p ) for any prime p and the theorem is correct. So we assume T ≥ 2, then N 1 = n( T − 1) log2 D . We first claim that there exists f

at least one p j in p 1 , p 2 , . . . , p 4N 1 such that C( D , p ) < j

t . 4

We prove it by contradiction. Assume for f

f

j = 1, 2, . . . , 4N 1 , C( D , p ) ≥ j

 12 C( D , p ) 

t . 4

Then by Lemma 3.3, for all j in [1, 4N 1 ], p j

j

divides A, where

A is defined in Lemma 3.3. Since p j , j = 1, 2, . . . , 4N 1 are different primes, then divides A. Now

4N 1

j =1

f

 12 C( D , p ) 

pj

j

≥

4N 1

j =1

contradicts to the inequality A ≤ ( D n − 1) By Corollary 3.2, we have

f C( D , p j ) 0

number of no collision terms of f in

 1 · 14 t 

pj2

t (t −1) 2

t

1

1

f

 12 C( D , p ) 

p j =1 j

j

1

≥ 24N 1 8 ≥ 2 2 t N 1 ≥ 2 2 tn(T −1) log2 D = D 2 nt (T −1) , which

. We proved the claim.

f

f

≤ 2C( D , p j ) , j = 1, 2, . . . , 4N 1 . So C( D , p j

f (mod D , p j0 )

4N 1

is > t −

t 2

=

1 t. 2

0

)

<2·

t 4

We proved the theorem.

= 12 t. So the 2

4. Deterministic univariate interpolation In this section, we consider the interpolation of a univariate polynomial f with deg( f ) < D. The algorithm works as follows. First, we use Theorem 3.4 to find an “ok” prime p such that at least half mod of the terms of f are not collisions in f (mod D , p ) . Second, we use f ( D , ppk ) , k = 1, 2, . . . , K D to find a set S containing these non-collision terms of f by the Chinese Remaindering Theorem, where pk is the k-th prime and K D is the smallest number such that p 1 p 2 . . . p K D ≥ D. Finally, we use Theorem 2.4 to pick up the terms of f from S. 4.1. Recovering terms from module residue In this section, let f be a univariate polynomial in R[x]. We will give an algorithm to recover mod those terms of f from f (mod D , p ) , which are not collisions in f ( D , p ) . 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 − 1). Write d1 d2 dr f (mod D , p ) = a 1 x + a 2 x + · · · + ar x

f (mod D , ppk )

= f k,1 + f k,2 + · · · + f k,r + gk ,

(2) k = 1, 2, . . . , K D

where d1 < · · · < dr , pk is the k-th prime, K D is the smallest number such that p 1 p 2 . . . p K D ≥ D, and f k,i mod (x p − 1) = ai xdi and gk mod (x p − 1) = 0. f (mod D , pp ) can be written as the above form, because k

mod p f (mod D , ppk ) mod (x − 1) = f ( D , p ) . We now introduce the following key notation

f

U D , p = {ai xe i |i ∈ [1, r ], ai is from (2), e i ∈ [0, D − 1], and U1 : f k,i = ai xbk,i , k = 1, 2, . . . , K D .

(3)

U2 : mod(e i , pk ) = mod(bk,i , pk ), k = 1, 2, . . . , K D .} f

The following lemma gives the geometric meaning of U D , p . f

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 D , p . f

Proof. It suffices to show that cm satisfies the conditions of the definition of U D , p . Assume m = xe . mod di Since cm is not a collision in f (mod for some i and di = mod(e , p ), where D , p ) , assume cm( D , p ) = ai x bk ,i ai xdi is defined in (2). By Lemma 2.2, cm is also not a collision in f (mod for bk,i = D , ppk ) and f k,i = ai x mod(e , ppk ), k = 1, 2, . . . , K D . Since mod(e , pk ) = mod(mod(e , ppk ), pk ) = mod(bk,i , pk ), conditions U1 and U2 are satisfied and the lemma is proved. 2

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.8 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

8 f

f

Note that U D , p may contain terms not in M f . The following algorithm computes the set U D , p . Algorithm 4.2 (UTerms). mod Input: Univariate polynomials f (mod D , p ) , f ( D , ppk ) , primes p k , for k = 1, 2, . . . , K D ; a prime p; a degree bound D > deg( f ). f

Output: U D , p . mod Step 1: Write f (mod D , p ) and f ( D , ppk ) as in form (2):

d1 d2 dr f (mod D , p ) = a 1 x + a 2 x + · · · + ar x

f (mod D , ppk ) = f k,1 + f k,2 + · · · + f k,r + gk where k = 1, . . . , K D . Step 2: Let U = {}. For i = 1, 2, . . . , r a: If for k = 1, 2, . . . , K D , one of #( f k,i ) = 1 or one of the coefficients of f k,i is not ai , then break. b: For k = 1, 2, . . . , K D , assume f k,i = ai xbk,i and let ek = mod(bk,i , pk ). c: Let β = Chinese Remainder ([e 1 , e 2 , . . . , e K D ], [ p 1 , p 2 , . . . , p K D ]). d: If β < D then let U = U {ai xβ }. Step 3: Return U . Remark 4.3. In c of Step 2, Chinese Remainder([e 1 , . . . , e K D ], [ p 1 , . . . , p K D ]) means to find an integer 0 ≤ w < D such that mod( w , p i ) = e i , i = 1, 2, . . . , K D . Since p 1 p 2 · · · 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 ) operations in R and O ∼ ( T log D log p ) bit operations. Proof. In Step 1, we need to do a traversal for the terms of f (mod D , ppk ) . Since K D is the smallest number such that p 1 p 2 · · · p K D ≥ D, then 2 K D −1 ≤ p 1 p 2 · · · p K D −1 < D. So K D ≤ log2 D + 1, which is O (log D ). In order to write f (mod D , pp ) as the form (2), we need to perform the modular operation mod p on k

every exponent of f (mod D , ppk ) , then use the merge sort method to write their terms in ascending order according to the degree. Since f (mod D , ppk ) has no more than T terms and k = 1, 2, . . . , K D , it needs O ∼ ( T log D ) arithmetic operations in Z. 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 ) bit operations. ∼ In a of Step 2, since # f (mod D , ppk ) ≤ T , it totally needs O ( T log D ) bit operations to determine whether #( f k,i ) = 1. To compare the coefficients of f k,i , it needs O ( T log D ) arithmetic operations in R. In b of Step 2, since the height of the data is O (log( pp K D )), it needs O ∼ ( T log D log p ) bit operations. In c, for each computing of β , by von zur Gathen and Gerhard (2013, Theorem 10.25), the cost of the Chinese remaindering algorithm is O ∼ (log D ) bit operations, since β is O ( D ). Since we call the Chinese remaindering at most T times, it needs O ∼ ( T log D ) bit operations. So the complexity of Step 2 is O ∼ ( T log D log p ) bit operations and O ( T log D ) arithmetic operations in R. 2 4.2. Interpolation algorithm for univariate polynomials 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 that represents a polynomial function. It takes as input a vector (a1 , . . . , an ) and outputs a vector (b1 , . . . , b L )

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.9 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

9

by way of a series of instructions ≤ L of the form i : bi ← α i β , where i is an operation i : 1 ≤ i  + ,  − or  × , and α , β ∈ R {a , . . . , a } {b , . . . , b 1 n 0 i −1 }. The inputs and outputs may belong to R or a ring extension of R. We say that an SLP computes a multivariate polynomial f ∈ R[x1 , . . . , xn ] if it sets b L to be f (a1 , . . . , an ). We now give the interpolation algorithm for univariate polynomials, where f ∗ in the input is introduced 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 1 ≥ #( f − f ∗ ), ( T ≥ T 1 ), D > max{deg( f ), deg( f ∗ )}. Output: The exact form of f − f ∗ . Step 1: Step 2: Step 3: Step 4:

Let N 1 = max{1, ( T 1 − 1) log2 D }, N 2 = max{1,  T 1 log2 D }, N = max{4N 1 , N 1 + N 2 − 1}. Find the first N primes p 1 , . . . , p N . Compute the smallest K D such that p 1 · · · p K D ≥ D. mod ∗mod ∗ For j = 1, 2, . . . , N, probe f (mod D , p ) . Let f j = f ( D , p ) − f ( D , p ) and h = 0. j

j

j

Step 5: Loop 5.1: Let α = max{# f j | j = 1, 2, . . . , N } and j 0 the smallest number such that # f j 0 = α . 5.2: If α = 0, then return h∗ . mod ∗mod ∗mod 5.3: For k = 1, 2, . . . , K D , probe f (mod D , p p ) and let gk = f ( D , p p ) − f ( D , p p ) − h ( D , p p ) . f − f ∗ −h∗

5.4: Let U D , p j0

j0

j0

k

k

j0

k

j0

k

:= UTerms( f j 0 , g 1 , g 2 , . . . , g K D , p j 0 , D ). f − f ∗ −h∗

5.5: Let h = 0. For each u ∈ U D , p

j0

, if

#{ j | #( f j − u mod ( D , p j ) ) < #( f j ), j = 1, . . . , N 1 + N 2 − 1} ≥ N 2 then h := h + u. 5.6: Let h∗ = h∗ + h, T 1 = T 1 − #h, N 1 = max{1, ( T 1 − 1) log2 D }, N 2 = max{1,  T 1 log2 D }, N = max{4N 1 , N 1 + N 2 − 1}. 5.7: For j = 1, 2, . . . , N, let f j = f j − hmod (D,p ) . j

Theorem 4.7. Algorithm 4.6 returns f − f ∗ using O ∼ ( LT 2 log2 D + LT log3 D ) 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 obtain at least half of the terms of f − f ∗ − h∗ . Then, the algorithm will return the correct f − f ∗ by running at most log2 T 1 times of the loop in Step 5. In Step 5.1, Theorem 3.4 is used to find an okay prime p j 0 . f − f ∗ −h∗ In Step 5.4, by Lemma 4.1 and Theorem 3.4, at least half of the terms of f − f ∗ − h∗ are in U D , p . j0 f − f ∗ −h∗

In Step 5.5, Theorem 2.4 is used to select the elements of M f − f ∗ −h∗ from U D , p . In summary, j0 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 analyze the complexity of the algorithm, which comes from Step 4 and Step 5. The complexity of other steps is lower than that in these two steps. In Step 2, since the bit complexity of finding the first N primes is O ∼ ( N ) by von zur Gathen and Gerhard (2013, p. 500, Thm.18.10) and N is O ∼ ( T log D ), the bit complexity of Step 2 is O ∼ ( T log D ). mod ∼ In Step 4, we probe N univariate polynomials f (mod D , p ) and probing f ( D , p ) costs O ( Lp j ), since j

j

S f is of length L and the univariate polynomials in the procedure have degrees < p j . Since p j is

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.10 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

10

2 ∼ 2 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 ) operations in R plus similarly many bit operations. It needs O ∼ ( T 2 log D ) ring operations and O ∼ ( T 2 log2 D ) bit 2 ∼ 2 operations to obtain f (∗Dmod , p ) . Then the total complexity of Step 4 is O ( LT log D ) operation in R j

plus similarly many bit operations. We now consider Step 5. We first consider the complexity of each loop of this step. In Step 5.1, ∼ 2 since N is of O ( T log D ) and the terms of ( f − f ∗ − h∗ )mod ( D , p ) is no more than T , it needs O ( T log D ) j

bit operations. In Step 5.3, since K D is of O (log D ) and p j 0 pk is of O ∼ ( T log2 D ), we need O ∼ ( LT log3 D ) arith∗ metic operations in R and similarly many bit operations to obtain f (mod D , p p ) . We need O ( K D # f ) j0

k

ring operations in R and O ∼ ( K D # f ∗ max{log deg( f ∗ ), log( p j 0 pk )}) bit operations to obtain f (∗Dmod , p j 0 pk ) . Since # f ∗ ≤ T , deg( f ∗ ) ≤ D, pk is of O (log D ), and p j 0 is of O ( T log D ), the cost is O ( T log D ) operations in R and O ∼ ( T log D max{log D , log T + log log D }) = O ∼ ( T log2 D ) bit operations. Since #h∗ ≤ T 1 and deg(h∗ ) < D, it also needs O ( T log D ) operations in R and O ∼ ( T log2 D ) bit operations 3 ∼ to obtain h∗( Dmod , p p ) . Then, the total complexity of this step is O ( LT log D ) arithmetic operations in j0

k

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 log q) bit operations. In Step 5.5, in order to determine whether #( f j − u mod ( D , p ) ) < #( f j ), we just need to determine j

whether u mod ( D , p j ) is a term of f j . We sort the terms of f j such that they are in ascending order according to their degrees, which costs O ∼ (( N 1 + N 2 ) T ) = O ∼ ( T 2 log D ) bit operations, since N 1 + N 2 is O ∼ ( T log D ). To find whether f j has a term with degree deg(u mod ( D , p ) ), we need O (log T ) comparisons. j

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 in R. So it totally needs O (log T log D ) bit operations and O (1) arithmetic operation in R to compare #( f j − u mod ( D , p ) ) with #( f j ). Hence, the total complexity of j

f−f∗

f−f∗

Step 5.5 is O ∼ (#U D , p ( N 1 + N 2 ) log T log D + ( N 1 + N 2 ) T ) bit operations and O (#U D , p ( N 1 + N 2 )) j0

f−f∗

arithmetic operations in R. Since #U D , p

j0

j0

≤ T and N 1 + N 2 − 1 is of O ( T log D ), the total complexity

is O ∼ ( T 2 log D ) arithmetic operations in R and O ∼ ( T 2 log2 D ) bit operations. So the total complexity of each loop of Step 5 is O ∼ ( LT log3 D + T 2 log D ) arithmetic operations in R and O ∼ ( LT log3 D + T 2 log2 D ) bit operations, which comes from Step 5.3 and Step 5.5. 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 log3 D + T 2 log D ) operations in R and O ∼ ( LT log3 D + T 2 log2 D ) bit operations. Combining with the complexity of Step 4, the total complexity of the algorithm is O ∼ ( LT 2 log2 D + LT log3 D + T 2 log2 D ) = O ∼ ( LT 2 log2 D + LT log3 D ) operations in R and similarly many bit operations, which comes from Step 4 and Step 5.3. The theorem is proved. 2

For an n-variate polynomial f of degree < D, we can use the Kronecker substitution (Kronecker, 1882) to reduce the interpolation of f to that of a univariate polynomial of degree D n , which can be computed with Algorithm UIPoly( S f , 0, T , T , D ). By Theorem 4.7, we have

Corollary 4.8. For an SLP S f to compute f ∈ R[X] with T ≥ # f and D < deg( f ), we can find f using O ∼ ( Ln2 T 2 log2 D + Ln3 T log3 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 has better complexity.

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.11 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

11

5. Deterministic multivariate polynomial interpolation 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 Algorithm 4.6 and works as follows. First, we use Theorem 3.4 to find an “ok” prime p for f . Second, we use a modified Kronecker substitution to obtain a set S of terms, which contains at least half of the terms of f . Finally, we use Theorem 2.4 to identify the terms of f from S. The multivariate interpolation algorithm will call Algorithm 4.6. 5.1. Recovering terms from module residue Let f be a multivariate polynomial, M f the set of terms in f , t = # f , D > deg( f ), T ≥ # f , and p ∈ N>0 . Consider the modified Kronecker substitutions:

f ( D , p ) = f (x, xmod( D , p ) , . . . , xmod( D mod( D , p )

f ( D , p ,k) = f (x, x

,...,x

n −1

, p)

(4)

)

p +mod( D k−1 , p )

mod( D n−1 , p )

,...,x

(5)

) i −1

where k ∈ {1, 2, . . . , n}, f ( D , p ) comes from the substitutions xi = xmod( D , p ) , i = 1, 2, . . . , n, and i −1 k −1 f ( D , p ,k) comes from the substitutions xi = xmod( D , p ) , i = 1, 2, . . . , n, i = k, xk = x p +mod( D , p ) . Note that when n = 1, f ( D , p ) = f ( D , p ,k) = f (x). Substitution (4) was introduced in Klivans and Spielman (2001) and substitution (5) is introduced in this paper. We have

deg( f ( D , p ) ) < Dp and deg( f ( D , p ,k) ) < 2Dp .

(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 exists an aw ∈ M f \{cm} such that m( D , p ) = w ( D , p ) or m( D , p ,k) = w ( D , p ,k) . We now show how to construct S f ( D , p) and S f ( D , p,k) , which are SLPs to compute f ( D , p ) and f ( D , p ,k) , respectively. Lemma 5.1. Let S f be an SLP procedure to compute f , which has length L. Then we can design a new SLP S f ( D , p) (S f ( D , p,k) ) to compute f ( D , p ) ( f ( D , p ,k) ), which has length O ( L log p ). The design procedure costs O ( L log D + L log n log p + L log2 p ) bit operations. Proof. Assume S f consists of the operations

 i : b i ← α i  i β i , i = 1, 2, . . . , L with input {a1 , a2 , . . . , an }. First observe that for a monomial xk , we can construct an SLP to represent it with length O (log k) with bit complexity O ∼ (log2 k) via multiplication explicitly. Define an SLP S f ( D , p) for f ( D , p ) as follows. As f ( D , p ) comes from the substitutions xi = i −1

i −1

xmod( D , p ) , i = 1, 2, . . . , n, when we replace ai in S f by xmod( D , p ) , the new procedure is used to compute f ( D , p ) . This new procedure is not an SLP because it includes exponentiations. So we rei −1 i −1 place each xmod( D , p ) with another SLP that computes xmod( D , p ) to obtain an SLP. As the exponent i −1 mod( D , p ) < p, the new SLP has length at most O ( L log p ).

i −1 Now we analyze the complexity of the construction. For each xmod( D , p ) , it needs O ∼ (log D + i −1 i −1 log n log p ) bit operations to compute mod( D , p ). As mod( D , p ) < p, it needs at most O (log2 p )

bit operations to construct the SLP for xmod( D complexity is O ∼ ( L log D

i −1

, p ) . Since there are L instructions in

S f , the total bit

2

+ L log n log p + L log p ).

The construction of S f ( D , p,k) is similar to that for S f ( D , p) . The only difference is that when j = k, then replace ak by xmod( D

k−1

, p )+ p .

2

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.12 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

12

Let d1 d2 dr f (mod D , p ) = a1 x + a2 x + · · · + ar x (d1 < d2 < · · · < dr )

(7)

Since f ( D , p ) mod (x p − 1) = f ( D , p ,k) mod (x p − 1) = f (mod D , p ) , for k = 1, 2, . . . , n, we can write

f ( D , p) = f 1 + f 2 + · · · + f r + g

(8)

f ( D , p ,k) = f k,1 + f k,2 + · · · + f k,r + gk where f i mod (x p − 1) = f k,i mod (x p − 1) = ai xdi , g mod (x p − 1) = gk mod (x p − 1) = 0. Similar to (3), we define the following key notation e i ,1

f

e i ,n

M D , p = {ai x1 · · · xn |ai is from (7) for some i ∈ [1, r ], and M1 : f i = ai xu i , f k,i = ai xbk,i , k = 1, 2, . . . , n. bk ,i − u i ∈ N, k = 1, 2, . . . , n. M2 : e i ,k = p

(9)

M3 : u i = e i ,1 + e i ,2 mod( D , p ) + · · · + e i ,n mod( D n−1 , p ). M4 :

n 

e i , j < D .}

j =1

Lemma 5.2. Let f =

t

i =1 c i mi

f

∈ R[X], D > deg( f ). If c i mi is not a collision in f (mod D , p ) , then c i mi ∈ M D , p . f

Proof. It suffices to show that c i mi satisfies the conditions of the definition of M D , p . Assume e e e mod ds mi = x11 x22 · · · xnn . Since c i mi is not a collision in f (mod and d s = D , p ) , assume (c i mi )( D , p ) = a s x

n

j −1

, p ), where as xds is defined in (7). It is easy to see that c i mi is also not a collin i −1 s for u = sion in f ( D , p ) and in f ( D , p ,k) . Hence, f s = as xu , p ); bk,s = u s + pek . Clearly, s i =1 e i mod( D n M1, M2 and M3 are correct. Since deg(mi ) = j =1 e i , j < D, M4 is correct. 2 mod(

j =1 e j D

f

Now we give the following algorithm to compute M D , p , whose correctness is obvious. Algorithm 5.3 (MTerms). Input: Univariate polynomials f (mod D , p ) , f ( D , p ) , f ( D , p ,k) , where k = 1, 2, . . . , n, a prime p, D > deg( f ). f

Output: M D , p . Step 1: Write f (mod D , p ) , f ( D , p ) , and f ( D , p ,k) as forms (7) and (8) d1 d2 dr f (mod D , p ) = a 1 x + a 2 x + · · · + ar x

f ( D , p) = f 1 + f 2 + · · · + f r + g f ( D , p ,k) = f k,1 + f k,2 + · · · + f k,r + gk , (k = 1, 2, . . . , n). Step 2: Let S = {}. For i = 1, 2, . . . , r do a: If one of f i , f 1,i , . . . , f n,i is not of the following form: f i = ai xu i , f k,i = ai xbk,i , then break. b: Let e i ,k =

b k , i −u i p

for k = 1, 2, . . . , n. If e i ,k ∈ / N , then break.

n−1 c: If u , p ), then break; i = e i ,1 + e i ,2 mod( D , p ) + · · · + e i ,n mod( D n d: If e ≥ D, then break; i , j j =1

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.13 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

e: Let S = S



e

13

e

{ai x1i,1 · · · xni,n }.

Step 3 Return S. Lemma 5.4. Algorithm 5.3 needs O (nT ) arithmetic operations in R and O ∼ (nT log( p D )) bit operations. 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 three steps. First, we perform the modular operation mod p on every degree of f ( D , p ) , f ( D , p ,k) , which costs O ∼ (nT log( p D )) bit operations, since each of f ( D , p ) , f ( D , p ,k) has no more than T terms and the height of the degree is O (log( p D )). Second, we sort the terms of f (mod D , p ) , f ( D , p ) , f ( D , p ,k) into ascending order according to the new degree modulo p, which costs O ∼ (nT log p ) bit operations, since the degrees are < p. In order to check whether f i mod (x p − 1) = f k,i mod (x p − 1) = ai xdi , we need O (nT ) operations in R. Finally, f i , f k,i can be obtained with O (nT ) comparisons of the degrees, which costs O ∼ (nT log p ) bit operations and O (nT ) R-operations. So, the total complexity of Step 1 is O ∼ (nT log( p D )) bit operations and O (nT ) operations in R. 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( p D )), Steps b, c, and d cost O (n log( p D )) bit operations. Since we ∼ have at most T loops, the total complexity is O (nT log( p D )) bit operations. In a of Step 2, since # f ( D , p ) , # f ( D , p ,k) ≤ T , it totally needs O ∼ (nT ) bit operations to determine whether # f i , #( f k,i ) = 1. To compare the coefficients of f i , f k,i , it needs O (nT ) arithmetic operations in R. We proved the lemma. 2 5.2. The interpolation algorithm In this section, we give the interpolation algorithm for multivariate polynomial. We first give a sub-algorithm, which computes f ( D , p ) , and f ( D , p ,k) efficiently. Algorithm 5.5 (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, 2, . . . , n. e i ,1 e i ,2

Step 1: Assume f = c 1 m1 + c 2 m2 + · · · + ct mt , where mi = x1 x2 Step 2: Let u 1 = 1; For i = 2, 3, . . . , n do u i = mod(u i −1 D , p ). Step 3: Let h0 = 0. For i = 1, 2, . . . , n, let h i = 0; Step 4: For i = 1, 2, . . . , t do a: b: c: d:

e

· · · xni,n , i = 1, 2, . . . , t.

Let d = 0. For k = 1, 2, . . . , n, let d = d + e i ,k uk . h0 = h0 + c i xd ; For k = 1, 2, . . . , n, let hk := hk + c i xd+e i,k p .

Step 5: Return h0 , h i , i = 1, 2, . . . , n; Lemma 5.6. Algorithm 5.5 is correct. The complexity is O ∼ (nt log p + nt log D ) bit operations and O (nt ) arithmetic operations in R. Proof. In Step 2, u i = mod( D i −1 , p ). In b of Step 4, d is the degree of mi ( D , p ) , so h0 is f ( D , p ) after finishing Step 4. In d, since deg(mi ( D , p ,k) ) = deg(mi ( D , p ) ) + pe i ,k , hk is f ( D , p ,k) after finishing Step 4. So the correctness is proved. Now we analyze 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(mi ( D , p ) ) is O ( p · deg( f )) ≤ O ( p D ),

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.14 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

14

the bit operations is O ∼ (nt log( p D )). In c and d, it needs O ∼ (nt log( p D )) bit operations and O (nt ) arithmetic operations in R. 2 We now give the interpolation algorithm. Algorithm 5.7 (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 1 = max{1, n( T − 1) log2 D }, N 2 = max{1, nT log2 D }, N = max{4N 1 , N 1 + N 2 − 1}. Step 2: Find the first N different primes p 1 , p 2 , . . . , p N . mod = Step 3: For j = 1, 2, . . . , N, construct the SLP S f ( D , p ) with Lemma 5.1 and probes f (mod D , p ) . Let f j j

j

f (mod D,p ). j

Step 4: Let h = 0 and T 1 = T . Step 5: Loop 5.1: 5.2: 5.3: 5.4:

Let α = max{# f jmod | j = 1, 2, . . . , N } and j 0 the smallest number such that # f jmod = α . 0 If α = 0 return h. ∗ ∗ ∗ { f , f 1 , . . . , f n } = Substitution(h, p j 0 , D ). Let f j 0 = UIPoly(S f ( D , p ) , f ∗ , T , T 1 , Dp j 0 ). j0

5.5: For k = 1, 2, . . . , n, construct the SLP S f ( D , p UIPoly(S f ( D , p 5.6: Let

f −h M D,p j0

5.7: Let s = 0.

j 0 ,k)

,

f ∗, T , T k

j 0 ,k)

with Lemma 5.1 and let

gk =

1 , 2Dp j 0 ).

:= MTerms( f jmod , f j 0 , g 1 , g 2 , . . . , gn , 0 f −h For each u ∈ M D , p , if j0

p j 0 , D ).

mod #{ j | #( f jmod − u mod ), j = 1, . . . , N 1 + N 2 − 1} ≥ N 2 , ( D , p j ) ) < #( f j

then s := s + u. 5.8: Let h = h + s, T 1 = T 1 − #s, N 1 = max{1, n( T 1 − 1) log2 D }, N 2 = nT 1 log2 D , N = max{4N 1 , N 1 + N 2 − 1}. 5.9: For j = 1, 2, . . . , N, let f jmod = f jmod − smod (D,p ) . j

Theorem 5.8. Algorithm 5.7 finds f using O ∼ ( Ln2 T 2 log2 D + LnT log3 D ) operations in R and a similar number of bit operations. Proof. The algorithm is quite similar to the univariate interpolation Algorithm 4.6. So, we will only give the sketch of the proof and give detailed proof only for those steps which are essentially different from that in Algorithm 4.6. By Theorem 3.4 and Lemma 5.2, at least half of terms of f − h are in f −h

M D,p

j0

f −h M D,p . j0

obtained in Step 5.6. In Step 5.7, Theorem 2.4 is used to select the elements of M f −h from So at least half of the terms of f − h will be found in each loop of Step 5. Then the correctness

of the algorithm is proved. We now analyze 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 log2 N log log N ) by von zur Gathen and Gerhard (2013, p. 500, Thm. 18.10) and N is O ∼ (nT log D ), the bit complexity of Step 2 is O ∼ (nT log D ). In Step 3, we construct O (nT log D ) SLPs, since p j is O ∼ (nT log D ), by Lemma 5.1, it needs O ∼ ( LnT log2 D ) bit operations. We probe SLPs for O (nT log D ) times, so the cost of probes is O ∼ ( Ln2 T 2 log2 D ) operations in R and similarly many bit operations.

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.15 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

15

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 jmod ≤ T , the bit complexity is O ∼ (nT 2 log D ). 0 In Steps 5.3, by Lemma 5.6, since p i is O ∼ (nT log D ), the complexity is O ∼ (n log D + nT log p j 0 + nT log D ) = O ∼ (nT log D ) bit operations and O (nT ) arithmetic operations in R. In Steps 5.4 and 5.5, since O ( p j 0 ) is O ∼ (nT log D ), by Lemma 5.1, it needs O ∼ (nL log D + nL log( T log D ) + nL log2 ( T log D )) bit operations to construct the SLPs S f ( D , p ,k) . As the length of

S f (D,p j

0

)

and S f ( D , p

j0

is O ∼ ( L log(nT log D )) and the numbers of terms and degrees of ( f −

j 0 ,k)

h)( D , p j ) , ( f − h)( D , p j ,k) are bounded by O ( T ) and O ∼ (nT D ), by Theorem 4.7, the complexity is 0 0 O ∼ ( LnT 2 log2 D + LnT log3 D ) operations in R and similarly many bit operations. In Step 5.6, by Theorem 5.4, the complexity is O ∼ (nT log D ) bit operations and O (nT ) ring operations in R. f −h ∼ In Step 5.7, to all the u mod ( D , p ) , we need O (n( N 1 + N 2 )#M D , p log( Dp )) bit operations. The proof j0

j

for rest of this step is similar to that of Step 5.5 of Algorithm 4.6. The complexity is O ∼ (n( N 1 + f −h

N 2 )#M D , p

j0

f −h

log( Dp j ) + #M D , p

j0

f −h

log T ( N 1 + N 2 ) log( Dp j )) bit operations and O (#M D , p ( N 1 + N 2 )) j0

f −h

ring operations. Since #M D , p ≤ T , p j = O ∼ (nT log D ), and N 1 + N 2 = O (nT log D ), it needs j0

O ∼ (n2 T 2 log2 D ) bit operations and O ∼ (nT 2 log D ) ring operations in R. mod In Step 5.9, we need n#s operations in Z to obtain smod ( D , p j ) . Subtract s( 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 in R. Since the sum of #s is t, it total costs O ∼ (nT 2 log D ) ring operations in R and O ∼ (n2 T 2 log2 D ) bit operations. The total complexity of one loop of Step 5 is O ∼ ( LnT 2 log2 D + LnT log3 D ) ring operations in R and O ∼ ( LnT 2 log2 D + LnT log3 D + n2 T 2 log2 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 2 log2 D + LnT log3 D ) ring operations in R and O ∼ ( LnT 2 log2 D + LnT log3 D + n2 T 2 log2 D ) bit operations. Plus the complexity of Step 3, the complexity of the algorithm is O ∼ ( Ln2 T 2 log2 D + LnT log3 D ) operations in R and similarly many bit operations, which is from Step 3 and Step 5.5. 2 6. Monte Carlo algorithm for multivariate polynomials In this section, we give Monte Carlo interpolation algorithms for multivariate polynomials, which could be considered as probabilistic versions of Algorithm 5.7. The following theorem shows how to use a probabilistic method to obtain a p such that the number of collision terms of f in f (mod D , p ) is very small, which is a probabilistic version of Theorem 3.4.

t

∈ R[X], T ≥ # f , D > deg( f ), N = max{1,  32 n( T − 1) log2 D }, and 5 p 1 , p 2 , . . . , p 2N be 2N different primes, 0 < ε < 1, k = log2 ( 1ε ). If j 1 , j 2 , . . . , jk are randomly chosen from mod mod mod [1, 2N ] and j 0 ∈ { j 1 , . . . , jk } is the integer such that # f (mod D , p j ) = max{# f ( D , p j ) , # f ( D , p j ) , . . . , # f ( D , p j ) }. Theorem 6.1. Let f =

i =1 c i mi

0

f

1

2

k

5 Then with probability ≥ 1 − ε , C( D , p ) ≤  16 T . j0

Proof. If T = 1, f has only one term, so there cannot be any collisions in f for any prime, the theorem is correct. So we assume T ≥ 2, then N =  32 n( T − 1) log2 D . First we claim that if randomly 5 f 1 5 , C( D , p ) < 32 T . It suffices to show 2 j f 5 that C( D , p ) < 32 T . We prove this by contraj f 5 integers in [1, 2N ] such that C( D , p ) ≥ 32 T. αi

choose an integer j in [1, 2N ], then with probability at least that there exist at least N integers j in [1, 2N ] such diction. Assume

α1 , α2 , . . . , αN +1 are N + 1 different f

C( D , p α ) /2

By Lemma 3.3, p αi

i

divides A, where A is defined in Lemma 3.3. Since p αi are differ-

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.16 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

16 f

ent and C( D , p ) ≥ αi

5 T, 32

5

5

( p α1 p α2 · · · p αN +1 ) 64 T  divides A. Now we have ( p α1 p α2 · · · p αN +1 ) 64 T  ≥

5 32 5 1 1 2( N +1) 64 T  ≥ 2 5 n( T −1) log2 D  64 T  ≥ 2 2 nT ( T −1) log2 D = D 2 nT ( T −1) , which contradicts to A ≤ ( D n − 1

1) 2 t (t −1) . So in { p 1 , p 2 , · · · , p 2N }, there are at least N primes p i such that C( D , p ) < i proved the claim. f

f

If there exists at least one j i in { j 1 , j 2 , . . . , jk } such that C( D , p ) < ji

f C( D , p j ) 0

≤

f 2C( D , p ) ji

<

5 T. 16

So

f C( D , p j ) 0

≥

5 T 16

only when all

just proved, the probability for this to happen is at most f

f C( D , p j ) i

( 12 )k

≤

≥

5 T, 32

5 T,i 32

1 ( 12 )log2 ε

5 T. 32

We have

then by Corollary 3.2,

= 1, 2, . . . , k, by the claim f

= ε . Since C( D , p j

f

0

)

<

5 T 16

5 5 T , the probability of C( D , p ) ≤  16 T is at least 1 − ε . The theorem is implies that C( D , p ) ≤  16 j0 j0 proved. 2

Remark 6.2. In Theorem 6.1, each j i is randomly chosen from [1, 2N ] independently, so j 1 , . . . , jk may contain repeated numbers. If ε is small enough such that k > 2N, j 1 , . . . , jk must contain repeated numbers. The theorem is still valid in these cases, since each choice is an independent repeated trial. f

5 Note that the result C( D , p ) ≤  16 T of Theorem 6.1 is different with that of Theorem 3.4. To j0

find a p such that at least  2t  of the terms of f are not collisions in f (mod D , p j 0 ) is not enough for our probabilistic algorithm. Lemma 6.3. If ε ∈ (0, 1) and a ≥ 1, then (1 − ε )a ≥ 1 − aε . Proof. By Taylor expansion, we have (1 + x)a = 1 + ax + a(a − 1) a−2

let x = −ε , then (1 − ε )a = 1 − aε + a(a − 1) (1−θ2ε) a(a − 1)

(1−θ ε )a−2 2

(1+θ x)a−2 2

. Since θ ∈ (0, 1),

, where θ ∈ (0, 1). Now we

ε ∈ (0, 1) and a ≥ 1, we have

≥ 0. So we have (1 − ε )a ≥ 1 − aε . 2

We first consider interpolation over an arbitrary computable ring. For the univariate interpolation algorithm, we use the following algorithm given in Arnold et al. (2013), which is the fastest known probabilistic algorithm over arbitrary rings. Theorem 6.4. (Arnold et al., 2013) 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 log3 D + LT log D log ν1 ) 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 and on any invocation returns the correct answer with probability greater than 1 − ν , for a user-supplied tolerance 0 < ν < 1. We use PUniPoly1 (S f , f ∗ , T 1 , D , ν ) to denote the algorithm in Theorem 6.4, where f ∗ is a current approximation to f and #( f − f ∗ ) ≤ T 1 , # 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 1 ≥ #( f − f ∗ ), D > max{deg( f ), deg( f ∗ )}, a tolerance ν such that 0 < ν < 1. Output: With probability ≥ 1 − ν , return a polynomial f ∗∗ such that #( f − f ∗ − f ∗∗ ) ≤  T21 . ν , and k = log 1 . Find the first 2N primes Step 1: Let N = max{1,  32 n( T 1 − 1) log2 D }, ε = n+ 2 ε 5 2 p 1 , p 2 , . . . , p 2N . Step 2: Let j 1 , j 2 , . . . , jk be randomly chosen from [1, 2N ]. mod ∗mod Step 3: For i = 1, 2, . . . , k, construct SLP S f ( D , p ) and probe f (mod D , p ) . Let g i = f ( D , p ) − f ( D , p ) .

Step 4: Let

ji

ji

ji

ji

α = max{#g i |i = 1, 2, . . . , k} and j 0 satisfying #g j0 = α . If α ≥ T , then return failure.

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.17 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

17

Step 5: Let {q, q1 , . . . , qn } = Substitution( f ∗ , p j 0 , D ). Step 6: Construct the SLPs S f ( D , p ) and S f ( D , p ,i) , i = 1, . . . , n with Lemma 5.1. Step 7: For i = 1, . . . , n, let

j0

j0

η = PUniPoly1 (S f ( D , p j ) , q, T 1 , p j0 D , ε), ηi = PUniPoly1 (S f ( D , p j 0

2p j 0 D , ε ). If η or ηi is failure, then return failure. Step 8: Let M = MTerms  ( g j 0 , η, η1 , η2 , . . . , ηn , D , p j 0 ). Step 9: Return f ∗∗ = s∈ M s.

0

,i )

, qi , T 1 ,

Lemma 6.6. Algorithm 6.5 computes f ∗∗ such that #( f − f ∗ − f ∗∗ ) ≤  T21 with probability at least 1 − ν . The algorithm costs O ∼ ( LnT log3 D + LnT log D log ν1 ) ring operations in R and similarly many bit operations.

Proof. We first show that Algorithm 6.5 returns the polynomial f ∗∗ such that #( f − f ∗ − f ∗∗ ) ≤ f−f∗

 T21 with probability 1 − ν . In Step 4, by Theorem 6.1, with probability 1 − ε , C( D , p j f−f∗

0

)

5 ≤  16 T 1 . If

5 j 0 satisfies C( D , p ) ≤  16 T 1 and η = ( f − f ∗ )( D , p j ) , ηi = ( f − f ∗ )( D , p j ,i ) , i = 1, 2, . . . , n, then by j0 0 0 5 Lemma 5.2, f − f ∗ contains at most  16 T 1 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 of f ∗∗ not in f f ∗ − f ∗∗ ) ≤  1 T

f ∗ . So #( f

f∗

f ∗∗ )

5  16 T 1

1 5  T 2 16 1

15 T 32 1

1 T 2 1

1 5  T 2 16 1

terms

− − − ≤ + ≤ < and we have #( f − . In Step 7, by Theorem 6.4, the probability of obtaining the correct ( f − f ∗ )( D , p j0 ) 1 2 ∗ ∗ ∗ (( f − f )( D , p j0 ,i ) ), or η = ( f − f )( D , p j0 ) (ηi = ( f − f )( D , p j0 ,i ) ), is ≥ 1 − ε , so the probability of obtaining correct polynomials ( f − f ∗ )( D , p j ) and (( f − f ∗ )( D , p j ,i ) ), i = 1, . . . , n is ≥ (1 − ε )n+2 . By 0 0 ν , (1 − ε )n+2 ≥ 1 − ν . Hence, with probability Lemma 6.3, (1 − ε )n+2 ≥ 1 − (n + 2)ε . Since ε = n+ 2 ≥ 1 − ν , we obtain an f ∗∗ satisfying #( f − f ∗ − f ∗∗ ) ≤  T21 . The correctness of the lemma is proved. Now we analyze the complexity. Since the bit complexity of finding the first 2N primes is O ( N log2 N log log N ) by von zur Gathen and Gerhard (2013, p. 500, Thm. 18.10) and N is O ∼ (nT log D ), the bit complexity of Step 1 is O ∼ (nT log D ). In Step 2, since probabilistic machines flip coins to decide binary digits, each of these random choices can be simulated with a machine with complexity O (log 2N ). So the complexity of Step 2 is O (log 1ε log(nT log D )) bit operations. Since O (log ε1 ) is O (log n + log ν1 ), the bit complexity of Step 2 is O (log n log(nT log D ) + log(nT log D ) log ν1 ).

In Step 3, we construct k = log2 ε1  SLPs, since p j i is O ∼ (nT log D ), by Lemma 5.1, it needs + L log n log(nT log D ) + L log2 (nT log D )) bit operations. We probe f for k times and the cost of the probes is O ∼ ( LnT log D log 1ε ) ring operations in R and a similar many bit operations. So O ∼ ( L log D

the complexity of step 3 is O ∼ ( LnT log D log ν1 ) arithmetic operations in R and a similar number of bit operations. 1 ∼ In Step 4, we find the integer j 0 . Since #( f − f ∗ )mod ( D , p ) ≤ T , it needs at most O ( T log ε ) bit ji

∼ operations to compute all #( f − f ∗ )mod ( D , p j i ) , i = 1, 2, . . . , k. Find j 0 needs O ( T ) bit operations. So the

bit complexity of Step 4 is O ∼ ( T log n + T log ν1 ). In Step 5, by Lemma 5.6, it needs O ∼ (nT log D ) bit operations and O (nT ) ring operations. In Step 6, since p j 0 is O ∼ (nT log D ), by Lemma 5.1, it needs O ∼ ( Ln log D + Ln log( T log D ) + Ln log2 ( T log D )) bit operations. In Step 7, we call n + 1 times Algorithm PUniPoly1 . Since the lengths of SLPs S f ( D , p

j0 )

and S f ( D , p

j 0 ,i )

are O ( L log( p j 0 )) and the terms and degrees of ( f − f ∗ )( D , p j ) , ( f − f ∗ )( D , p j ,k) are respectively 0 i

bounded by T and 2p j 0 D, by Theorem 6.4, the complexity of Step 7 is O ∼ ( L log( p j 0 )nT log3 ( p j 0 D ) + L log( p j 0 )nT log( p j 0 D ) log ε1 )) arithmetic operations in R and plus a similar number of bit operaν , p ∼ ∼ tions. Since ε = n+ j 0 is O (nT log D ), and 2p j 0 D is O (nT D ), the complexity of step 6 is 2 O ∼ ( LnT log3 D + LnT log D log ν1 )) ring operations in R and similarly many bit operations. In Step 8, by Theorem 5.4, the complexity is O (nT ) arithmetic operation in ring R and O ∼ (nT log p D ) bit operations.

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.18 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

18

It is easy to see that the complexity is dominated by Step 3 and Step 6. The theorem is proved.

2

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 Output: With probability at least 1 − μ, return f .

μ ∈ (0, 1).

μ

Step 1: Let h = 0, T 1 = T , ν = log T +1 . 2 Step 2: While T 1 > 0 do b: Let g = HalfPoly(S f , h, T , T 1 , D , ν ). If g = f ailure, then return failure. c: Let h = h + g, T 1 =  T21 . Step 3: Return h. Theorem 6.8. Algorithm 6.7 computes f with probability ≥ 1 − μ. The algorithm costs O ∼ ( LnT log3 D + 1 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 ) ≤  T21 with probability ≥ 1 − ν . By Lemma 6.3, Step 2

will run at most k = log2 T  + 1 times and return the correct f with probability ≥ (1 − ν )k ≥ 1 − μ. The correctness of the theorem is proved. Now we analyze 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 Lemma 6.6, the complexity of Step 2 is O ∼ ( LnT log3 D + LnT log D log ν1 ) ring operations in R and a similar number ∼ log2 T +1 , the complexity of Step 2 is O ( LnT log D

μ

3

+ LnT log D log μ1 ) ring operations in R and a similar number of bit operations. The theorem is proved. 2 of bit operations. Since

Remark 6.9. Set

ν=

μ = 1/4. Then Algorithm 6.7 computes f with probability at least 34 . The algorithm

costs O ∼ ( LnT log3 D ) ring operations in

R and a similar number of bit operations.

We now consider interpolation over finite fields. In Arnold et al. (2014), Arnold, Giesbrecht & Roche gave a new univariate interpolation algorithm for the finite field with better complexities. Theorem 6.10. (Arnold et al., 2014) Let f ∈ Fq [x] with at most T non-zero terms and degree at most D, and let 0 < η ≤ 12 . Suppose we are given an SLP S f of length L that computes f . Then there exists an algorithm that interpolates f , with probability at least 1 − η , with a cost of O ∼ ( LT log2 D (log D + log q) log η1 ) bit operations.

We use PUniPoly2 ( f , f ∗ , T 1 , D , η) to denote the algorithm in Theorem 6.10, where f ∗ is a current approximation to f and #( f − f ∗ ) ≤ T 1 , # f ∗ ≤ T , deg( f ∗ ) < D , deg( f ) < D. Replacing Algorithm PUniPoly1 with Algorithm PUniPoly2 in Step 7 of Algorithm HalfPoly, we obtain a multivariate interpolation algorithm for finite fields. We assume S f can evaluate in an extension field of Fq and have the following result. Theorem 6.11. Let f ∈ Fq [X] be given as an SLP, with at most T non-zero terms and degree at most D, and let 0 ≤ μ < 1/2. Then we can interpolate f , with probability at least 1 − μ, with a cost of O ∼ ( LnT log2 D (log D + 1 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.10 instead of Theorem 6.4. 2

JID:YJSCO

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.19 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

19

7. Maple code We implement all the algorithms proposed in this paper in Maple. The Maple codes can be found in http://www.mmrc.iss.ac.cn/~xgao/software/slppoly.zip The codes also include some test examples and their running times. The data are collected on a desktop with Windows system, 3.60 GHz Core i7-4790 CPU, and 8 GB RAM memory. We randomly construct five polynomials, then regard them as SLP polynomials and reconstruct them with Algorithms 4.6, 5.7 and 6.7. We just collect the time of recovering f from the univariate polynomials f (mod D , p ) , f ( D , p ) , and f ( D , p ,k) , which means the running time does not include the (dominant) cost of probing the SLP. The average times for the five examples are collected, which are in accordance with the complexity bounds given in this paper. For instance, the experimental results show that the runtime for Algorithm 5.7 is roughly quadratic in n and T and cubic in log D, respectively. 8. Conclusion In this paper, we give a new deterministic interpolation algorithm and two Monte Carlo interpolation algorithms for SLP sparse multivariate polynomials. Our deterministic algorithm has better complexity than existing deterministic interpolation algorithms in most cases. Our Monte Carlo interpolation algorithms are the first algorithms whose complexities are linear in nT and polynomial in log D. The algorithms are based on several new ideas. In order to have a deterministic algorithm, we give a criterion for checking whether a term belongs to a polynomial. We also give a deterministic method to find an “ok” prime p in the sense that at least half of the terms in f are not collisions in f (mod D , p ) . Finally, a new Kronecker type substitution is given to reduce multivariate polynomial interpolations to univariate polynomial interpolations. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement We thank the anonymous referee for suggesting a better bound in Lemma 4.4. References Alon, N., Mansour, Y., 1995. Epsilon-discrepancy sets and their application for interpolation of sparse polynomials. Inf. Process. Lett. 54 (6), 337–342. Arnold, A., 2016. Sparse Polynomial Interpolation and Testing. PhD thesis. University of Waterloo. http://hdl.handle.net/10012/ 10307. Arnold, A., Roche, D.S., 2014. Multivariate sparse interpolation using randomized Kronecker substitutions. In: International Symposium on Symbolic & Algebraic Computation, pp. 35–42. Arnold, A., Giesbrecht, M., Roche, D.S., 2013. Faster sparse interpolation of straight-line programs. In: Computer Algebra in Scientific Computing - 15th International Workshop, CASC 2013, Proceedings. Berlin, Germany, September 9–13, 2013, pp. 61–74. Arnold, A., Giesbrecht, M., Roche, D.S., 2014. Sparse interpolation over finite fields via low-order roots of unity. In: International Symposium on Symbolic and Algebraic Computation, ISSAC ’14. Kobe, Japan, July 23-25, 2014, pp. 27–34. Arnold, A., Giesbrecht, M., Roche, D.S., 2016. Faster sparse multivariate polynomial interpolation of straight-line programs. J. Symb. Comput. 75, 4–24. https://doi.org/10.1016/j.jsc.2015.11.005. Avendaño, M., Krick, T., Pacetti, A., 2006. Newton-Hensel interpolation lifting. Foundations of Computational Mathematics 6 (1), 82–120. https://doi.org/10.1007/s10208-005-0172-3. Ben-Or, M., Tiwari, P., 1988. A deterministic algorithm for sparse multivariate polynomial interpolation. In: Proceedings of the Twentieth Annual ACM Symposium on Theory of Computing. ACM, pp. 301–309. Garg, S., Schost, É., 2009. Interpolation of polynomials given by straight-line programs. Theor. Comput. Sci. 410 (27–29), 2659–2662. https://doi.org/10.1016/j.tcs.2009.03.030. Giesbrecht, M., Roche, D.S., 2011. Diversification improves interpolation. In: Symbolic and Algebraic Computation, International Symposium, ISSAC 2011 (co-located with FCRC 2011), Proceedings. San Jose, CA, USA, June 7-11, 2011, pp. 7–11. Huang, Q.-L., Gao, X.-S., 2019. Deterministic sparse interpolation of black-box multivariate polynomials using Kronecker type substitutions (in Chinese). Accepted by Sci. Sin. Math. https://doi.org/10.1360/N012019-00018.

JID:YJSCO

20

AID:1968 /FLA

[m1G; v1.261; Prn:24/10/2019; 10:23] P.20 (1-20)

Q.-L. Huang, X.-S. Gao / Journal of Symbolic Computation ••• (••••) •••–•••

Kaltofen, E., 1985. Computing with polynomials given by straight-line programs I: greatest common divisors. In: Proceedings of the 17th Annual ACM Symposium on Theory of Computing. Providence, Rhode Island, USA, May 6–8, 1985, pp. 131–142. Kaltofen, E., Lakshman, Y.N., 1988. Improved sparse multivariate polynomial interpolation algorithms. In: Symbolic and Algebraic Computation, International Symposium ISSAC’88, Proceedings. Rome, Italy, July 4–8, 1988, pp. 467–474. Klivans, A.R., Spielman, D.A., 2001. Randomness efficient identity testing of multivariate polynomials. In: Proceedings on 33rd Annual ACM Symposium on Theory of Computing. Heraklion, Crete, Greece, July 6–8, 2001, pp. 216–223. Kronecker, L., 1882. Grundzüge einer arithmetischen theorie der algebraischen grössen (abdruck einer festschrift zu herrn ee kummers doctor-jubiläum, 10. September 1881), J. Reine Angew. Math. 92, 1–122. Mansour, Y., 1995. Randomized interpolation and approximation of sparse polynomials. SIAM J. Comput. 24 (2), 357–368. https://doi.org/10.1137/S0097539792239291. von zur Gathen, J., Gerhard, J., 2013. Modern Computer Algebra, 3 ed.. Cambridge University Press. Zippel, R., 1990. Interpolating polynomials from their values. J. Symb. Comput. 9 (3), 375–403. https://doi.org/10.1016/S07477171(08)80018-1.