A fast algorithm for solving linearly recurrent sequences
aa r X i v : . [ c s . S C ] J un A fast algorithm for solving linearly recurrent sequences
Seung Gyu Hyun, Stephen Melczer, Catherine St-Pierre
Abstract
We present an algorithm which computes the D th term of a sequence satisfying a linearrecurrence relation of order d over a field K in O ( M ( ¯ d ) log( D ) + M ( d ) log( d )) operations in K ,where ¯ d ≤ d is the degree of the squarefree part of the annihilating polynomial of the recurrenceand M is the cost of polynomial multiplication in K . This is a refinement of the previouslyoptimal result of O ( M ( d ) log( D )) operations, due to Fiduccia. Overview.
Consider a sequence ( a i ) i ≥ with entries in a field K that is generated by the recurrence a i + d = d − X j =0 c j a i + j (1)for all i ≥
0, where c = 0. Given initial conditions a j , for 0 ≤ j ≤ d −
1, along with the annihilatingpolynomial P = x d − P d − j =0 c j x j , we are interested in the complexity of computing one term a D ofthe sequence for some index D ≫ a , . . . , a D , but one can do much better. It hasbeen known since at least Fiduccia’s work [2] that computing a D can be reduced to multiplicationmodulo P . Explicitly, define A = K [ x ] /P , together with the K -linear form ℓ : A → K given by ℓ ( x i ) = a i for i = 0 , . . . , d −
1. Then, since the residue class of x in A is a root of P , its powers ( x i ) in A satisfy (1), and so do the values ( ℓ ( x i )) i ≥ ; this implies that ℓ ( x i ) = a i holds for all values of i ≥ a D it is enough to compute R = x D mod P as R = r + · · · + r d − x d − ,since then a D = r a + · · · + r d − a d − .Letting M denote a function such that polynomials of degree n can be multiplied in M ( n )operations (under the assumptions of [6, Chapter 8]), computing R costs O ( M ( d ) log( D )) operationsin K ; we can then deduce a D in O ( d ) steps. A new algorithm.
In this note, we present an improvement over this previous result that isuseful when P has multiple factors of high multiplicities.First, we reduce to the case where P has the form P = Q m , for squarefree Q . To that end,let P = Q i Q m i i be the squarefree factorization of P , with pairwise distinct m , . . . , m s and Q i squarefree of degree f i for all i ; note d = P i f i m i . We can compute x D mod P by applying theChinese Remainder Theorem to the quantities x D mod Q m i i , giving the following algorithm.1 lgorithm 1 Computing one element in a linear recurrent sequence
Input: • P : characteristic polynomial of the sequence • v = [ a , . . . , a d − ]: vector of initial conditions • D : index Output: D th element of the sequence ( a i ) as in (1)1. compute the squarefree factorization of P as P = Q i Q m i i
2. for i = 1 , . . . , n , compute C i = x D mod Q m i i
3. compute R = x D mod P by CRT as R = r + · · · + r d − x d −
4. return a D = r a + · · · + r d − a d − .To compute the C i ’s efficiently, we will use bivariate computations. Indeed, for i = 1 , . . . , s define A i = K [ X ] /Q m i i ; then there exists a K -algebra isomorphism π i : A i = K [ X ] /Q m i i → K [ y, x ] / h Q i ( y ) , ( x − y ) m i i . Following van der Hoeven and Lecerf [5], we will call π i the operation of untangling (this is aconversion from a univariate representation to a bivariate one) and its inverse tangling . (To beprecise, van der Hoeven and Lecerf consider a mapping to K [ y, x ] / h Q i ( y ) , x m i i , which is isomorphicto K [ y, x ] / h Q i ( y ) , ( x − y ) m i i through the shift x x + y ).van der Hoeven and Lecerf prove that for a given index i , untangling can be done in O ( M ( f i m i ) log( m i ))operations in K ; note that input and output sizes are f i m i in this case. They also give an algorithmfor tangling of cost O ( M ( f i m i ) log ( m i ) + M ( f i ) log( f i )); we give below a Las Vegas algorithm ofcost O ( M ( f i m i ) log( f i m i )) for this task. Taking the existence of such algorithms for granted, wewrite C i = x D mod Q m i i = π − i ( δ i ) , with δ i = x D mod h Q i ( y ) , ( x − y ) m i i . The following allows us to compute δ i efficiently. Define coefficients e , . . . , e m i − by x D mod ( x − m i = e + e x + · · · + e m i − x m i − , and define S i ( x ) = ( yx ) D mod ( x − m i ∈ K [ y ][ x ], where y is seen as an element of K [ y ] /Q i ( y ).Then δ i = x D mod ( x − y ) m i = S i (cid:18) xy (cid:19) , where we note y is invertible in K [ y ] /Q i ( y ) as c = 0. Now, S i = y D x D mod ( x − m i = y D ( e + e x + · · · + e m i − x m i − ) , so that δ i = y D e + y D − e x + · · · + y D − ( m i − e m i − x m i − . In this algorithm, we first need to compute coefficients e , . . . , e m i − ; assuming that 2 , . . . , m i − K , they can be obtained in O (log( D ) + M ( m i )) operations in K . The powers of y we2eed are computed modulo Q i , in time O ( M ( f i ) log( D ) + m i M ( f i )). Altogether, we obtain δ i using O ( M ( f i ) log( D ) + M ( f i m i ))) operations in K . As said above, we can deduce C i from δ i in Las Vegastime O ( M ( f i m i ) log( f i m i )) . Taking all i ’s into account and using the super-linearity of M , the totaltime to compute C , . . . , C s is thus O ( M ( ¯ d ) log( D ) + M ( d ) log( d )) , where ¯ d = P i f i ≤ d is the degree of the squarefree part of P .Computing the squarefree factorization of P and Chinese remaindering both cost O ( M ( d ) log d )[6, Corollary 10.23], so the overall runtime is O ( M ( ¯ d ) log( D ) + M ( d ) log( d )). This is to be comparedwith the cost O ( M ( d ) log( D )) of Fiduccia’s algorithm. Tangling and untangling.
We conclude by sketching our new algorithm for tangling. van derHoeven and Lecerf reduce the tangling operation to untangling by means of a divide-and-conquerprocess; we propose a direct reduction that uses transposition, inspired by an algorithm from [4]that applies in univariate situations.In what follows we use the same notation as in the previous paragraphs, but we drop thesubscript i for clarity. In particular, we write f = deg( Q ) and n = f m for the degree of Q m ,that is, the input and output size. Given δ in K [ y, x ] / h Q ( y ) , ( x − y ) m i , we want to find C = c + · · · + c n − x n − such that π ( C ) = δ ; this simply means that C mod h Q ( y ) , ( x − y ) m i = δ. Choose a random linear form λ : K [ y, x ] / h Q ( y ) , ( x − y ) m i → K . For j ≥
0, multiply the formerequality by x j and apply λ ; this gives c λ ( x j ) + · · · + c n − λ ( x j + n − ) = λ ( x j δ ) . Taking j = 0 , . . . , n −
1, we can collect these equalities in a linear system HA = L , with H = λ (1) . . . λ ( x n − )... . . . ... λ ( x n − ) . . . λ ( x n − ) A = c ... c n − L = λ ( δ )... λ ( x n − δ ) . Once H and L are known, we can recover coefficients A in O ( M ( n ) log( n )) operations in K , sincethe linear system is Hankel (for a generic choice of λ , matrix H has full rank n ). Hence, the mainquestion is the efficient computation of the entries of matrices H and L . Both are instances of thesame problem: given a K -linear form λ over K [ y, x ] / h Q ( y ) , ( x − y ) m i , and β in K [ y, x ] / h Q ( y ) , ( x − y ) m i , compute the values λ ( β ) , . . . , λ ( x n − β ).As already recognized by Shoup in the univariate case, this question is the transpose of theuntangling map π − . As a result, using the so-called transposition principle [4, 3, 1], we can deducean algorithm of cost O ( M ( n ) log( n )) = O ( M ( f m ) log( f m )) for tangling, by transposition of van derHoeven and Lecerf’s untangling algorithm. References [1] A. Bostan, G. Lecerf, and ´E. Schost. Tellegen’s principle into practice. In
ISSAC’03 , pages37–44. ACM, 2003. 32] C. M. Fiduccia. An efficient formula for linear recurrences.
SIAM Journal on Computing ,14(1):106–112, 1985.[3] E. Kaltofen. Challenges of symbolic computation: my favorite open problems.
J. Symb. Comp. ,29(6):891–919, 2000.[4] V. Shoup. Fast construction of irreducible polynomials over finite fields.
Journal of SymbolicComputation , 17(5):371–391, 1994.[5] J. van der Hoeven and G. Lecerf. Composition modulo powers of polynomials. In
ISSAC ’17 ,pages 445–452. ACM, 2017.[6] J. von zur Gathen and J. Gerhard.