Computation of the Similarity Class of the p-Curvature
aa r X i v : . [ c s . S C ] M a y Computation of the Similarity Class of the p-Curvature ∗ Alin Bostan
Inria (France) [email protected]
Xavier Caruso
Université Rennes 1 (France) [email protected]
Éric Schost
Univ. of Waterloo (Canada) [email protected]
ABSTRACT
The p -curvature of a system of linear differential equationsin positive characteristic p is a matrix that measures howfar the system is from having a basis of polynomial solu-tions. We show that the similarity class of the p -curvaturecan be determined without computing the p -curvature itself.More precisely, we design an algorithm that computes the in-variant factors of the p -curvature in time quasi-linear in √ p .This is much less than the size of the p -curvature, which isgenerally linear in p . The new algorithm allows to answera question originating from the study of the Ising model instatistical physics. CCS Concepts • Computing methodologies → Algebraic algorithms;
Keywords differential equations; p -curvature; algebraic complexity
1. INTRODUCTION
Differential equations in positive characteristic p are im-portant and well-studied objects in mathematics [22, 32, 33].The main reason is arguably one of Grothendieck’s (still un-solved) conjectures [26,27,1], stating that a linear differentialequation with coefficients in Q ( x ) admits a basis of algebraicsolutions if and only if its reductions modulo (almost) allprimes p admit a basis of polynomial solutions modulo p .Another motivation stems from the fact that the reductionsmodulo prime numbers yield useful information about thefactorization of differential operators in characteristic zero.To a linear differential equation in fixed characteristic p ,or more generally to a system of such equations, is attacheda simple yet very useful object, the p -curvature . Let F q be the finite field with q = p a elements. The p -curvatureof a system of linear differential equations with coefficients ∗ We warmly thank the referees for their very helpful comments,and for pointing out a mistake in an earlier version of the article.We thank M. Giesbrecht, G. Labahn and A. Storjohann for usefuldiscussions. The third author is supported by NSERC.
ISSAC ’16, July 19–22, 2016, Waterloo, ON, Canada.
ACM ISBN 978-1-4503-2138-9.DOI: http://dx.doi.org/10.1145/2930889.2930897 in F q ( x ) is a matrix with entries in F q ( x ) that measures theobstructions for the given system to possess a fundamen-tal matrix of polynomial solutions in F q [ x ]. Its definitionis remarkably simple, especially at a higher level of gen-erality: the p -curvature of a differential module ( M, ∂ ) ofdimension r over F q ( x ) is the “differential-Frobenius-map” ∂ p = ∂ ◦ · · · ◦ ∂ ( p times). When applied to the differentialmodule canonically attached with the system Y ′ = A ( x ) Y ,the p -curvature materializes into the p -th iterate ∂ pA of themap ∂ A : F q ( x ) r → F q ( x ) r that sends v to v ′ − Av , or moreconcretely, into the matrix A p ( x ) of this map with respectto the canonical basis of F q ( x ) r . It is given as the term A p of the sequence ( A i ) i of matrices in M r ( F q ( x )) defined by A = − A and A i +1 = A ′ i − A · A i for i ≥ . From a computer algebra perspective, many effectivityquestions naturally arise. They primarily concern the al-gorithmic complexity of various operations and propertiesrelated to the p -curvature: How fast can one compute A p ?How fast can one decide its nullity? How fast can one de-termine its minimal and characteristic polynomial? Apartthe fundamental nature of these questions from the algebraiccomplexity theory viewpoint, there are concrete motivationsfor the efficient computation of the p -curvature, coming fromvarious applications, notably in enumerative combinatoricsand statistical physics [7, 8, 2].We pursue the algorithmic study of the p -curvature, ini-tiated in [9, 3, 4]. In those articles, several questions wereanswered satisfactorily, but a few other problems were leftopen. In summary, the current state of affairs is as follows.First, the p -curvature A p can be computed in time O (log p )when r = 1 and O ˜( p ) when r >
1. The soft-O notation O ˜( ) indicates that polylogarithmic factors in the argumentof O ( ) are deliberately not displayed. These complexitiesmatch, up to polylogarithmic factors, the generic size of A p .Secondly, one can decide the nullity of A p in time O ˜( p ) andcompute its characteristic polynomial in time O ˜( √ p ). Itis not known whether the exponent 1 / ± , × , ÷ ) inthe ground field F q , and the dependence is expressed in themain parameter p only. Nevertheless, precise estimates arealso available in terms of the other parameters of the input.In the present work, we focus on the computation of allthe invariant factors of the p -curvature, and show that theycan also be determined in time O ˜( √ p ). Previously, thiswas unknown even for the minimal polynomial of A p or fortesting the nullity of A p . The fact that a sublinear costcould in principle be achievable, although A p itself has aotal arithmetic size linear in p , comes from the observationthat the coefficients of the invariant factors of A p lie in thesubfield F q ( x p ) of F q ( x ), in other words they are very sparse.To achieve our objective, we blend the methods used inour previous works [3] and [4]. The first key ingredient isthe construction, for any point a in the algebraic closureof F q that is not a pole of A ( x ), of a matrix Y a with entriesin ℓ = F q ( a ) which is similar to the evaluation A p ( a ) of the p -curvature at the point a . This construction comes from [4]and ultimately relies on the existence of a well-suited ring, ofso-called Hurwitz series in x − a , for which an analogue of theCauchy–Lipschitz theorem holds for the system Y ′ = A ( x ) Y around the (ordinary) point x = a . The matrix Y a is the p -th coefficient of the fundamental matrix of Hurwitz seriessolutions of Y ′ = A ( x ) Y at x = a .The second key ingredient is a baby step / giant step al-gorithm that computes Y a in O ˜( √ p ) operations in ℓ via fastmatrix factorials. Finally, we recover the invariant factorsof A p from those of the matrices Y a , for a suitable numberof values a . The main difficulty in this interpolation pro-cess is that there exist badly behaved points a for which theinvariant factors of A p ( a ) are not the evaluations at a ofthe invariant factors of A p ( x ). The remaining task is thento bound the number of unlucky evaluation points a . Thekey feature allowing a good control on these points, inde-pendent of p , is the fact that the invariant factors of A p ( x )have coefficients in F q ( x p ). Relationship to previous work.
There exists a largebody of work concerning the computation of so-called Frobe-nius forms of matrices (that is, the list of their invariantfactors, possibly with corresponding transformation matri-ces), and the related problem of Smith forms of polynomialmatrices. The specificities of our problem prevent us fromapplying these methods directly; however, our work is re-lated to several of these previous results.Let ω be a feasible exponent for matrix multiplication.The best deterministic algorithm known so far for the com-putation of the Frobenius form of an n × n matrix over afield k is due to Storjohann [30]. This algorithm has run-ning time O ( n ω log( n ) log log( n )) operations in k . We willuse it to compute the invariant factors of the matrices Y a above. Las Vegas algorithms were given by Giesbrecht [19],Eberly [14] and Pernet and Storjohann [28], the latter havingexpected running time O ( n ω ) over sufficiently large fields.The case of matrices with integer or rational entries hasattracted a lot of attention; this situation is close to ours,with the bit size of integers playing a role similar to the de-gree of the entries in the p -curvature. Early work goes backto algorithms of Kaltofen et al. [23,24] for the Smith form ofmatrices over Q [ x ], which introduced techniques used in sev-eral further algorithms, such as the Las Vegas algorithm byStorjohann and Labahn [31]. Giesbrecht’s PhD thesis [18]gives a Las Vegas algorithm with expected cost O ˜( n ω +2 d )for the Frobenius normal form of an n × n matrix with integerentries of bit size d ; Storjohann and Giesbrecht substantiallyimproved this result in [20], with an algorithm of expectedcost O ˜( n d + n d ). The best Monte Carlo running timeknown to us is O ˜( n . d ), by Kaltofen and Villard [25].In the latter case of matrices with integer coefficients, acommon technique relies on reduction modulo primes, anda main source of difficulty is to control the number of “un-lucky” reductions. We pointed out above that this is thecase in our algorithm as well. In general, the number of un- lucky primes is showed to be O ˜( n d ) in [18]; in our case, thedegree d of the entries grows linearly with p , but as we saidabove, we can alleviate this issue by exploiting the proper-ties of the p -curvature. Storjohann and Giesbrecht provedin [20] that a candidate for the Frobenius form of an integermatrix can be verified using only O ˜( nd ) primes; it wouldbe most interesting to adapt this idea to our situation. Structure of the paper.
In Section 2, we recall themain theoretical properties of the invariant factors of a poly-nomial matrix, and study their behavior under specializa-tion. We obtain bounds on bad evaluation points, and usethem to design (deterministic and probabilistic) evaluation-interpolation algorithms for computing the invariant factorsof a polynomial matrix. Section 3 is devoted to the de-sign of our main algorithms for the similarity class of the p -curvature, with deterministic and probabilistic versions forboth the system case and the scalar case. Finally, Section 4presents an application of our algorithm, that allows to an-swer a question coming from theoretical physics. Complexity basics.
We use standard complexity notation,such as ω for the exponent of matrix multiplication. Thebest known upper bound is ω < . d in k [ x ] can be performed in O ˜( d ) operations in the field k :addition, multiplication, shift, interpolation, etc , the key tothese results being fast polynomial multiplication [29,11,21].A general reference for these questions in [17].
2. COMPUTING INVARIANT FACTORS OFSPECIAL POLYNOMIAL MATRICES2.1 Definition and classical facts
We recall here some basic facts about invariant factors ofmatrices defined over a field. We fix for now a field K , anda matrix M ∈ M n ( K ). For a monic polynomial P = T d − P d − i =0 a i T i ∈ K [ T ], let M P denote its companion matrix: M P = a a .. . ...1 a d − . A well-known theorem [16, Th. 9, Ch. VII] asserts thatthere exist a unique sequence of monic polynomials I , . . . , I n for which I j divides I j +1 for all j and M is similar to a blockdiagonal matrix whose diagonal entries are M I , . . . , M I n .The I j ’s are called the invariant factors of M . We empha-size that, with our convention, there are always n invariantfactors but some of them may be equal to 1, in which casethe corresponding companion matrix is the empty one. Un-der this normalization, the j -th invariant factor I j can beobtained as I j = G j /G j − , where G j is the gcd of the mi-nors of size j of the matrix T I n − M , where I n stands for theidentity matrix of size n . The invariant factors are closelyrelated to the characteristic polynomial; indeed, we have I · I · · · I n = G n = det( T I n − M ) . (1)Given some irreducible polynomial P in K [ T ], we considerthe sequence (of integers): e d P,e = dim K ker P e ( M )deg P . (2)t turns out that this sequence completely determines the P -adic valuation of the invariant factors. Indeed, denotingby v j the P -adic valuation of I j , we have the relations: d P,e = n X j =1 min( e, v j ) , (3) d P,e − d P,e − = Card { j | v j ≥ e } (4)from which the v j ’s can be recovered without ambiguitysince they form a nondecreasing sequence. It also followsfrom the above formula that the sequence e d P,e is con-cave and eventually constant. Its final value is the dimen-sion of the characteristic subspace associated to P and it isreached as soon as e is greater than or equal to v n . Let k be a perfect field of characteristic p . We consider amatrix M ( x ) with coefficients in k [ x ]. For an element a lyingin a finite extension ℓ of k , we denote by M ( a ) the image of M ( x ) under the mapping k [ x ] → ℓ , x a . Our aim is tocompare the invariant factors of M ( x ) and those of M ( a ).We introduce some notation. Let I ( x, T ) , . . . , I n ( x, T )be the invariant factors of M ( x ). It follows from the re-lation (1) that they all lie in k [ x, T ]. We can thereforeevaluate them at x = a for each element a ∈ ℓ as aboveand get this way univariate polynomials with coefficientsin ℓ . Let I ( a, T ) , . . . , I n ( a, T ) be these evaluations. Wealso consider the invariant factors of M ( a ) and call them I ,a ( T ) , . . . , I n,a ( T ). We furthemore define G j ( x, T ) = I ( x, T ) · I ( x, T ) · · · I j ( x, T )and G j,a ( T ) = I ,a ( T ) · I ,a ( T ) · · · I j,a ( T ) . The characterization of the G j ’s in term of minors yields: Lemma For all a ∈ ℓ and all j ∈ { , . . . , n } , the poly-nomial G j ( a, T ) divides G j,a ( T ) in ℓ [ T ] . Let P ( x, T ) , . . . , P s ( x, T ) be the irreducible factors of thecharacteristic polynomial χ ( x, T ) of M ( x ), and let us write χ sep ( x, T ) for P ( x, T ) · · · P s ( x, T ). For all 1 ≤ i ≤ s and1 ≤ j ≤ n , let e i,j be the multiplicity of P i ( x, T ) in I j ( x, T ). Proposition We assume χ sep ( a, T ) is separable and dim k ( x ) ker P i ( x, M ( x )) e i,j +1 = dim ℓ ker P i ( a, M ( a )) e i,j +1 for all i and for all j < n . Then I j ( a, T ) = I j,a ( T ) for all j . Proof.
The equality of dimensions is also true for j = n ,as their sum on both sides equals n (using separability) andthese dimensions can only increase by specialization. Let d P i ,e be the sequence defined by Eq. (2) with respect to theirreducible polynomial P i ( x, T ) and the matrix M ( x ). Wedefine similarly for each irreducible factor P ( T ) of P i ( a, T )the sequence d P,e corresponding to the polynomial P ( T ) andthe matrix M ( a ). We claim that it is enough to prove that d P i ,e = d P,e for all e , i and all irreducible divisors P ( T ) of P i ( a, T ). Indeed, by Eq. (4), such an equality would imply: v P ( T ) ( I j,a ( T )) = e i,j (5)provided that P ( T ) is an irreducible divisor of P i ( a, T ), andwhere v P ( T ) denotes the P ( T )-adic valuation. On the otherhand, still assuming that P ( T ) is an irreducible divisor of P i ( a, T ), it follows from the definition of the e i,j ’s that: v P ( T ) ( I j ( a, T )) ≥ e i,j (6) and that the equality holds if and only if P ( T ) does not di-vide any of the P i ′ ( a, T ) for i ′ = i . Comparing characteristicpolynomials, we know moreover that P nj =1 v P ( T ) ( I j,a ( T )) = P nj =1 v P ( T ) ( I j ( a, T )). Combining this with (5) and (6), wefind that the P i ( a, T )’s are pairwise coprime and finally get I j ( a, T ) = I j,a ( T ) for 1 ≤ j ≤ n , as wanted.Until the end of the proof, we fix the index i and reservethe letter P to denote an irreducible divisor of P i ( a, T ). Fora fixed integer e , denote by j the greatest index j for which v P ( T ) ( I j,a ( T )) < e and observe that Eq. (3) can be rewritten d P,e = e · (cid:0) n − j (cid:1) + v P ( T ) (cid:0) G j ,a ( T ) (cid:1) . Using Lemma 1, wederive d P,e ≥ e · (cid:0) n − j (cid:1) + v P ( T ) (cid:0) G j ( a, T ) (cid:1) ≥ d P i ,e forall P and e . Eq. (4) now implies that the indices e forwhich d P i ,e − d P i ,e − > d P i ,e +1 − d P i ,e are exactly the e i,j ’s(1 ≤ j ≤ n ). Using concavity, we then observe that it isenough to check that d P i ,e = d P,e for indices e of the form e i,j + 1. For those e , we have by assumption: P P deg P · d P,e = dim ℓ ker P i ( a, M ( a )) e = dim k ( x ) ker P i ( x, M ( x )) e = deg T P i · d P i ,e = P P deg P · d P i ,e and thus d P,e = d P i ,e for all P because the inequalities d P,e ≥ d P i ,e are already known. Let M ( x ) be a square matrix of size n with coefficientsin k [ x ]. We set X = x p and assume that:(i) the entries of M ( x ) have degree at most pm (for a m ∈ N ),(ii) M ( x ) is similar to a matrix with coefficients in k ( X ).We are going to bound the number of values of a forwhich the invariant factors of M ( x ) do not specialize cor-rectly at x = a . Similar discussions appear is Section 4 ofGiesbrecht’s thesis [18] in the (more complicated) case ofinteger matrices. Our treatment is nevertheless rather dif-ferent in many places. The basic bound.
By assumption (ii), the characteristicpolynomial χ ( x, T ) lies in the subring k [ X, T ] of k [ x, T ]. Lemma The invariant factors I j ( x, T ) all belong to k [ X, T ] . Their degree with respect to X is at most mn . Proof.
By assumption (i), χ ( x, T ) is a polynomial in x of degree at most pmn . It then follows from Eq. (1) thatthe I j ( x, T )’s are polynomials in x of degree at most pmn as well. Now, the assumption (ii) ensures that the I j ( x, T )’sactually lie in k ( X )[ T ]. This completes the proof. Lemma We assume that p > n . There are at most deg X χ ( x, T ) · (2 n − points a ∈ k such that at least one ofthe P i ( a, T ) ’s is not separable. Proof.
We have that deg X χ sep ( x, T ) ≤ deg X χ ( x, T )and deg T χ sep ( x, T ) ≤ n , since χ sep divides χ . Denote by D ( x ) the discriminant of χ sep ( x, T ) with respect to T . Itsdegree in X is at most deg X χ ( x, T ) · (2 n − p > n implies that D ( x ) is not identically zero. Forany a ∈ k such that D ( a p ) = 0, the polynomial χ sep ( a p , T )is separable, and the same holds for the P i ( a p , T )’s. Notingthat k is perfect, the conclusion holds. Proposition We assume p > n . Let a , . . . , a N beelements in a separable closure of k which are pairwise nononjugate over k . We assume that for each i ∈ { , . . . , N } ,there exists j ∈ { , . . . , n } with I j ( a i , T ) = I a i ,j ( T ) . Then: N X i =1 deg( a i ) ≤ mn · ( n −
1) + mn · (2 n − where deg( a i ) denotes the algebraicity degree of a i over k . Proof.
We use the criteria of Proposition 2. We startby putting away the values of a for which at least one ofthe P i ( a, T )’s is not separable. By Lemma 4, there are atmost mn · (2 n −
1) such values. We then have to bound fromabove the values of a such that the equalities:dim k ( x ) ker P i ( x, M ( x )) e = dim ℓ ker P i ( a, M ( a )) e may fail for some i and some exponent e = e i,j +1 for some j .Let us fix such a pair ( i, e ). Set N ( x ) = P i ( x, M ( x )) e for simplicity. By assumption (i), the entries of N ( x ) havedegree at most pm i,e with m i,e = e · (cid:0) m deg T P i + deg X P i (cid:1) .On the other hand, we deduce from assumption (ii) that the P i ( x, T )’s all lie in k [ X, T ] and, as a consequence, that N ( x )is similar to a matrix with coefficients in k ( X ). Define d =dim k ( x ) ker N ( x ). The equality dim ℓ ker N ( a ) = d then failsif and only if the minors of N ( x ) of size n − d all vanish at x = a , i.e., if and only if the gcd ∆( x ) of these minors is divisibleby the minimal polynomial of a over k , say π a ( x ). Notingthat ∆( x ) ∈ k [ X ], the latter condition is also equivalent tothe fact that π a ( x ) p divides ∆( x ) in the ring k [ X ]. Thiscan be possible for at most deg X ∆( x ) ≤ ( n − d ) m i,e ≤ ( n − m i,e values of a .Therefore, if a , . . . , a N are pairwise non-conjugate “un-lucky values” of a , the sum appearing in the statement ofthe proposition is bounded from above by:( n − P i,e m i,e = m ( n − P i,e e deg T P i + ( n − P i,e e deg X P i . We notice that, when i remains fixed, the number of expo-nents of the form e i,j + 1 (1 ≤ j < n ) is bounded from aboveby e i,n + 1. The sum of these exponents is then at most: (cid:0) P n − j =1 e i,j (cid:1) + e i,n + 1 = e i + 1 ≤ e i , where e i denotes the multiplicity of the factor P i ( x, T ) in thecharacteristic polynomial χ ( x, T ). Our bound then becomes2 m ( n −
1) deg T χ + 2( n −
1) deg X χ . Using deg T χ = n anddeg X χ ≤ mn yields the bound. A refinement.
For the applications we have in mind, weshall need a refinement of Proposition 5 under the followinghypothesis depending on a parameter µ ∈ N :( H µ ): the polynomial χ has degree at most pµ w.r.t x .We observe that ( H µ ) is fulfilled when M ( x ) is a companionmatrix whose entries are polynomials of degree at most pµ . Proposition Under the assumptions of Prop. 5 andthe additional hypothesis ( H µ ) , we have: N X i =1 deg( a i ) ≤ µ · (2 n −
1) + µ · (2 n − . Proof.
Let P ( x, T ) be any bivariate polynomial with co-efficients in k . Set N ( x ) = P ( x, M ( x )) and let δ ( x ) denote the gcd of the minors of size s (for some integer s ) of N ( x ).We claim that:deg x δ ( x ) ≤ pµ · deg T P + s · deg x P (7)To prove the claim, we consider the Frobenius normal form˜ M ( x ) of M ( x ) and set ˜ N ( x ) = P ( x, ˜ N ( x )). Observe thatany minor of ˜ M ( x ) vanishes or has the shape ± c ( x ) · · · c n ( x )where c j ( x ) is a coefficient of I j ( x, T ) for all j . Noting thatdeg x I + · · · + deg x I n = deg x χ ≤ pµ , we derive that all the minors of ˜ M ( x ) have degree at most pµ . Now write P ( x, T ) = P deg T Pj =0 a j ( x ) T j where the a i ( x )’s lie in k [ x ]. Let˜ f denote the k [ x ]-linear endomorphism of k [ x ] n attached tothe matrix ˜ M ( x ). Set ˜ g = P ( x, ˜ f ); it clearly corresponds to˜ N ( x ). Given a vector space E and s linear endomorphisms u , . . . , u s of E , let us agree to define u ∧ · · · ∧ u s as E ⊗ s → V s Ex ⊗ · · · ⊗ x s u ( x ) ∧ · · · ∧ u s ( x s ) . where V s E is here defined as a quotient of E ⊗ s . Expandingthe exterior product V s ˜ g , we get: V s ˜ g = deg T P X i ,...,i s =0 a i ( x ) · · · a i s ( x ) · ˜ f i ∧ · · · ∧ ˜ f i s . (8)Moreover, assuming for simplicity that i ≤ i ≤ · · · ≤ i s and letting i = 0 by convention, we can write:˜ f i ⊗ · · · ⊗ ˜ f i s = (cid:13) sj =0 (cid:2) ( N j id) ⊗ ( N s − j ˜ f ) i j − i j − (cid:3) , where (cid:13) denotes the composition of the above (pairwisecommuting) maps. We get that the entries of the matrix (inthe canonical basis) of ˜ f i ∧ · · · ∧ ˜ f i s all have degree at most pµ · i s . The same argument demonstrates that the degreesof the entries of the above matrix are not greater than: pµ · max( i , . . . , i s ) ≤ pµ · deg T P when we no longer assume that the i j ’s are sorted by nonde-creasing order. Therefore, back to Eq. (8), we find that theentries of V s ˜ N ( x ) have degree at most pµ · deg T P + s · deg x P .It is then also the case of its trace, which is the same as thetrace of V s N ( x ) since N ( x ) and ˜ N ( x ) are similar. This fi-nally implies the claimed inequality (7) because δ ( x ) has todivide this trace.The Proposition now follows by inserting the above inputin the proof of Proposition 5. We keep the matrix M ( x ) satisfying the assumptions (i)and (ii) of §2.3. From now on, we assume that the onlyaccess we have to the matrix M ( x ) passes through a blackbox invariant_factors_at M ( x ) that takes as input an ele-ment a lying in a finite extension ℓ of k and outputs instantlythe invariant factors I j,a ( T ) of the matrix M ( a ). Our aim isto compute the invariant factors of M ( x ). We will proposetwo possible approaches: the first one is deterministic butrather slow although the second one is faster but probabilis-tic and comes up with a Monte-Carlo algorithm which maysometimes output wrong answers.Throughout this section, the letter D refers to a priori up-per bound on the X -degree of the characteristic polynomialof M ( x ). One can of course always take D = mn but bet-ter bounds might be available in particular cases. Similarlywe reserve the letter F for an upper bound on the sum ofegrees of “unlucky evaluation points”. Proposition 5 tellsus that mn (6 n −
5) is always an acceptable value for F . Re-member however that this value can be lowered to 3 µ (2 n − H µ ) thanks to Proposition 6. We willalways assume that F ≥ D .For simplicity of exposition, we assume from now on that k = F q is a finite field of cardinality q (it is more difficultand the case of most interest for us). Deterministic.
The discussion of §2.3 suggests the follow-ing algorithm whose correctness follows directly from thedefinition of F together with the assumption F ≥ D . Algorithm invariant_factors_deterministic
Input: M ( x ) satisfying (i) and (ii), D , F with F ≥ D Output:
The invariant factors of M ( x )1. Construct an extension ℓ of F q of degree F + 1and pick an element a ∈ ℓ such that ℓ = F q [ a ] Cost: O ˜( F ) operations in F q I ,a ( T ) , . . . , I n,a ( T ) = invariant_factors_at M ( x ) ( a )3. for j = 1 , . . . , n
4. Find I j ( x, T ) of degree ≤ D s.t. I j ( a, T ) = I j,a ( T )5. return I ( x, T ) , . . . , I n ( x, T ) Proposition The algorithm above requires only onecall to the black box invariant_factors_at M ( x ) with an in-put of degree exactly F + 1 . Probabilistic.
We now present a Monte-Carlo algorithm:
Algorithm invariant_factors_montecarlo
Input: M ( x ) s.t. (i) and (ii), ε ∈ (0 , D , F with F ≥ D Output:
The invariant factors of M ( x )1. Find the smallest integer s such that:2 · ( D + s +1) s ( q s − F ) + 12 · (cid:16) Fq s (cid:17) ( D − /s ≤ ε (9)and set K = ⌈ Ds ⌉ and k = ⌈ D +1 s ⌉ .2. for i = 1 , . . . , K
3. pick at random a i ∈ F q s s.t. F q s = F q [ a i ] Cost: O ˜( s ) operations in F q I ,i ( T ) , . . . , I n,i ( T ) = invariant_factors_at M ( x ) ( a i )5. for j = 1 , . . . , n d j = max i deg( I ,i ( T ) · I ,i ( T ) · · · I j,i ( T ))7. select I ⊂ { , . . . , K } of cardinality k s.t.(i) deg( I ,i ( T ) · I ,i ( T ) · · · I j,i ( T )) = d j for all i ∈ I (ii) the a i are pairwise non conjugate for i ∈ I Remark: if such I does not exist, raise an error8. compute I j ∈ F q [ X, T ] of X -degree ≤ D s.t. I j ( a i , T ) = I j,i ( T ) for all i ∈ I Cost: O ˜( D ) operations in F q return I ( x, T ) , . . . , I n ( x, T ) Proposition We have s ∈ O (log F Dε ) . Moreover: • Correctness:
Algorithm invariant_factors_montecarlo fails or returns a wrong answer with probability at most ε . • Complexity:
It performs ⌈ Ds ⌉ calls to the black box withinputs of degree s and O ˜( n ( D + log Fε )) operations in F q . Proof.
The first assertion is left to the reader. Let A be the set of elements a of F q s such that F q [ a ] = F q s . It isan easy exercise to prove that A has at least q s elements(the bound is not sharp). Let C , . . . , C C be the conjugacyclasses (under the Galois action) in A . Remark that each C i has by definition s elements, so that C ≥ q s s . We say that aconjugacy class is bad if it contains one element a for which I j ( a, T ) = I a,j ( T ) for some j . Otherwise, we say that it is good . Let B (resp. G ) be the number of bad (resp. good)classes. We have B + G = C and B ≤ Fs by definition of F .The algorithm invariant_factors_montecarlo succeedsif there exist at least k indices i for which the corresponding a i ’s lie in pairwise distinct good classes. This happens withprobability at least:1 C K · K X j = k (cid:0) Kj (cid:1) · G ( G − · · · ( G − k + 1) · G j − k · B K − j . (The above formula gives the probability that the first k good classes are pairwise distinct, which is actually strongerthan what we need.) The above quantity is at least equal to (cid:16) − kG (cid:17) k · − k − X j =0 (cid:0) Kj (cid:1) · (cid:16) GC (cid:17) j · (cid:16) BC (cid:17) K − j ! . Moreover for j ≤ k −
1, we have: (cid:16) GC (cid:17) j · (cid:16) BC (cid:17) K − j ≤ (cid:16) BGC (cid:17) j · (cid:16) BC (cid:17) K − j ≤ j · (cid:16) Fq s (cid:17) K − j ≤ K · (cid:16) Fq s (cid:17) K − j ≤ K · (cid:16) Fq s (cid:17) ( D − /s . Therefore the probability of success is at least: (cid:16) − kG (cid:17) k · (cid:18) − · (cid:16) Fq s (cid:17) ( D − /s (cid:19) . Using k ≤ D + s +1 s and G ≥ q − F s , we find that the probabil-ity of failure is at most the LHS of Eq. (9). The correctnessis proved. As for the complexity, the results are obvious.
3. COMPUTING INVARIANT FACTORSOF THE P-CURVATURE
Throughout this section, we fix a finite field k = F q ofcardinality q and characteristic p . We endow the field ofrational functions k ( x ) with the natural derivation f f ′ . We recall that a differential module over k ( x ) is k ( x )-vector space M endowed with an additive map ∂ : M → M satisfying the following Leibniz rule: ∀ f ∈ k ( x ) , ∀ m ∈ M, ∂ ( fm ) = f ′ · m + f · ∂ ( m ) . The p -curvature of a differential module M is the mapping ∂ p = ∂ ◦ · · · ◦ ∂ ( p times). Using the fact that the p -thderivative of any f ∈ k ( x ) vanishes, we derive from theLeibniz relation above that ∂ p is k ( x )-linear endomorphismof M . It follows moreover from [4, Remark 4.5] that ∂ p isdefined over k ( x p ), in the sense that there exists a k ( x )-basisof M in which the matrix of ∂ p has coefficients in k ( x p ). Inparticular, all the invariant factors of the p -curvature havetheir coefficients in k ( x p ). tatement of the main Theorem. From now on, we fixa differential module (
M, ∂ ). We assume that M is finitedimensional over k ( x ) and let r denote its dimension. Wepick ( e , . . . , e r ) a basis of M and let A denote the matrix of ∂ with respect to this basis. We write A = f A ˜ A where f A and the entries of ˜ A all lie in k [ x ]. Let d be an upper boundon the degrees of all these polynomials. The aim of thissection is to design fast deterministic and probabilistic algo-rithms for computing the invariant factors of the p -curvatureof ( M, ∂ ). The following Theorem summarizes our results.
Theorem We assume p > r .1. There exists a deterministic algorithm that computes theinvariant factors of the p -curvature of ( M, ∂ ) within O ˜ (cid:0) d ω + r ω +2 √ p (cid:1) operations in k = F q .2. Let ε ∈ (0 , . There exists a Monte-Carlo algorithm com-puting the invariant factors of the p -curvature of ( M, ∂ ) in O ˜ (cid:0) d ω + r ω · ( dr − log ε ) · √ p (cid:1) operations in k = F q . This algorithm returns a wrong an-swer with probability at most ε . In what follows, we will use the notation A p ( x ) for thematrix of the p -curvature of ( M, ∂ ) with respect to the dis-tinguished basis ( e , . . . , e r ). Given an element a lying ina finite extension ℓ of k , we denote by A p ( a ) ∈ M r ( ℓ ) thematrix deduced from A p by evaluating it at x = a . The similarity class of A p ( a ) . Let S be an irreduciblepolynomial over k . Set ℓ = k [ u ] /S and let a denote theimage of the variable u in ℓ . We assume that S does notdivide f A , i.e., f A ( a ) = 0. The first ingredient we need isthe construction of an auxiliary matrix which is similar to A p ( a ). This construction comes from our previous paper [4].Let us recall it briefly. We define the ring ℓ [[ t ]] dp of Hurwitzseries whose elements are formal infinite sums of the shape: a + a γ ( t ) + a γ ( t ) + · · · + a n γ n ( t ) + · · · (10)and on which the addition is straightforward and the multi-plication is governed by the rule γ i ( t ) · γ j ( t ) = (cid:0) i + ji (cid:1) γ i + j ( t ).(The symbol γ i ( t ) should be thought of as t i i ! .) We moreoverendow ℓ [[ t ]] dp with the derivation defined by γ i ( t ) ′ = γ i − ( t )(with the convention that γ ( t ) = 1) and the projection mappr : ℓ [[ t ]] dp → ℓ sending the series given by Eq. (10) toits constant coefficient a . We shall often use the alterna-tive notation f (0) for pr( f ). If f ∈ ℓ [[ t ]] dp is given by theseries (10), we then have a n = f ( n ) (0) for all nonnegativeintegers n . We have a homomorphism of rings: ψ S : k [ x ][ f A ] → ℓ [[ t ]] dp , f ( x ) p − X i =0 f ( i ) ( a ) γ i ( t ) . It is easily checked that ψ S commutes with the derivation.We can then consider the differential module over ℓ [[ t ]] dp obtained from ( M, ∂ ) by scalar extension. By definition, itcorresponds to the differential system Y ′ = ψ S ( A ) · Y .The benefit of working over ℓ [[ t ]] dp is the existence of ananalogue of the well-known Cauchy–Lipschitz Theorem [4,Proposition 3.4]. This notably implies the existence of afundamental matrix of solutions, i.e., an r × r matrix Y S with entries in ℓ [[ t ]] dp , and satisfying: Y ′ S = ψ S ( A ) · Y S and Y S (0) = I r (11) with I r the identity matrix of size r . Moreover, as explainedin more details later, the construction of Y S is effective.For any integer n ≥
0, we let Y ( n ) S denote the matrixobtained from Y S by taking the n -th derivative entry-wise.The next proposition is a consequence of [4, Proposition 4.4]. Proposition
The matrices A p ( a ) and − Y ( p ) S (0) aresimilar over ℓ . Fast computation of Y ( p ) S (0) . We recall that Y S is definedas the solution of the system (11). Remembering that wehave written A = f A ˜ A , we obtain the relation: ψ S ( f A ) · Y ′ S = ψ S ( ˜ A ) · Y S . (12)Write f A = P di =0 f i · ( x − a ) i and ˜ A = P di =0 ˜ A i · ( x − a ) i wherethe f i ’s lie in ℓ and the A i ’s are square matrices of size r withentries in ℓ . Remark that f does not vanish because it isequal to f A ( a ). Note moreover that the f i ’s can be computedfor a cost of O ˜( d ) operations in k using divide-and-conquertechniques. Given a fixed pair of indices ( i ′ , j ′ ), the samediscussion applies to the collection of the ( i ′ , j ′ )-entries ofthe A i ’s. The total cost for computing the decompositionsof f A and ˜ A is then O ˜( dr ). Now, coming back to thedefinitions, we find that ψ S ( f A ) = P di =0 i ! f i · γ i ( t ) and ψ S ( ˜ A ) = P di =0 i ! ˜ A i · γ i ( t ). Eq. (12) yields the recurrence: Y ( n +1) S (0) = min( n,d ) X i =0 B i ( n ) · Y ( n − i ) S (0) (13)where the B i ∈ M r ( ℓ [ u ]) are defined by: f B i = u ( u − · · · ( u − i +1) · (cid:0) ˜ A i − ( u − i ) f i +1 · I r (cid:1) (14)with the convention that f d +1 = 0. Now setting: Z n = Y ( n − d ) S (0) Y ( n − d +1) S (0)... Y ( n ) S (0) , B = I r I r . .. I r B d · · · · · · · · · B (with the convention Y ( i ) s (0) = 0 when i < Z n +1 = B ( n ) · Z n . Hence, we obtain Z p = B ( p − · B ( p − · · · B (0) · Z from what we finally getthat Y ( p ) S (0) is the ( r × r )-matrix located at the bottom rightcorner of B ( p − · B ( p − · · · B (0). The computation of theformer matrix factorial can be performed efficiently using avariation of the Chudnovskys’ algorithm [12, 6]. Combiningthis with Proposition 10, we end up with the following. Proposition
The invariant factors of A p ( a ) can becomputed in O ˜( d ω r ω √ dp ) operations in the field ℓ . Proof.
Note that B is a square matrix of size ( d +1) r .Moreover coming back to (14), we observe that the entriesof B all have degree at most d . By [5, Theorem 2] the matrixfactorial − B ( p − · B ( p − · · · B (0) can then be computedfor the cost of O ˜( d ω r ω √ dp ) operations in ℓ . By [30], theinvariant factors of its submatrix − Y ( p ) S (0) can be obtainedfor an extra cost of O ˜( r ω ) operations in ℓ (which is negli-gible compared to the previous one). Using Proposition 10these invariant factors are also those of A p ( a ), and we aredone. onclusion. Proposition 11 yields an acceptable primitive invariant_factors_at A p ( x ) . Plugging it in the algorithm invariant_factors_deterministic and using the parame-ters D = dr and F = 6 dr ( r − A p ( x ) forthe cost of one unique call to invariant_factors_at A p ( x ) with an input lying in an extension ℓ/k of degree F + 1 (cf.Proposition 7). By Proposition 11, we find that the totalcomplexity of the obtained algorithm is O ˜ (cid:0) d ω + r ω +2 √ p (cid:1) operations in F q . The first part of Theorem 9 is then estab-lished. The second part is obtained in a similar fashion us-ing the algorithm invariant_factors_montecarlo togetherwith Proposition 8 for correctness and complexity results. The ring of differential operators k ( x ) h ∂ i is the ring ofusual polynomials over k ( x ) in the variable ∂ except thatthe multiplication is ruled by the relation ∂ · f = f · ∂ + f ′ .We define similarly the ring k [ x ] h ∂ i . We say that L ∈ k [ x ] h ∂ i has bidegree ( d, r ) if it has degree d with respect to x anddegree r with respect to ∂ .If L is a differential operator in k ( x ) h ∂ i , one easily checksthat the set k ( x ) h ∂ i L of left multiples of L is a left ideal of k ( x ) h ∂ i . The quotient M L = k ( x ) h ∂ i /k ( x ) h ∂ i L is then avector space over k ( x ). It is moreover endowed with a map ∂ : M L → M L given by the left multiplication by ∂ . Thismap turns M L into a differential module.We shall prove in this section that the complexities an-nounced in Theorem 9 can be improved in the case of differ-ential modules coming from differential operators. Below isthe statement of our precise result. Theorem
Let L ∈ k [ x ] h ∂ i be a differential operatorof bidegree ( d, r ) . We assume p > r .1. There exists a deterministic algorithm that computes theinvariant factors of the p -curvature of M L within O ˜ (cid:0) ( d + r ) ω +1 d r · √ p (cid:1) operations in k = F q .2. Let ε ∈ (0 , . There exists a Monte-Carlo algorithm thatcomputes the invariant factors of the p -curvature of M L in O ˜ (cid:0) ( d + r ) ω d · ( d − log ε ) · √ p (cid:1) operations in k = F q . This algorithm returns a wrong an-swer with probability at most ε . Better bounds.
From now on, we fix a differential opera-tor L ∈ k ( x ) h ∂ i of bidegree ( d, r ). We denote by A p ( x ) thematrix of the p -curvature of M L with respect to the canon-ical basis (1 , ∂, . . . , ∂ r − ). If a r ( x ) is the leading coefficientof L (with respect to ∂ ), it follows from [13, Proposition 3.2]that A p ( x ) has the form A p ( x ) = a r ( x ) p · ˜ A p ( x ) where ˜ A p ( x )is a matrix with polynomial entries of degree at most pd . Proposition
The matrix ˜ A p ( x ) satisfies the hypoth-esis ( H r + d ) (introduced just before Proposition 6). Proof.
This is a direct consequence of Lemma 3.9 andTheorem 3.11 of [3].
The similarity class of A p ( a ) . We now revisit Proposi-tion 11 when the differential module comes from the differen-tial operator L . We fix an irreducible polynomial S ∈ k [ x ] and assume that S is coprime with the leading coefficient a r ( x ) of L . We set ℓ = k [ x ] /S and let a denote the image of x is ℓ . We define t = x − a ∈ ℓ [ x ] and consider the ring ofdifferential operators ℓ [ x ] h ∂ i . The latter acts on ℓ [[ t ]] dp byletting ∂ act as the derivation. Let Y S be the fundamentalsystem of solutions of the equation Y ′ S = ψ S ( A ) · Y S where A is the companion matrix which gives the action of ∂ on M L .It takes the form: Y S = y y · · · y r − y ′ y ′ · · · y ′ r − ... ... ... y ( r − y ( r − · · · y ( r − r − where y j ∈ ℓ [[ t ]] dp is the unique solution of the differen-tial equation Ly j = 0 with initial conditions y ( n ) j (0) = δ j,n (where δ · , · is the Kronecker symbol) for 0 ≤ n < r .We introduce the Euler operator θ = t · ∂ ∈ ℓ [ x ] h ∂ i . Usingthe techniques of [3, Section 4.1], one can write L · ∂ d = P d + ri =0 b i ( θ ) ∂ i within O ˜(( r + d ) d ) operations in ℓ . Here the b i ’s are polynomials with coefficients in ℓ of degree at most d .One can check moreover that the polynomial b d + r is constantequal to a r ( a ); in particular, it does not vanish thanks to ourassumption on S . For all j , define z j = P ∞ n =0 y ( n ) j (0) γ n + d ( t ).Clearly ∂ d z j = y j , so that we have (cid:16)P d + ri =0 b i ( θ ) ∂ i (cid:17) · z i = 0for all i . Noting that θ acts on γ n ( t ) by multiplication by n ,we get the recurrence relation: ∀ n ≥ , d + r X i =0 b i ( n ) · y ( n + i − d ) j (0) = 0with the convention that y ( n ) j = 0 when n <
0. Letting: Z n = y ( n − d )0 (0) · · · y ( n − d ) r − (0) y ( n − d +1)0 (0) · · · y ( n − d +1) r − (0)... ... Y ( n + r − (0) · · · y ( n + r − r − (0) ∈ M d + r,r ( ℓ )and B = − a r ( a ) · b b · · · b d + r − ∈ M d + r,d + r ( ℓ )the above recurrence rewrites Z n +1 = B ( n ) Z n . Solving therecurrence, we get Z p = B ( p − · · · B (0) · Z , and we derivethat Y ( p ) S (0) is the ( r × r ) matrix located at the bottom rightcorner of B ( p − · B ( p − · · · B (0). Using Proposition 10and [5, Theorem 2], we end up with the following proposition(compare with Proposition 11). Proposition
The invariant factors of A p ( a ) can becomputed in O ˜(( d + r ) ω √ dp ) operations in the field ℓ . Conclusion.
The final discussion is now similar to the onewe had in the case of differential modules. Proposition 14provides the primitive invariant_factors_at A p ( x ) . Usingit in the algorithms invariant_factors_deterministic and invariant_factors_montecarlo with the parameters D = d and F = 3 d (2 r −
1) (coming from the combination ofPropositions 6 and 13), we respectively end up with de-terministic and Monte-Carlo algorithms whose complexitiesagree with the ones announced in Theorem 12.t is instructive to compare the methods and results of thissection with those of our previous paper [3]. We remark thatthe matrix factorial considered above is nothing but the spe-cialization at θ = 0 of the matrix factorial in [3]. Althoughthe theoretical approaches of the two papers are definitivelydifferent, they lead to very similar computations. However,each of them has its own advantages and disadvantages. Onthe one hand, the methods of [3] deal with characteristicpolynomials only and cannot see invariant factors. On theother hand, they do not require the assumption a r ( a ) = 0(that is why we always took a = 0 in [3]) and can handleat the same time the local computations at the point a and around it, i.e., they provide roughly speaking a frameworkwhich allows to work modulo ( x − a ) pn for some integer n fixed in advance (not just n = 1) without increasing thecost with respect to p . The practical consequence is thatthe methods of the current paper end up with algorithmswhose cost is weakened by a factor √ d compared to whatwe might have expected at first. It would be interesting tofind a general theoretical setting unifying the two approachesdiscussed above and allowing the benefits of both of them.
4. SOLVING A PHYSICAL APPLICATION
In [10], a globally nilpotent differential operator φ (6) H wasintroduced in order to model the 6-particle contribution tothe square lattice Ising model. As shown in loc. cit. , this op-erator factors as a product of differential operators of smallerorders. The factor which is the least understood is called L and has order 23. Actually L has not been computed sofar because its size is too large. Nevertheless there exists amultiple of L which has a more reasonable size: its bide-gree is (140 , L , hasbeen determined modulo several prime numbers. Based onthis computation and using the strategy developed in thispaper, we were able to study a bit further the factorizationof L , answering a question of the authors of [10]. Proposition
The operator L cannot be factorizedas a product L · L where L is an operator of order , and L is an operator of order whose differential module isisomorphic to a symmetric product Sym n M for an integer n > and a differential module M . Proof.
We argue by contradiction by assuming that sucha factorization exists. This would imply that, for all p thematrix A ,p of the p -curvature of L mod p decomposes: A ,p = (cid:18) A ,p ⋆ A ,p (cid:19) (15)where A ,p (resp. A ,p ) is the square matrix of size 2(resp. 21) and A ,p is similar to a symmetric product of a d × d matrix A d,p . We now pick p = 32647 for which L mod p is known. Using Proposition 11, we were able to determinethe invariant factors of the p -curvature of A ,p (15). Thecomputation ran actually rather fast: just a few minutes.We observed that the generalized eigenspace of A ,p (15)for the eigenvalue 0 has dimension 23. Combining this withthat fact that L is a factor of L whose p -curvature isnilpotent, we deduce that the restriction of A ,p (15) to thischaracteristic space is similar to A ,p (15). Arguing simi-larly, we determine the Jordan type of A ,p (15): m ≥ (cid:0) A ,p (15) m (cid:1)
23 17 11 6 3 0 Moreover the writing (15) would imply that for all m :0 ≤ rank (cid:0) A ,p (15) m (cid:1) − rank (cid:0) A ,p (15) m (cid:1) ≤ (cid:0) A ,p (15) m (cid:1) = (cid:18) n − (cid:0) A d,p (15) m (cid:1) n (cid:19) . There is only one way to satisfy these numerical constraintswhich consists in taking n = 2 and: m ≥ (cid:0) A d,p (15) m (cid:1) (cid:0) A d,p (15) m (cid:1) − rank (cid:0) A d,p (15) m +1 (cid:1) has to be non-increasing, this is impossible.
5. REFERENCES [1] Y. André. Sur la conjecture des p -courbures de Grothendieck-Katz et unproblème de Dwork. In Geometric aspects of Dwork theory. Vol. I, II , pages55–112. Walter de Gruyter GmbH & Co. KG, Berlin, 2004.[2] A. Bostan, S. Boukraa, S. Hassani, J.-M. Maillard, J.-A. Weil, andN. Zenine. Globally nilpotent differential operators and the square Isingmodel.
J. Phys. A , 42(12):125206, 50, 2009.[3] A. Bostan, X. Caruso, and E. Schost. A fast algorithm for computing thecharacteristic polynomial of the p -curvature. In ISSAC’14 , pages 59–66.ACM, New York, 2014.[4] A. Bostan, X. Caruso, and E. Schost. A fast algorithm for computing the p -curvature. In ISSAC’15 , pages 69–76. ACM, New York, 2015.[5] A. Bostan, T. Cluzeau, and B. Salvy. Fast algorithms for polynomialsolutions of linear differential equations. In
ISSAC’05 , pages 45–52. ACMPress, 2005.[6] A. Bostan, P. Gaudry, and É. Schost. Linear recurrences with polynomialcoefficients and application to integer factorization and Cartier-Maninoperator.
SIAM Journal on Computing , 36(6):1777–1806, 2007.[7] A. Bostan and M. Kauers. Automatic classification of restricted latticewalks. In
FPSAC’09 , DMTCS Proc., AK, pages 201–215. 2009.[8] A. Bostan and M. Kauers. The complete generating function for Gesselwalks is algebraic.
Proc. Amer. Math. Soc. , 138(9):3063–3078, 2010. With anappendix by Mark van Hoeij.[9] A. Bostan and É. Schost. Fast algorithms for differential equations inpositive characteristic. In
ISSAC’09 , pages 47–54. ACM, New York, 2009.[10] S. Boukraa, S. Hassani, I. Jensen, J.-M. Maillard, and N. Zenine.High-order Fuchsian equations for the square lattice Ising model: χ (6). J.Phys. A , 43(11):115201, 22, 2010.[11] D. G. Cantor and E. Kaltofen. On fast multiplication of polynomials overarbitrary algebras.
Acta Inform. , 28(7):693–701, 1991.[12] D. V. Chudnovsky and G. V. Chudnovsky. Approximations and complexmultiplication according to Ramanujan. In
Ramanujan revisited(Urbana-Champaign, 1987) , pages 375–472. Academic Press, Boston, 1988.[13] T. Cluzeau. Factorization of differential systems in characteristic p . In ISSAC’03 , pages 58–65. ACM Press, 2003.[14] W. Eberly. Asymptotically efficient algorithms for the Frobenius form.Technical report, 2000.[15] F. L. Gall. Powers of tensors and fast matrix multiplication. In
ISSAC’14 ,pages 296–303, 2014.[16] F. R. Gantmacher.
The theory of matrices. Vols. 1, 2 . Translated by K. A.Hirsch. Chelsea Publishing Co., New York, 1959.[17] J. von zur Gathen and J. Gerhard.
Modern Computer Algebra . CambridgeUniversity Press, second edition, 2003.[18] M. Giesbrecht.
Nearly optimal algorithms for canonical matrix forms . PhD thesis,University of Toronto, 1993.[19] M. Giesbrecht. Nearly optimal algorithms for canonical matrix forms.
SIAM Journal on Computing , 24(5):948–969, 10 1995.[20] M. Giesbrecht and A. Storjohann. Computing rational forms of integermatrices.
Journal of Symbolic Computation , 34(3):157–172, 2002.[21] D. Harvey, J. van der Hoeven, and G. Lecerf. Faster polynomialmultiplication over finite fields. http://arxiv.org/abs/1407.3361, 2014.[22] T. Honda. Algebraic differential equations. In
Symposia Mathematica, Vol.XXIV , pages 169–204. Academic Press, London, 1981.[23] E. Kaltofen, M. S. Krishnamoorthy, and B. D. Saunders. Fast parallelcomputation of Hermite and Smith forms of polynomial matrices.
SIAMJournal on Matrix Analysis and Applications , 8(4):683–690, 1987.[24] E. Kaltofen, M. S. Krishnamoorthy, and B. D. Saunders. Parallelalgorithms for matrix normal forms.
Linear Algebra and its Applications ,136:189–208, 1990.[25] E. Kaltofen and G. Villard. On the complexity of computing determinants.
Comput. Complexity , 13(3-4):91–130, 2004.[26] N. M. Katz. Algebraic solutions of differential equations ( p -curvature andthe Hodge filtration). Invent. Math. , 18:1–118, 1972.[27] N. M. Katz. A conjecture in the arithmetic theory of differential equations.
Bull. Soc. Math. France , (110):203–239, 1982.[28] C. Pernet and A. Storjohann. Frobenius form in expected matrixmultiplication time over sufficiently large fields, 2007.[29] A. Schönhage. Schnelle Multiplikation von Polynomen über Körpern derCharakteristik 2.
Acta Informatica , 7:395–398, 1977.[30] A. Storjohann. Deterministic computation of the Frobenius form. In
FOCS’01 , pages 368–377. IEEE Computer Society Press, 2001.[31] A. Storjohann and G. Labahn. A fast Las Vegas algorithm for computingthe Smith normal form of a polynomial matrix.