Power series solutions of (q)-differential equations

Truncated power series are a fundamental class of objects of computer algebra. Fast algorithms are known for a large number of operations starting from addition, derivative, integral and product and extending to quotient, powering and several more. The main open problem is composition: given two power series f and g, with g(0) = 0, known mod x N, the best known algorithm computing f(g) mod x N has a cost which is roughly that of N √ products in precision N; it is not known whether quasi-linear (i.e., linear up to logarithmic factors) complexity is possible in general. Better results are known over finite fields [Ber98, KU11] or when more information on f or g is available. Quasi-linear complexity has been reached when g is a polynomial [BK78], an algebraic series [Hoe02], or belongs to a large class containing for instance the expansions of exp (x) − 1 and log (1 + x) [BSS08]. One motivation for this work is to deal with the case when f is the solution of a given differential equation. Using the chain rule, a differential equation for f(g) can be derived, with coefficients that are power series. We focus on the case when this equation is linear, since in many cases linearization is possible [BCO+07]. When the order n of the equation is larger than 1, we use the classical technique of converting it into a first-order equation over vectors, so we consider equations of the form x k δ(F) = A F + C, (4.1) where A is an n × n matrix over the power series ring k[[x]] (k being the field of coefficients), C and the unknown F are size n vectors over k[[x]] and for the moment δ denotes the differential operator d/d x. The exponent k in (4.1) is a nonnegative integer that plays a key role for this equation.By solving such equations, we mean computing a vector F of power series such that (4.1) holds modulo x N. For this, we need only to compute F polynomial of degree less than N + 1 (when k = 0) or N (otherwise). Conversely, when (4.1) has a power series solution, its first N coefficients can be computed by solving (4.1) modulo x N (when k  0) or x N −1 (otherwise). If k = 0 and the field k has characteristic 0, then a formal Cauchy theorem holds and (4.1) has a unique vector of power series solution for any given initial condition. In this situation, algorithms are known that compute the first N coefficients of the solution in quasi-linear complexity [BCO+07]. Also the relaxed algorithm of [Hoe02] applies to this case. In this article, we extend the results of [BCO+07] in three directions: Singularities We deal with the case when k is positive. A typical example is the computation of the composition F = f(g) when f is Gauss’ 2F1 hypergeometric series. Although f is a very nice power series.

Positive characteristic 

Even when k = 0, Cauchy’s theorem does not hold in positive characteristic and Equation (4.1) may fail to have a power series solution (a simple example is F ′ = F). However, such an equation may have a solution modulo x N. Efficient algorithms for finding such a solution are useful in conjunction with the Chinese remainder theorem. Other motivations for considering algorithms that work in positive characteristic come from applications in number-theory based cryptology or in combinatorics [BMSS08, BSS08, BS09]. Our objectives in this respect are to overcome the lack of a Cauchy theorem, or of a formal theory of singular equations, by giving conditions that ensure the existence of solutions at the required precisions. More could probably be said regarding the p-adic properties of solutions of such equations (as in [BGVPS05, LS08]), but this is not the purpose of this chapter.

Notation and complexity model

We adopt the convention that uppercase letters denote matrices or vectors while lowercase letters denote scalars. The set of n × m matrices over a ring R is denoted Mn,m(R); when n=m, we write Mn(R). If f is in k[[x]], its degree i coefficient is written fi ; this carries over to matrices. The identity matrix is written Id (the size will be obvious from the context). To avoid any confusion, the entry (i, j) of a matrix M is denoted M(i,j) . Our algorithms are sometimes stated with input in k[[x]], but it is to be understood that we are given only truncations of A and C and only their first N coefficients will be used. The costs of our algorithms are measured by the number of arithmetic operations in k they use. We let M:N→N be such that for any ring R, polynomials of degree less than n in R[x] can be multiplied in M(n) arithmetic operations in R. We assume that M(n) satisfies the usual assumptions of [GG03, §8.3]; using Fast Fourier Transform, M(n) can be taken in O(n log (n) loglog (n)) [CK91, SS71]. We note ω ∈ (2, 3] a constant such that two matrices in Mn(R) can be multiplied in O(n ω ) arithmetic operations in R. The current best bound is ω < 2.3727 ([VW11] following [CW90, Sto10]). Our algorithms rely on linear algebra techniques; in particular, we have to solve several systems of non-homogeneous linear equations. For U in Mn(k) and V in Mn,1(k), we denote by LinSolve(U X = V ) a procedure that returns ⊥ if there is no solution, or a pair F , K, where F is in Mn,1(k) and satisfies U F = V , and K ∈Mn,t(k), for some t ≤n, generates the nullspace of U. This can be done in time O(n ω ). In the pseudo-code, we adopt the convention that if a subroutine returns ⊥, the caller returns ⊥ too (so we do not explicitly handle this as a special case).

Main results

Equation (4.3) is linear, non-homogeneous in the coefficients of F, so our output follows the convention mentioned above. We call generators of the solution space of Eq. (4.3) at precision N either ⊥ (if no solution exists) or a pair F , K where F ∈ Mn,1(k[x]) and K ∈ Mn,t(k[x]) with t ≤ n N, such that for G∈Mn,1(k[x]), with deg (G)< N, x k δ(G) =A σ(G) +C modx N if and only if G can be written G = F + K B for some B ∈Mt,1(k). Seeing Eq. (4.3) as a linear system, one can obtain such an output using linear algebra in dimension n N. While this solution always works, we give algorithms of much better complexity, under some assumptions related to the spectrum SpecA0 of the constant coefficient A0 of A. First, we simplify our problem: we consider the case k = 0 as a special case of the case k = 1. Indeed, the equation δ(F) = A σ(F) +C modx N is equivalent to x δ(F)=P σ(F) + Q modx N +1, with P =xA and Q = xC. Thus, in our results, we only distinguish the cases k = 1 and k > 1.