A Fast Algorithm for Computing the p-Curvature
aa r X i v : . [ c s . S C ] J un A Fast Algorithm for Computing the p-Curvature
Alin Bostan
Inria (France) [email protected]
Xavier Caruso
Université Rennes 1 [email protected]
Éric Schost
Western University [email protected]
ABSTRACT
We design an algorithm for computing the p -curvature of adifferential system in positive characteristic p . For a systemof dimension r with coefficients of degree at most d , its com-plexity is O ˜( pdr ω ) operations in the ground field (where ω denotes the exponent of matrix multiplication), whereasthe size of the output is about pdr . Our algorithm is thenquasi-optimal assuming that matrix multiplication is ( i.e. ω = 2). The main theoretical input we are using is the exis-tence of a well-suited ring of series with divided powers forwhich an analogue of the Cauchy–Lipschitz Theorem holds. Categories and Subject Descriptors:
I.1.2 [
Computing Methodologies ]: Symbolic and Alge-braic Manipulation –
Algebraic Algorithms
Keywords:
Algorithms, complexity, differential equations, p -curvature.
1. INTRODUCTION
We study in this article algorithmic questions related tolinear differential systems in positive characteristic. Let k be an arbitrary field of prime characteristic p , and A bean r × r matrix with entries in the field k ( x ) of rationalfunctions over k . A simple-to-define, yet very importantobject attached to the differential system Y ′ = AY is itsso-called p -curvature . It is the p -th iterate ∂ pA of the map ∂ A : k ( x ) r → k ( x ) r that sends v to v ′ − Av . It turns outthat it is k ( x )-linear. It is moreover classical that its matrixwith respect to the canonical basis of k ( x ) r is equal to theterm A p of the recursive sequence ( A i ) i defined by A = − A and A i +1 = A ′ i − A · A i for i ≥ . (1)In all what follows, we will thus deliberately identify thematrix A p with the p -curvature of Y ′ = AY . The aboverecurrence yields an algorithm for computing it, sometimesreferred to as Katz’s algorithm .The p -curvature is related to solutions; it measures towhat extent the usual Cauchy–Lipschitz theorem applies in Publication rights licensed to ACM. ACM acknowledges that this contribution wasauthored or co-authored by an employee, contractor or affiliate of a national govern-ment. As such, the Government retains a nonexclusive, royalty-free right to publish orreproduce this article, or to allow others to do so, for Government purposes only.
ISSAC’15,
July 6–9, 2015, Bath, United Kingdom..Copyright is held by the owner/author(s). Publication rights licensed to ACM.ACM 978-1-4503-3435-8/15/07 ...$15.00.http://dx.doi.org/10.1145/2755996.2756674. characteristic p . More precisely, at an ordinary point, thesystem Y ′ = AY admits a fundamental matrix of powerseries solutions in k [[ x ]] if and only if the p -curvature A p vanishes. In this case, the system Y ′ = AY even admits afundamental matrix of solutions which are rational functionsin k ( x ). More generally, the dimension of the kernel of A p is equal to the dimension of the space of rational functionsolutions of Y ′ = AY .The primary importance of the notion of p -curvature re-lies in its occurrence in one of the versions of the celebratedGrothendieck–Katz conjecture [19, 20, 12, 30]. This con-jecture, first formulated by Alexandre Grothendieck in thelate 1960s, is a local-global principle for linear differentialsystems, which states that a linear differential system withrational function coefficients over a function field admits afundamental matrix of algebraic solutions if and only if its p -curvatures vanish for almost all primes p .In computer algebra, p -curvature has been introduced byvan der Put [22, 23], who popularized it as a tool for fac-toring differential operators in characteristic p . Cluzeau [13]generalized the approach to the decomposition of differen-tial systems over k ( x ). The p -curvature has also been usedby Cluzeau and van Hoeij [14] as an algorithmic filter forcomputing exponential solutions of differential operators incharacteristic zero.Computing efficiently the p -curvature is in itself a chal-lenging problem, especially for large values of p . Our initialmotivation for studying this question emerged from concreteapplications, in lattice path combinatorics [6, 7] and in sta-tistical physics [3]. In this article, we address the question ofthe computation of A p in good complexity, with respect tothree parameters: the dimension r of the system Y ′ = AY ,the maximum degree d of the rational function entries of A , and the characteristic p of the ground field. In terms ofthese quantities, the arithmetical size of A p is genericallyproportional to pdr if r > Previous work.
Cluzeau [13, Prop. 3.2] observed thatthe direct algorithm based on recurrence (1) has complexity O ˜( p dr ω ), where ω is the matrix multiplication exponentand the soft-O notation O ˜( ) hides polylogarithmic factors.Compared to the size of the p -curvature, this cost is goodwith respect to r and d , but not to p . The first subquadraticalgorithm in p , of complexity O ˜( p ω/ ), was designed in [9, § O ˜( p ) for certain systems of order r = 2. However, the ques-tion of designing a general algorithm for computing A p withquasi-linear complexity in p remained open. In a related, butifferent direction, the article [4] proposed an algorithm forcomputing the characteristic polynomial of the p -curvaturein time essentially linear in √ p , without computing A p itself. Contribution.
We prove that the p -curvature A p can becomputed in quasi-linear time with respect to p . More pre-cisely, our main result (Theorem 4.2) states that O ˜ (cid:0) pdr ω )operations in k are sufficient for this task. This complexityresult is quasi-optimal not only with respect to the main pa-rameter p , but also to d ; with respect to the dimension r , itis as optimal as matrix multiplication. Moreover the algo-rithm we obtain is highly parallelizable by design. The keytools underlying the proof of Theorem 4.2 are the notion ofdivided power rings in characteristic p , and a new formulafor the p -curvature (Propositions 4.3 and 4.4) in terms ofdivided power series. Crucial ingredients are the fact thata Cauchy–Lipschitz theorem for differential systems holdsover divided power rings (Proposition 3.4) and the fact thatNewton iteration can be used to efficiently compute (trunca-tions of) fundamental matrices of divided power solutions. Structure of the paper.
In Section 2, we recall the maintheoretical properties of the basic objects used in this article.Section 3 is devoted to the existence and the computationof solutions of differential systems in divided power rings.In Section 4, we move to the main objective of the arti-cle, the computation of the p -curvature: after relating A p to the framework of divided powers, we describe our mainalgorithm for A p , of complexity O ˜( pdr ω ). We conclude inSection 5 by describing the implementation of our algorithmand some benchmarks. Complexity measures.
Throughout this article, we es-timate the cost of our algorithms by counting arithmeticoperations in the base ring or field at unit cost.We use standard complexity notations. The letter ω refersto a feasible exponent for matrix multiplication ( i.e. thereexists an algorithm for multiplying n × n matrices over aring A with at most O ( n ω ) operations in A ); the best knownbound is ω < . O ˜( · )indicates that polylogarithmic factors are omitted; in par-ticular, we will use the fact that many arithmetic operationson univariate polynomials of degree d can be done in O ˜( d )operations: addition, multiplication, Chinese remaindering, etc , the key to these results being fast polynomial multiplica-tion [27, 26, 11, 18]. A general reference for these questionsis [17].
2. THEORETICAL SETTING
We introduce and briefly recall the main properties of thetheoretical objects we are going to use in this article. Allthe material presented in this section is classical; a generalreference is [24].
Definitions and notations.
Let A be a commutative ringwith unit. We recall that a derivation on A is an additivemap ′ : A → A , satisfying the Leibniz rule ( fg ) ′ = f ′ g + fg ′ for all f, g ∈ A . The image f ′ of f under the derivation iscalled the derivative of f . From now on, we assume that A is equipped with a derivation. A differential system withcoefficients in A is an equation of the form Y ′ = AY where A is a given r × r matrix with coefficients in A (for a certainpositive integer r ), the unknown Y is a column vector oflength r and Y ′ denotes the vector obtained from Y by tak-ing the derivative component-wise. The integer r is called the dimension of the system. We recall briefly that a lineardifferential equation: a r y ( r ) + · · · + a y ′ + a y = 0 (with a i ∈ A ) (2)can be viewed as a particular case of a differential system.Indeed, defining the companion matrix C = − a a r − a a r . .. ...1 − a r − a r (3)and A = t C , the solutions of the system Y ′ = AY areexactly the vectors of the form t ( y, y ′ , . . . , y ( r − ) where y is a solution of (2). In this correspondence, the order ofthe differential equation agrees with the dimension of theassociated differential system. Differential modules. A differential module over A is apair ( M, ∂ ) where M is an A -module and ∂ : M → M is anadditive map satisfying a Leibniz-like rule, which is: ∀ f ∈ A , ∀ x ∈ M, ∂ ( fx ) = f ′ · x + f · ∂ ( x ) . (4)There exists a canonical one-to-one correspondence betweendifferential systems and differential modules ( M, ∂ ) for which M = A r for some r : to a differential system Y ′ = AY of di-mension r , we attach the differential module ( A r , ∂ A ) where ∂ A : A r → A r is the function mapping X to X ′ − AX . Underthis correspondence, the solutions of Y ′ = AY are exactlyvectors in the kernel of ∂ A .To a differential equation as (2), one can associate the differential operator L = a r ∂ r + a r − ∂ r − + · · · + a ∂ + a ;it lies in the non-commutative ring A h ∂ i , endowed with theusual addition of polynomials and a multiplication ruled bythe relation ∂ · f = f · ∂ + f ′ for all f ∈ A (note that, as oftenin the literature, we are using ∂ to denote either the struc-ture map of a differential module, and a non-commutativeindeterminate).Then, if a r is a unit in A , one can further associate to L the quotient A h ∂ i / A h ∂ i L ≃ A r . The differential struc-ture inherited from A h ∂ i makes it a differential module withstructure map X ∈ A r X ′ + CX , where C is the compan-ion matrix defined above; in other words, this is the module( A h ∂ i / A h ∂ i L, ∂ − C ), with the previous notation. Scalar extension.
Let A and B be two rings equippedwith derivations and let ϕ : A → B be a ring homomorphismcommuting with derivation. From a given differential system Y ′ = AY with coefficients in A , one can build a differentialsystem over B by applying ϕ : it is Y ′ = ϕ ( A ) Y , where ϕ ( A )is the matrix obtained from A by applying ϕ entry-wise.This operation admits an analogue at the level of dif-ferential modules: to a differential module ( M, ∂ ) over A ,we attach the differential module ( M B , ∂ B ) over B where M B = B ⊗ ϕ, A M and ∂ B : M B → M B is defined by: ∀ f ∈ B , ∀ x ∈ M, ∂ B ( f ⊗ x ) = f ′ ⊗ x + f ⊗ ∂ ( x ) . It is easily seen that if (
M, ∂ ) is a differential module asso-ciated to the system Y ′ = AY then ( M B , ∂ B ) is that asso-ciated to the system Y ′ = ϕ ( A ) Y . The p -curvature. Let k be any field of characteristic p . Weassume here that A is the field k ( x ) — consisting of rationalfunctions over k — equipped with the standard derivation.he p -curvature of a differential module ( M, ∂ ) over k ( x ) isdefined as the mapping ∂ p : M → M . It follows from theLeibniz rule (4) and the fact that the p -th derivative of anyelement of k ( x ) vanishes that the p -curvature is k ( x )-linear.This definition extends to differential systems as follows:the p -curvature of the system Y ′ = AY is the k ( x )-linearmap ∂ pA : M A → M A where ( M A , ∂ A ) is the correspondingdifferential module. One can check that the matrix of ∂ pA (inthe canonical basis of M A ) is the p -th term of the recursivesequence ( A i ) defined in (1).Considering again a differential operator L and the as-sociated differential module ( A h ∂ i / A h ∂ i L, ∂ − C ), for the as-sociated companion matrix C , we obtain the usual recur-rence A = C and A i +1 = A ′ i + C · A i . The p -curvature of A h ∂ i / A h ∂ i L will simply be called the p -curvature of L .
3. SERIES WITH DIVIDED POWERS
In all this section, we let ℓ be a ring in which p vanishes.We recall the definition of the divided power ring over ℓ ,and its main properties — mainly, a Cauchy–Lipschitz the-orem that will allow us to compute solutions of differen-tial systems. We show how some approaches that are well-known for power series solutions carry over without signifi-cant changes in this context. Most results in this section arenot new; those from § § § ℓ [[ t ]] dp Let ℓ [[ t ]] dp be the ring of formal series of the form: f = a + a γ ( t ) + a γ ( t ) + · · · + a i γ i ( t ) + · · · (5)where the a i ’s are elements of ℓ and each γ i ( t ) is a symbolwhich should be thought of as t i i ! . The multiplication on ℓ [[ t ]] dp is defined by the rule γ i ( t ) · γ j ( t ) = (cid:0) i + ji (cid:1) · γ i + j ( t ). Remark
The ring ℓ [[ t ]] dp is not the PD-envelope inthe sense of [1, 2] of ℓ [[ t ]] with respect to the ideal ( t ) but itscompletion for the topology defined by the divided powers ide-als. Taking the completion is essential to have an analogueof the Cauchy–Lipschitz Theorem ( cf Proposition 3.4).
Invertible elements of ℓ [[ t ]] dp are easily described: they areexactly those for which the“constant” coefficient a is invert-ible in ℓ . The ring ℓ [[ t ]] dp is moreover endowed with a deriva-tion defined by f ′ = P ∞ i =0 a i +1 γ i ( t ) for f = P ∞ i =0 a i γ i ( t ).It then makes sense to consider differential systems over ℓ [[ t ]] dp . A significant difference with power series is theexistence of an integral operator: it maps f as above to R f = P ∞ i =0 a i γ i +1 ( t ) and satisfies ( R f ) ′ = f for all f . Divided power ideals.
For all positive integers N , we de-note by ℓ [[ t ]] dp ≥ N the ideal of ℓ [[ t ]] dp consisting of series of theform P i ≥ N a i γ i ( t ). The quotient ℓ [[ t ]] dp /ℓ [[ t ]] dp ≥ N is a free ℓ -module of rank N and a basis of it is (1 , γ ( t ) , . . . , γ N − ( t )).In particular, for N = 1, the quotient ℓ [[ t ]] dp /ℓ [[ t ]] dp ≥ is iso-morphic to ℓ : in the sequel, we shall denote by f (0) ∈ ℓ thereduction of an element f ∈ ℓ [[ t ]] dp modulo ℓ [[ t ]] dp ≥ . On thewriting (5), it is nothing but the constant coefficient a inthe expansion of f .We draw the reader’s attention to the fact that ℓ [[ t ]] dp ≥ N is not stable under derivation, so the quotients ℓ [[ t ]] dp /ℓ [[ t ]] dp ≥ N do not inherit a derivation. Relationship with ℓ [ t ] . There exists a natural map ε : ℓ [ t ] → ℓ [[ t ]] dp taking a polynomial P i a i t i to P p − i =0 a i i ! · γ i ( t ).The latter sum stops at i = p − i ! becomes divisibleby p after that. Clearly, the kernel of ε is the principal idealgenerated by t p . Hence ε factors through ℓ [ t ] /t p as follows: ℓ [ t ] pr −→ ℓ [ t ] /t p ι −→ ℓ [[ t ]] dp (6)where pr is the canonical projection taking a polynomial toits reduction modulo t p . We observe moreover that the ideal t p ℓ [ t ] is stable under derivation and, consequently, that thequotient ring ℓ [ t ] /t p inherits a derivation. Furthermore, thetwo mappings in (6) commute with the derivation. It turns out that the γ n ( t )’s can all be expressed in termsof only few of them, resulting in a more flexible descriptionof the ring ℓ [[ t ]] dp . To make this precise, we set t i = γ p i ( t )and first observe that t ni = n ! · γ np i ( t ) for all i and n ; this isproved by induction on n , using the equalities t n +1 i = n ! · γ np i ( t ) · γ p i ( t ) = n ! · (cid:0) ( n +1) p i p i (cid:1) γ ( n +1) p i ( t ) , since Lucas’ Theorem shows that (cid:0) ( n +1) p i p i (cid:1) ≡ n +1 (mod p ).In particular t pi = 0 for all i . Proposition
Let n be a positive integer and n = P si =0 n i p i its writing in basis p . Then: γ n ( t ) = γ n ( t ) · γ n p ( t ) · · · γ n s p s ( t ) = t n n ! · t n n ! · · · t n s s n s ! . Proof.
The first equality is proved by induction on s using the fact that if n = a + bp with 0 ≤ a < p , then γ a γ bp = γ n , since (cid:0) a + bpa (cid:1) ≡ p ). The second equalitythen follows from the relations t n i i = n i ! · γ n i p i ( t ).A corollary of the above proposition is that elements of ℓ [[ t ]] dp can be alternatively described as infinite sums ofmonomials a n ,...,n s · t n · t n · · · t n s s where the n i ’s are inte-gers in the range [0 , p ) and the coefficient a n ,...,n s lies in ℓ .The product in ℓ [[ t ]] dp is then the usual product of seriessubject to the additional rules t pi = 0 for all i .More precisely, restricting ourselves to some given preci-sion of the form N = np s , we deduce from the above discus-sion the following corollary. Corollary
For N = np s , with s ∈ N and n ∈{ , . . . , p } , there is a canonical isomorphism of ℓ -algebras: ℓ [[ t ]] dp /ℓ [[ t ]] dp ≥ N ≃ ℓ [ t , . . . , t s ] / ( t p , . . . , t ps − , t ns ) . For instance, if we take s = 0 and N = n in { , . . . , p } , weobtain the isomorphism ℓ [[ t ]] dp /ℓ [[ t ]] dp ≥ N ≃ ℓ [ t ] /t N .In terms of complexity, the change of bases between left-and right-hand sides can both be done in O ˜( N ) operationsin ℓ : all the factorials we need can be computed once and forall for O (min( N, p )) operations; then each monomial con-version takes O ( s ) = O (log( N )) operations, for a total of O ( N log( N )) = O ˜( N ).The previous corollary is useful in order to devise a mul-tiplication algorithm for divided powers, since it reducesthis question to multivariate power series multiplication (ad-dition takes linear time in both bases). To multiply in ℓ [ t , . . . , t s ] / ( t p , . . . , t ps − , t ns ), one can use a direct algorithm:multiply and discard unwanted terms. Using for instanceronecker’s substitution and FFT-based univariate arith-metic, we find that a multiplication in ℓ [[ t ]] dp at precision N ( i.e. modulo ℓ [[ t ]] dp ≥ N ) can be performed with O ˜(2 log p N N )operations in k . A solution that leads to a cost N ε forany ε > A nice feature of the ring ℓ [[ t ]] dp — which does not holdfor ℓ [[ t ]] notably — is the existence of an analogue of theclassical Cauchy–Lipschitz theorem. This property will havea fundamental importance for the purpose of our paper; seefor instance [21, Proposition 4.2] for similar considerations. Proposition
Let Y ′ = AY be a differential systemof dimension r with coefficients in ℓ [[ t ]] dp . For all initialdata V ∈ ℓ r (considered as a column vector) the followingCauchy problem has a unique solution in ℓ [[ t ]] dp : ( Y ′ = A · YY (0) = V. Proof.
Let us write the expansions of A and Y : A = ∞ X i =0 A i γ i ( t ) and Y = ∞ X i =0 Y i γ i ( t )where the A i ’s and Y i ’s have coefficients in ℓ . The Cauchyproblem translates to Y = V and Y n +1 = P ni =0 (cid:0) ni (cid:1) · A i · Y n − i . It is now clear that it has a unique solution.Of course, Proposition 3.4 extends readily to the casewhere the initial data V is any matrix having r rows. Inparticular, taking V = I r (the identity matrix of size r ), wefind that there exists a unique r × r matrix Y with coeffi-cients in ℓ [[ t ]] dp such that Y (0) = I r and Y ′ = A · Y . Thismatrix is often called a fundamental system of solutions . Finding solutions using Newton iteration.
In charac-teristic zero, it is possible to compute power series solutionsof a differential system such as Y ′ = A · Y using Newtoniteration; an algorithm for this is presented on [5, Fig. 1].One can use this algorithm to compute a fundamentalsystem of solutions in our context. For this, we first needto introduce two notations. Given an element f ∈ ℓ [[ t ]] dp written as f = P i a i γ i ( t ) together with an integer m , weset ⌈ f ⌉ m = P m − i =0 a i γ i ( t ). Similarly, if M is a matrix withcoefficients in ℓ [[ t ]] dp , we define ⌈ M ⌉ m and R M by applyingthe corresponding operations entry-wise. Algorithm fundamental_solutions
Input: a differential system Y ′ = AY , an integer N Output: the fund. system of solutions modulo ℓ [[ t ]] dp ≥ N Y = I r + t A (0); Z = I r ; m = 22. while m ≤ N/ Z = Z + (cid:6) Z ( I r − Y Z ) (cid:7) m Y = Y − l Y (cid:0) R Z · ( Y ′ − ⌈ A ⌉ m − Y ) (cid:1)m m m = 2 m return Y Correction is proved as in the classical case [5, Lemma 1].Let us take n ∈ { , . . . , p } and s ∈ N such that n − N written in basis p , and s the corresponding exponent; then, we have ( n − p s ≤ N < np s . Since we areonly interested in costs up to logarithmic factors, we mayassume that we do all operations at precision np s (a betteranalysis would take into account the fact that the precisiongrows quadratically).By Corollary 3.3 and the discussion that follows, arith-metic operations in ℓ [[ t ]] dp /ℓ [[ t ]] dp ≥ np s take time O ˜(2 log p N N ).This is also the case for differentiation and integration, inview of the formulas given in the previous subsection; trun-cation is free. The total complexity of Algorithm funda-mental_solutions is therefore O ˜(2 log p N Nr ω ) operationsin ℓ , where r is the dimension of the differential system. If N = p O (1) , which is what we need later on, this is O ˜( Nr ω ). The case of differential operators.
We now considerthe case of the differential system associated to a differentialoperator L = a r ∂ r + · · · + a ∂ + a ∈ ℓ [[ t ]] dp h ∂ i . We will workunder the following few assumptions: we assume that a r isinvertible, and that there exists an integer d < p such thatall a i ’s can be written a i = α i, + α i, γ ( t ) + · · · + α i,d γ d ( t )for some coefficients α i,j in ℓ ; thus, by assumption, α r, isa unit in ℓ . Our goal is still to compute a basis of solutionsup to precision N ; the algorithm is a direct adaptation of aclassical construction to the case of divided powers.In all that follows, we let f , . . . , f r − be the solutions of L in ℓ [[ t ]] dp , such that f i is the unique solution of the Cauchyproblem ( cf Proposition 3.4): L ( f i ) = 0 ; f ( j ) i (0) = δ ij for 0 ≤ j < r (7)where δ ij is the Kronecker delta. For f = P ∞ j =0 ξ j γ j ( t ) in ℓ [[ t ]] dp , a direct computation shows that the n -th coefficientof L ( f ) is P ri =0 P nj =0 (cid:0) nj (cid:1) α i,j ξ n + i − j . Assume L ( f ) = 0.Then, extracting the term in ξ n + r , and using that α i,j = 0for j > d , we get ξ n + r = − α r, P r − i =0 P dj =0 (cid:0) nj (cid:1) α i,j ξ n + i − j .Letting m = i − j , we find ξ n + r = P r − m = − d A m ( n ) ξ n + m with A m ( n ) = − α r, r − X i =0 (cid:0) ni − m (cid:1) α i,i − m = X ≤ i ≤ r − ≤ i − m ≤ d − α i,i − m α r, ( i − m )! n i − m and n i − m = n ( n − · · · ( n − ( i − m − A m is well-defined, since weassumed that d < p , and shows that A m is a polynomial ofdegree at most d .From this, writing the algorithm is easy. We need two sub-routines: from falling factorial ( F ), which computes theexpansion on the monomial basis of a polynomial of the form F = P ≤ j ≤ n f j n j , and eval ( F, N ), which computes the val-ues of a polynomial F at the N points { , . . . , N − } . Theformer can be done using the divide-and-conquer algorithmof [8, Section 3] in time O ˜( n ); the latter by the algorithmof [17, Chapter 10], in time O ˜(deg( F ) + N ). The previousdiscussion leads to the algorithm solutions_operator be-low. In view of the previous discussion, the cost analysis isstraightforward (at step 2., notice that all required factorialscan be computed in time O ( d )). The costs reported in thepseudo-code indicate the total amount of time spent at thecorresponding line. Algorithm solutions_operator
Input: a differential operator L ∈ ℓ [[ t ]] dp h ∂ i of bidegree( d, r ), with d < p ; an integer N Output: the solutions f , . . . , f r − at precision N . for m = − d, . . . , r − A m = P ≤ i ≤ r − , ≤ i − m ≤ d − α i,i − m α r, ( i − m )! x i − m Cost: O ( d ( r + d ))3. A m = from falling factorial ( ˆ A m ) Cost: O ˜( d ( r + d ))4. Store eval ( A m , N − r ) Cost: O ˜(( d + N )( r + d ))5. for i = 0 , . . . , r − f i = [0 , . . . , , , , . . . ,
0] ( i th unit vector of length r ) Cost: O ( r )7. for n = 0 , . . . , N − r − f i,n + r = P r − m = − d A m ( n ) f i,n + m Cost: O ( rN ( r + d ))9. return f , . . . , f r − Altogether, we obtain the following result, where we usethe assumption
N > d to simplify slightly the cost estimate.
Lemma
Suppose that p < d . Given a positive
N > d ,the classes of f , . . . , f r − modulo ℓ [[ t ]] dp ≥ N can be computedwith at most O ( rN ( r + d )) operations in ℓ . In particular, Algorithm solutions_operator has a bet-ter cost than fundamental_solutions when d = O ( r ω − ).
4. COMPUTING THE P-CURVATURE
In all this section, we work over a field k of characteristic p >
0. We consider a differential system Y ′ = AY of di-mension r and denote by A p the matrix of its p -curvature.We write A = f A ˜ A , where f A is in k [ x ] and ˜ A is a matrixwith polynomial entries. Let d = max(deg f A , deg ˜ A ), wheredeg ˜ A is the maximal degree of the entries of ˜ A . We recall([13, Prop. 3.2], [9, Lemma 1]) a bound on the size of A p .The bound follows from the recurrence (1), and it is tight. Lemma
The entries of the matrix f pA · A p are all poly-nomials of degree at most dp . The goal of this section is to prove the following theorem.
Theorem
There exists an algorithm (presented be-low) which computes the matrix of the p -curvature of thedifferential system Y ′ = AY in O ˜ (cid:0) pdr ω ) operations in k . It is instructive to compare this cost with the size of the out-put. By Lemma 4.1, the latter is an r × r matrix whose en-tries are rational functions whose numerator and denomina-tor have degree ≃ pd , so its size is roughly pdr elements of k .Our result O ˜ (cid:0) pdr ω ) is quasi-optimal if we assume that ma-trix multiplication can be performed in quasi-optimal time. Let A p denote the matrix of the p -curvature of the dif-ferential system Y ′ = AY (in the usual monomial basis).The expression of A p given at the very end of § A p using theframework of divided powers.In order to relate k ( x ) and a ring ℓ [[ t ]] dp , we pick a sepa-rable polynomial S ∈ k [ x ] which is coprime with f A and set ℓ = k [ x ] /S (which is thus not necessarily a field). Let a ∈ ℓ be the class of x . We consider the ring homomorphism: ϕ S : k [ x ] → ℓ [ t ] /t p f ( x ) f ( t + a ) mod t p . Regarding the differential structure, we observe that ϕ S commutes with the derivation when ℓ [ t ] /t p is endowed withthe standard derivation ddt . We furthermore deduce fromthe fact that S and f A are coprime that ϕ S extends to ahomomorphism of differential rings k [ x ][ f A ] → ℓ [ t ] /t p thatwe continue to denote by ϕ S . We set ψ S = ι ◦ ϕ S where ι is the canonical inclusion ℓ [ t ] /t p ֒ → ℓ [[ t ]] dp ( cf § ψ S commutes with the derivation. Finally, because S is separable, we can check that ϕ S is surjective and itskernel is the ideal generated by S p . Hence ϕ S induces anisomorphism: k [ x ] /S p = k [ x ][ f A ] /S p ∼ −→ ℓ [ t ] /t p . (8)Let Y S be a fundamental system of solutions of the dif-ferential system Y ′ = ψ S ( A ) · Y , i.e. Y S is an r × r ma-trix with coefficients in ℓ [[ t ]] dp such that Y S (0) = I r and Y ′ S = ψ S ( A ) · Y S . The existence of Y S is guaranteed byProposition 3.4. Moreover, the matrix Y S is invertible be-cause Y S (0) = I r is. Proposition
Keeping the above notations, we have: ϕ S ( A p ) = − Y ( p ) S · Y − S (9) where Y ( p ) S is the matrix obtained from Y S by taking the p -thderivative entry-wise. Proof.
We set Z S = Y − S and let ( M, ∂ ) denote thedifferential module over ℓ [[ t ]] dp associated to the differen-tial system Y ′ = ψ S ( A ) Y . Let y , . . . , y r denote the col-umn vectors of Y S . They are all solutions of the system Y ′ = ψ S ( A ) Y , meaning that ∂ ( y i ) = 0 for all i . Further-more, if ( e , . . . , e r ) is the canonical basis of ( ℓ [[ t ]] dp ) r , wehave the matrix relations: t Y S · e = y and e = t Z S · y where y (resp. e ) is the column vector whose coordinates are the vectors y i ’s (resp. the e i ’s). Applying ∂ to the above rela-tion, we find ∂ ( e ) = t Z ′ S · y + t Z S · ∂ ( y ) = t Z ′ S · y and iteratingthis p times, we deduce ∂ p ( e ) = t Z ( p ) S · y = t Z ( p ) S · t Y S · e .On the other hand, the matrix ψ S ( A p ) of the p -curvature isdefined by the relation ∂ p ( e ) = t ψ S ( A p ) · e . Therefore we get ψ S ( A p ) = Y S · Z ( p ) S . Now differentiating p times the relation Y S Z S = I r , we find Y ( p ) S Z S + Y S · Z ( p ) S = 0. Combining thiswith the above formula for ψ S ( A p ) concludes the proof.In our setting, the matrix A p has coefficients in k [ x ][ f A ]( cf Lemma 4.1), from which we deduce that ψ S ( A p ) hasactually coefficients in the subring ℓ [ t ] /t p of ℓ [[ t ]] dp . There-fore, using Eq. (9), one can compute ψ S ( A p ) knowing only Y S modulo the ideal ℓ [[ t ]] dp ≥ p .One can actually go further in this direction and establisha variant of Eq. (9) giving an expression of ψ S ( A p ) whichinvolves only the reduction of Y S modulo ℓ [[ t ]] dp ≥ p . To makethis precise, we need an extra notation. Given an integer i ∈ [0 , p ) and a polynomial f ∈ ℓ [ t ] /t p (resp. a matrix M with coefficients in ℓ [ t ] /t p ), we write Coeff( f, i ) (resp.Coeff( M, i )) for the coefficient in t i in f (resp. in M ). roposition Keeping the above notations, we have: ψ S ( A p ) = − ¯ Y S · Y ( p ) S (0) · ¯ Y − S = ¯ Y S · Coeff( A · ¯ Y S , p − · ¯ Y − S (10) where we have set ¯ Y S = Y S mod ℓ [[ t ]] dp ≥ p . Proof.
Differentiating p times the relation Y ′ S = ψ S ( A ) · Y S , we observe that Y ( p ) S is solution of the same differ-ential system Y ′ = ψ S ( A ) Y . Hence, thanks to unique-ness in Cauchy–Lipschitz Theorem, we have the relation Y ( p ) S = Y S · Y ( p ) S (0). The first part of the Proposition followsby plugging this in Eq. (9) and reducing the result modulo ℓ [[ t ]] dp ≥ p . To establish the second part, it is now enough tonotice that the relation Y ′ S = ψ S ( A ) · Y S implies: Y ( p ) S (0) = ( A · Y S ) ( p − (0) = − Coeff( A · ¯ Y S , p − p − ≡ − p ). Remark
We can rephrase Proposition 4.4 as fol-lows: letting y , . . . , y r denote the column vectors of Y S and ¯ y i ∈ ( ℓ [ t ] /t p ) r be the reduction of y i , the p -curvature of A modulo t p is the linear endomorphism of ( ℓ [ t ] /t p ) r whosematrix in the basis (¯ y , . . . , ¯ y r ) is Coeff( A · ¯ Y S , p − . It isworth remarking that the latter matrix has coefficients in thesubring ℓ of ℓ [ t ] /t p . Remembering Eq. (8), we conclude that Proposition 4.4allows us to compute the image of the p -curvature A p mod-ulo S p . The strategy of our algorithm now becomes clear:we first compute A p modulo S p for various polynomials S and, when we have collected enough congruences, we putthem together to reconstruct A p . The first step is detailedin § § In all this subsection, we fix a separable polynomial S ∈ k [ x ] and denote by m its degree. Our goal is to design analgorithm for computing the matrix A p modulo S p . AfterProposition 4.4, the main remaining algorithmic issue is theeffective computation of the isomorphism ϕ S and its inverse. Applying ϕ S and its inverse. We remark that ϕ S factorsas follows: k [ x ] /S p → k [ x, t ] / h S, ( t − x ) p i → k [ x, t ] / h S, t p i x t t + a. Applying the right-hand mapping, or its inverse, amountsto doing a polynomial shift in degree p with coefficients in k [ x ] /S . Using the divide-and-conquer algorithm of [16], thiscan be done in O ˜( p ) arithmetic operations in k [ x ] /S , whichis O ˜( pm ) operations in k . Thus, we are left with the left-hand factor, say ϕ ⋆S . Applying it is straightforward and canbe achieved in O ˜( pm ) operations in k . It then only remainsto explain how one can apply efficiently ϕ ⋆S − .We start by determining the image of x by ϕ ⋆S − ; call it y = ϕ ⋆S − ( x ); we may identify it with its canonical preim-age in k [ x ], which has degree less than pm . Write y = P ≤ i
1. The first equality implies that x p generates k [ x ] /S , so the fact that ζ hasdegree less than m implies that ζ is the unique polynomialwith this degree constraint such that ζ ( x p ) = x mod S . Theother equalities then imply that ζ i = 0 for i = 1 , . . . , p − ζ , we first compute ν = x p mod S , using O ˜( m log( p )) operations in k . Then, we have tofind the unique polynomial ζ of degree less than m suchthat ζ ( ν ) = x mod S . In general, one can compute ζ in O ( m ω ) operations in k by solving a linear system. In thecommon case where m < p , there exists a better solution.Indeed, denote by tr : k [ x ] /S → k the k -linear trace formand write t i = tr( ν i ) and t ′ i = tr( xν i ), for i = 0 , . . . , m −
1. Then formulas such as those in [25] allow us to recover ζ from t = ( t , . . . , t m − ) and t ′ = ( t ′ , . . . , t ′ m − ) in time O ˜( m ). These formulas require that m < p and that S ′ be invertible modulo S , which is ensured by our assumptionthat S is separable. To compute t and t ′ , we can use Shoup’spower projection algorithm [29], which takes O ( m ( ω +1) / )operations in k .Once ζ is known, to apply the mapping ϕ ⋆S − to an el-ement g ( x, t ), we proceed coefficient-wise in t . Write g = P ≤ i
In the case where S = x m − c (where c ∈ k and p does not divide m ), there actually exists a quite simpleexplicit formula for ϕ ⋆S − : it takes t to x and x to c q x pn where n and q are integers satisfying the B´ezout’s relation pn + qm = 1 . Using this, one can compute ϕ ⋆S − ( g ) in O ˜( pm ) operations in k in this special case. Conclusion.
Let us call phiS and phiS_inverse the twosubroutines described above for computing ϕ S and its in-verse respectively. Proposition 4.4 leads to the followingalgorithm for computing the p -curvature modulo S p . Algorithm local_p_curvature
Input: a polynomial S and a matrix A S ∈ M r ( k [ x ] /S p ) Output: the p -curvature of the system Y ′ = A S Y A S,ℓ = phiS ( A S ) Cost: O ˜( pr m ) operations in k (with m = deg S )2. compute a fund. system of solutions Y S ∈ M r ( ℓ [ t ] /t p )of the system Y ′ = A S,ℓ Y at precision p . Cost: O ˜( pr ω ) op. in ℓ using fundamental_solutions Remark:
Here ℓ = k [ x ] /S A p,ℓ = Y S · Coeff( AY S , p − · Y − S at precision O ( t p ) Cost: O ˜( pr ω ) operations in ℓ A p = phiS_inverse ( A p,ℓ ) Cost: O ˜( pr m ω ) operations in k in general O ˜( pr m ( ω +1) / ) operations in k if m < p return A p .o conclude with, it is worth remarking that implementingthe algorithm local_p_curvature can be done using usualpower series arithmetic: indeed, we only need to performcomputations in the quotient ℓ [[ t ]] dp /ℓ [[ t ]] dp ≥ p which is iso-morphic to ℓ [ t ] /t p by Corollary 3.3. Furthermore, we notethat if we are using the algorithm fundamental_solutions at line 2, then Y − S can be computed by performing an ex-tra loop in fundamental_solutions ; indeed the matrix Z we obtain this way is exactly Y − S . We recall that we have started with a differential system Y ′ = AY (with A = f A ˜ A ) and that our goal is to computethe matrix A p of its p -curvature. Lemma 4.1 gives boundson the size of the entries of A p . We need another lemma,which ensures that we can find enough small “evaluationpoints” (lying in a finite extension of k ). Let F p denote theprime subfield of k . Lemma
Given a positive integer D and a nonzeropolynomial f ∈ k [ x ] , there exist pairwise coprime polynomi-als S , . . . , S n ∈ F p [ x ] with n ≤ D such that: • P ni =1 deg S i ≥ D • for all i , the polynomial S i is coprime with f and hasdegree at most p ( D + deg f ) . Proof.
Let m be the smallest integer such that p m ≥ D +deg f . Clearly m ≤ q ( D +deg f ) ≤ p ( D +deg f ).Let F p m be an extension of F p of degree m and K be thecompositum of k and F p m . Let S , . . . , S t be the mini-mal polynomials over F p (without repetition) of all elementsin F p m ⊂ K which are not a root of f . We then havedeg S i ≤ m for all i and P ti =1 deg S i ≥ p m − deg f ≥ D .It remains now to define n as the smallest integer such that P ni =1 deg S i ≥ D . Minimality implies P n − i =1 deg S i < D and thus n ≤ D . Therefore S , . . . , S n satisfy all the re-quirements of the lemma.The above proof yields a concrete algorithm for producinga sequence S , . . . , S n satisfying the properties of Lemma4.7: we run over elements in F p m and, for each new element,append its minimal polynomial over F p to the sequence ( S i )unless it is not coprime with f . We continue this processuntil the condition P ni =1 deg S i ≥ D holds. Keeping in mindthe logarithmic bound on m , we find that the complexity ofthis algorithm is at most O ˜( D + deg f ) operations in k .Let us call generate_points the resulting routine: it takesas input the parameters f and D and return an admissiblesequence S , . . . , S n .We are now ready to present our algorithm for computingthe p -curvature: Algorithm p_curvature
Input: a matrix A written as A = f A · ˜ A Output: the p -curvature of the differential system Y ′ = AY S , . . . , S n = generate_points ( f A , d + 1) Cost: O ˜( d ) operations in k Remark: we have n = O ( d ) and deg S i = O (log d ), ∀ i for i = 1 , . . . , n : A i,p = local_p_curvature ( S i , A mod S pi ) Cost: O ˜( pdr ω ) operations in k
3. compute B ∈ M r ( k [ x ]) with entries of degree ≤ pd such that B ≡ f pA · B i (mod S pi ) for all i Cost: O ˜( pdr ) operations in k return f pA · B In view of the previous discussion and Lemma 4.1, the cor-rectness and the cost analysis of the algorithm p_curvature are both straightforward. Hence, Theorem 4.2 is proved.We conclude this subsection with three remarks. First,when applying Chinese Remainder Theorem (CRT) on line3 of Algorithm p_curvature , we notice that all moduli S pi are polynomials in x p . This allows the following optimiza-tion. Writing f pA · B i ≡ P p − j =0 B i,j ( x p ) x j (mod S pi ( x )) anddenoting by C j the unique solution of degree at most d tothe congruence system: B j ( x ) ≡ B i,j ( x ) (mod T i ( x )) where T i ( x p ) = S pi ( x )we have B = P p − j =0 B j ( x ) x j . This basically allows us to re-place one CRT with polynomials of degree dp by p CRT withpolynomials of degree d . We save this way the polynomialfactors in log( p ) in the complexity.Second, instead of working with n polynomials S i , onemay alternatively choose a unique polynomial S of the form S = X m − a where m ≥ d is an integer not divisible by p and a ∈ k are such that S and f A are coprime. This avoidsthe use of Chinese Remainder Theorem and the resultingcomplexity stays in O ˜( pdr ω ) provided that we use Remark4.6 in order to compute the inverse of ϕ S .Third, we observe that the algorithm p_curvature is veryeasily parallelizable. Indeed, each iteration of the main loop(on line 2) is completely independent from the others. Thus,they all can be performed in parallel. Moreover, accordingto the first remark (just above), the application of the Chi-nese Remainder Theorem (on line 3) splits into pr smallerindependent problems and can therefore be efficiently par-allelized as well. To conclude with, we would like to discuss the case of adifferential operator L = a r ∂ r + a r − ∂ r − + · · · + a ∂ + a with a i ∈ k [ x ] for all i , of maximal degree d .Recall that the p -curvature of L is that of the differen-tial module ( A h ∂ i / A h ∂ i L, ∂ − C ), where C is the compan-ion matrix associated to L as in (3). Applying directlythe formulas in Proposition 4.4 requires the knowledge ofthe solutions of the system Y ′ = − CY . It is in fact eas-ier to compute solutions for the system X ′ = t CX , sincewe saw that these solutions are the vectors of the form t ( y, y ′ , . . . , y ( r − ), where y is a solution of L . This is how-ever harmless: the p -curvatures A p and B p of the respectivesystems Y ′ = − CY and X ′ = t CX (which are so-called ad-joint) satisfy A p = − t B p . Thus, we can use the formulasgiven above to compute ϕ S ( B p ), and deduce ϕ S ( A p ) for anegligible cost. Equivalently, one may notice that the funda-mental matrices of solutions of our two systems are transposeof one another, up to sign.Moreover, instead of using the second formula of Propo-sition 4.4 to compute the local p -curvatures, we recommendusing the first one, which is ϕ S ( B p ) = − X S · X ( p ) S (0) · X − S where X S is a fundamental system of solutions of X ′ = t CX and ¯ X S denotes its reduction in M r ( ℓ [ t ] /t p ). If f , . . . , f r − are solutions of the system (7), the ( i, j )-th entry of X S is just f ( i ) j . Hence the matrices ¯ X S and X ( p ) S (0) can beobtained from the knowledge of the image of f i ’s modulo
157 281 521 983 1 811 3 433 6 421 12 007 d = 5, r = 5 0 .
39 s 0 .
71 s 1 .
22 s 2 .
34 s 4 .
41 s 8 .
93 s 18 . . .
26 s 0 .
76 s 2 .
69 s 9 .
05 s 32 . d = 5, r = 11 1 .
09 s 2 .
05 s 3 .
65 s 7 .
05 s 12 . . . .
25 s 3 .
70 s 12 . . − d = 5, r = 20 2 .
93 s 5 .
25 s 9 .
52 s 17 . . . .
29 s 12 . . − − d = 11, r = 20 6 .
89 s 13 . . . . . . − − − d = 20, r = 20 14 . . . . . . − − − − Running times obtained with Magma V2.19-4 on an AMD Opteron 6272 machine at 2GHz and 8GB RAM, running Linux.
Figure 1: Average running time on random inputs of various sizes ℓ [[ t ]] dp ≥ p + r just by reorganizing coefficients (and possibly mul-tiplying by some factorials depending on the representationof elements of ℓ [[ t ]] dp we are using).As for the f i ’s, they can be computed by the algorithm so-lutions_operator (provided its assumptions are satisfied).We need finally to compute X − S : since X S (0) is the iden-tity matrix, this can be done either using Newton iterator,a divide-and-conquer approach or a combination of both,which computes the inverse of X S at a small precision, anduses divide-and-conquer techniques for higher ones (the lat-ter being the most efficient in practice). All these remarksdo speed up the execution of our algorithms when d is nottoo large compared to r .Last but not least, we notice that, in the case of differ-ential operators, the matrix A p is easily deduced from itsfirst column. Indeed, writing A p = ( a i,j ) ≤ i,j 5. IMPLEMENTATION AND TIMINGS We implemented our algorithms in Magma in the case ofdifferential operators; the source code is available at https://github.com/schost .Figure 1 gives running times for random operators of degrees( d, r ) in k [ x ] h ∂ i and compares them with running times of(a fraction free version of) Katz’s algorithm which consistsin computing the recursive sequence ( A i ) until i = p . Ineach cell, the first line (resp. the second line) corresponds tothe running time obtained with our algorithm (resp. Katz’salgorithm); a dash indicates that the corresponding runningtime exceeded one hour. Our benchmarks rather well reflectthe predicted dependence with respect to p : quasi-linear forour algorithm and quadratic for Katz’s algorithm.Larger examples (than those presented in Fig. 1) are alsoreachable: for instance, we computed the first column of the p -curvature of a “small” multiple of the operator φ (5) H con-sidered in [10, Appendix B.3] modulo the prime 27449. Thisoperator has bidegree ( d, r ) = (108 , 6. REFERENCES [1] P. Berthelot. Cohomologie cristalline des sch´emas de caract´eristique p > Notes on crystalline cohomology . PrincetonUniversity Press, Princeton, N.J.; University of Tokyo Press, Tokyo, 1978.[3] 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.[4] 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.[5] A. Bostan, F. Chyzak, F. Ollivier, B. Salvy, ´E. Schost, and A. Sedoglavic.Fast computation of power series solutions of systems of differentialequations. In , pages1012–1021, 2007. New Orleans, January 2007.[6] A. Bostan and M. Kauers. Automatic classification of restricted latticewalks. In FPSAC’09 , DMTCS Proc., AK, pages 201–215. 2009.[7] 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.[8] A. Bostan and ´E. Schost. Polynomial evaluation and interpolation onspecial sets of points. J. Complexity , 21(4):420–446, 2005.[9] A. Bostan and ´E. Schost. Fast algorithms for differential equations inpositive characteristic. In ISSAC’09 , pages 47–54. ACM, New York, 2009.[10] S. Boukraa, S. Hassani, J.-M. Maillard, and N. Zenine. Singularities of n -fold integrals of the Ising class and the theory of elliptic curves. J. Phys.A , 40(39):11713–11748, 2007.[11] D. G. Cantor and E. Kaltofen. On fast multiplication of polynomials overarbitrary algebras. Acta Inform. , 28(7):693–701, 1991.[12] A. Chambert-Loir. Th´eor`emes d’alg´ebricit´e en g´eom´etrie diophantienne(d’apr`es J.-B. Bost, Y. Andr´e, D. & G. Chudnovsky). S´eminaire Bourbaki ,282(886):175–209, 2002.[13] T. Cluzeau. Factorization of differential systems in characteristic p . In ISSAC’03 , pages 58–65. ACM Press, 2003.[14] T. Cluzeau and M. van Hoeij. A modular algorithm for computing theexponential solutions of a linear differential operator. J. Symbolic Comput. ,38(3):1043–1076, 2004.[15] F. L. Gall. Powers of tensors and fast matrix multiplication. In ISSAC’14 ,pages 296–303, 2014.[16] J. von zur Gathen and J. Gerhard. Fast algorithms for Taylor shifts andcertain difference equations. In ISSAC’97 , pages 40–47. ACM, 1997.[17] J. von zur Gathen and J. Gerhard. Modern Computer Algebra . CambridgeUniversity Press, Cambridge, second edition, 2003.[18] D. Harvey, J. van der Hoeven, and G. Lecerf. Faster polynomialmultiplication over finite fields. http://arxiv.org/abs/1407.3361 , 2014.[19] N. M. Katz. Algebraic solutions of differential equations ( p -curvature andthe Hodge filtration). Invent. Math. , 18:1–118, 1972.[20] N. M. Katz. A conjecture in the arithmetic theory of differentialequations. Bull. Soc. Math. France , (110):203–239, 1982.[21] W. F. Keigher and F. L. Pritchard. Hurwitz series as formal functions. J.Pure Appl. Algebra , 146(3):291–304, 2000.[22] M. van der Put. Differential equations in characteristic p . CompositioMathematica , 97:227–251, 1995.[23] M. van der Put. Reduction modulo p of differential equations. Indag.Mathem. , 7(3):367–387, 1996.[24] M. van der Put and M. Singer. Galois theory of linear differential equations .Springer, 2003.[25] F. Rouillier. Solving zero-dimensional systems through the rationalunivariate representation. Appl. Algebra Engrg. Comm. Comput. , 9(5):433–461,1999.[26] A. Sch¨onhage. Schnelle Multiplikation von Polynomen ¨uber K¨orpern derCharakteristik 2. Acta Informatica , 7:395–398, 1977.[27] A. Sch¨onhage and V. Strassen. Schnelle Multiplikation großer Zahlen. Computing , 7:281–292, 1971.[28] ´E. Schost. Multivariate power series multiplication. In ISSAC’05 , pages293–300. ACM, 2005.[29] V. Shoup. Fast construction of irreducible polynomials over finite fields. Journal of Symbolic Computation , 17(5):371–391, 1994.[30] Y. Tang. Algebraic solutions of differential equations over the projectiveline minus three points. http://arxiv.org/abs/1412.7875http://arxiv.org/abs/1412.7875