An algorithm for computing invariants of differential Galois groups

Mark van Hoeij

Universiteit Nijmegen, hoeij@sci.kun.nl

Jacques-Arthur Weil^{1}

Rijksuniversiteit Groningen, J.A.Weil@math.rug.nl

This paper presents an algorithm to compute invariants of the differential Galois group of linear differential equations L(y) = 0: if V (L) is the vector space of solutions of L(y) = 0, we show how given some integer m, one can compute the elements of the symmetric power Sym^{m}(V (L)) that are left fixed by the Galois group. The bottleneck of previous methods is the construction of a differential operator called the ‘symmetric power of L’. Our strategy is to split the work into first a fast heuristic that produces a space that contains all invariants, and second a criterion to select all candidates that are really invariants. The heuristic is built by generalizing the notion of exponents. The checking criterion is obtained by converting candidate invariants to candidate dual first integrals; this conversion is done efficiently by using a symmetric power of a formal solution matrix and showing how one can reduce significantly the number of entries of this matrix that need to be evaluated.

⋆ Expanded version of a talk presented at the MEGA 96 conference (Eindhoven, June 4-8 1996). The work in this paper was prepared in the Dutch NWO-project AIDA. 1 Research supported by the Dutch NWO-MRI (project AIDA). The authors thank the universities of Groningen and Nijmegen as well as the École Polytechnique for their support and hospitality during the preparation of this paper. Some computations were done on the machines of the CNRS-GDR MEDICIS at the École Polytechnique.

Preprint submitted to Elsevier Preprint 1 October 1996

1 Introduction

Let C be a field of characteristic 0 and \overline{C} be its algebraic closure. Denote k = \overline{C}(x) with the derivation ∂ = \frac{d}{dx}. Let

L(y) = ∑_{i=0}^{n} a_{i}y^{(i)} = 0,  a_{n} ≠ 0,  a_{i} ∈ C[x]

denote a homogeneous linear differential equation of n-th order. For such differential equations, there is a differential Galois theory analogous to that for polynomial equations. Let y_{1}, …, y_{n} be a basis of the vector space V (L) of solutions. By adjoining the solutions y_{1}, …, y_{n} and all their derivatives to k, we get a differential field extension K ⊃ k (called a Picard-Vessiot extension); the differential Galois group G of L (over k) is then defined as the group of k-automorphisms of the differential field K (i.e. k-automorphisms of K that commute with the derivation). The group G acts faithfully on the vector space V (L), and so G can be viewed as a subgroup of GL(V (L)). It ‘measures’ the differential relations satisfied by the solutions of L(y) = 0 over k. One way to obtain information on G (and thus on the solutions) is to compute invariants:

Definition 1.1 An element v of some symmetric power Sym^{m}(V (L)) that is fixed by the differential Galois group G is called an invariant of G.

A standard method for computing invariants consists in building an operator L^{◯ s m} (for a definition see section 2.3) whose solution space is a G-homomorphic image of Sym^{m}(V (L)) and then search for rational solutions of this operator. Via differential Galois theory, one can (usually) reconstruct invariants from these rational solutions (see [SUl96] for more details). However, the computation of L^{◯ s m} can be complicated for computers. For this and for other reasons (cf. section 2.3) we will use the companion system of L, which we note Y′ = AY. It is then easy to construct a system Y′ = S^{m}(A)Y whose solution space is G-isomorphic with Sym(V (L)). Our algorithm consists in finding rational solutions of the latter system under two guidelines: we do not perform a costly conversion into an equation and for efficiency we use as much as possible the structure of the system (i.e. the fact that it is a symmetric power of a companion matrix system).

In section 2, we develop and motivate this approach and its links with the previous methods. Let F denote a rational (i.e. entries in k = \overline{C}(x)) solution of Y′ = S^{m}(A)Y. Such an F will be called a dual first integral. In section 3, we define generalized exponents of a local differential operator and show how to use these to compute bounds for the numerators of denominators of the entries of F. Let Sym^{m}(Û) denote the m-th symmetric power matrix (definitions follow later) of Û, where Û is a fundamental solution matrix of Y′ = AY. Note that Sym^{m}(Û) can be computed from a basis ŷ_{1}, …, ŷ_{n} of formal solutions of L. Using the bounds from section 3, we show in section 4 how the evaluation of a finite number of terms of the series in Sym^{m}(Û) (plus linear algebra) yields all rational solutions of Y′ = S^{m}(A)Y. Our strategy is first to design a fast heuristic to construct a space that contains all invariants (plus maybe some additional rubbish), then to convert these (candidate) invariants to (candidate) dual first integrals using the matrix Sym^{m}(Û). Then we check which candidate invariants are really invariants by checking which candidate dual first integrals are indeed dual first integrals. To do this conversion efficiently, we show how to reduce significantly the number of rows and columns of Sym^{m}(Û) that need to be evaluated. The algorithm is implemented in maple; an experimental code is available from the authors.

Acknowledgement

We would like to thank Elie Compoint, Marius van der Put, Felix Ulmer, and Michael Singer for useful remarks concerning this paper. We also thank Elie for swamping us with complicated examples that helped us with the implementation.

2 Invariants of differential Galois groups

In this section, we recall some basic facts and notations about the various ways to present the invariants of differential Galois groups. For more detailed introductions to differential Galois theory, unfamiliar readers could consult [Sin89,Lev89,MRa89,Ber86,Beu92].

2.1 Two presentations of the invariants

If y is a generic solution of L(y) = 0, we can form the vector

Y = (y, y′, y′′, …, y^{(n−1)})^{t} .

This vector satisfies a first order linear system Y′ = AY, where A is the companion matrix of L. Let y_{1}, …, y_{n} denote a basis, fixed once for all, of the solution space V (L) of L(y) = 0. Then the vectors Y_{i} = (y_{i}, y′_{i}, y′′_{i}, …, y_{i}^{(n−1)})^{t} form a fundamental system of solutions of Y′ = AY. This solution space is G-isomorphic with V (L). The n × n matrix U whose columns are the Y_{i} is called a fundamental solution matrix for Y′ = AY.

2.1.1 Polynomial invariants

It is well known ([Lan92]) that Sym(V (L)) can be identified with the polynomial ring \overline{C}[X_{1}, …, X_{n}], where X_{1}, …, X_{n} are variables on which G acts the same as on y_{1}, …, y_{n}. Under this identification, we will say that a homogeneous polynomial P that is fixed by the Galois group is a polynomial invariant. In the sequel, the coefficients (in \overline{C}) of such a P will be referred to as the vector of coefficients of the invariant, or the vector invariant. Let f = P (y_{1}, …, y_{n}) ∈ K. As P is an invariant, f is fixed by G. The differential Galois correspondence then implies that f ∈ k. We will call this f the value of P.

For an invariant in Sym^{m}(V (L)), the expression of P depends on the choice of the basis of V (L). But the value f of the invariant is independent of this choice. For some applications, one just needs this value (for example to compute algebraic solutions [SUl93b] or to solve second order equations [UWe94]), and there, ‘to compute an invariant’ means ‘to compute its value’. For other applications (to compute Liouvillian solutions [SUl96]), one needs the expression of the polynomial invariant, together with its value.

2.1.2 The symmetric power system

Let y denote again a generic solution of L(y) = 0, and let μ_{[m_{1},...,m_{n}]} := y^{m_{1}} · (y′)^{m_{2}} ⋯ (y^{(n−1)})^{m_{n}} (with ∑ m_{i} = m) denote a differential monomial of degree m in y. Then

μ′_{[m_{1},...,m_{n}]} = m_{1}μ_{[m_{1}−1,m_{2}+1,...,m_{n}]} + ⋯ + m_{n−1}μ_{[m_{1},...,m_{n−1}−1,m_{n}+1]}

+m_{n} (− ∑_{j=1}^{n} \frac{a_{j−1}}{a_{n}} μ_{[...,m_{j}+1,...,m_{n}−1]})

(with the convention that μ_{[...,−1,...]} = μ_{[...,m+1,...]} = 0) so this derivative is a k-linear combination of monomials of degree m in the y^{(i)}. As there are N = \binom{n+m−1}{n−1} such monomials, the vector

Y = (y^{m}, …, y^{(n−2)}(y^{(n−1)})^{m−1}, (y^{(n−1)})^{m})

of all such monomials satisfies an N × N system Y′ = S^{m}(A)Y. Note that the matrix S^{m}(A) is very sparse and that it is immediately given by the relations above. So it can be computed quickly, even for large n and m.

2.1.3 The symmetric power matrix

Let K be a field. The action of g ∈ Gl_{n}(K) on K^{n} induces an action, denoted by Sym^{m}(g), on the vector space Sym^{m}(K^{n}). In other words we have a group homomorphism Sym^{m} : Gl_{n}(K) → Gl(Sym^{m}(K^{n})). After having chosen a ordering on the monomials in X_{1}, …, X_{n} of degree m, we can identify the vector space Sym^{m}(K^{n}) with K^{N} (here N is the number of such monomials; N = \binom{n+m−1}{n−1}). This way a group homomorphism

Sym^{m} : Gl_{n}(K) → Gl_{N}(K). (1)

has been defined.

Remark 2.1 The above definition of Sym^{m}(g) (which from now on will be considered as an element Gl_{N}(K) instead of Gl(Sym^{m}(K^{n}))) depends on the ordering that was chosen for the monomials of degree m. It is irrelevant which ordering we choose, however, to have a consistent definition we must always use the same ordering. We will use the lexicographic ordering with X_{1} < ⋯ < X_{n}.

The matrix Sym^{m}(g) is called m-th symmetric power matrix of the matrix g. We use the same symbol Sym^{m} for symmetric powers of vector spaces as well. We use the symbol S^{m} for the symmetric power of a differential system (cf. section 2.1.2); Sym^{m} is not the same matrix construction as S^{m} and this is why we must use a different notation.

Remark 2.2 If g is the matrix (g_{ij}) then Sym^{m}(g) can be computed as follows: Put v_{i} = ∑_{j} c_{j}g_{ij} and Y is the vector^{2} of monomials in the v_{i}. Then Sym^{m}(g)_{r,s} is found from the r-th entry of Y by taking the coefficient of the s-th monomial in the c_{j}.

2.2 Dual first integrals

Proposition 2.3 Let A be the companion matrix of a differential operator L, let G be the differential Galois group and let W be the solution space of Y′ = S^{m}(A)Y. There exists a G-isomorphism

λ : Sym^{m}(V (L)) → W.

2 Defining a vector of monomials implies choosing an ordering on the monomials, we take the same ordering as in remark 2.1.

Let U be a fundamental solution matrix of Y′ = AY. Then the columns of Sym^{m}(U) form a basis of W.

PROOF. Let K be the Picard-Vessiot extension generated by the entries of U. From remark 2.2 one can verify that Sym^{m}(U) satisfies the equation Y′ = S^{m}(A)Y. As Sym^{m} is a group homomorphism, Sym^{m}(U) is an invertible matrix and hence the second statement follows.

The entries of Sym^{m}(U) are in K, so W ⊂ K^{N} and hence the Galois group G acts on W. Let g ∈ G. Because K is the base field in the construction of the homomorphism Sym^{m} in the previous section, automorphisms of K commute with Sym^{m}, i.e. g(Sym^{m}(U)) = Sym^{m}(g(U)).

The automorphism g acts on U as multiplication on the right with a matrix g̃ ∈ Gl_{n}(\overline{C}). Let W_{2} be the solution space of Y′ = AY. The columns of U form a basis of W_{2}. On this choice of basis, the action of g on W_{2} is given by the matrix g̃. The action of g on V (L) is also the matrix g̃, where U_{1,1}, …, U_{1,n} is chosen as a basis for V (L). Then by definition of the matrix Sym^{m}(g̃), the action of g on Sym^{m}(V (L)) is given by the matrix Sym^{m}(g̃).

g(Sym^{m}(U)) = Sym^{m}(g(U)) = Sym^{m}(U · g̃) = Sym^{m}(U) · Sym^{m}(g̃)

So g acts on W as the matrix Sym^{m}(g̃), where the columns of Sym^{m}(U) are chosen as basis for W. So the matrix of the action of g is the same for W as for Sym^{m}(V (L)), hence W is G-isomorphic with Sym^{m}(V (L)). □

We can describe λ more explicitly as follows. After choosing a basis y_{1}, …, y_{n} of V (L), or equivalently, after choosing a fundamental solution matrix U of Y′ = AY, an element of Sym^{m}(V (L)) can be represented as a homogeneous polynomial P in the variables X_{1}, …, X_{n} of degree m. Let C ∈ \overline{C}^{N} be the vector of coefficients of P. Then

λ(P) = Sym^{m}(U) · C ∈ W. (2)

Note that in fact both Sym^{m}(V (L)) and W are defined independently of a choice of basis y_{1}, …, y_{n}, and that λ : Sym^{m}(V (L)) → W is also independent of this choice.

Via λ, an invariant in Sym^{m}(V (L)) can be presented as an element F ∈ W whose entries are left fixed by the Galois group; this is equivalent with saying that F ∈ W ⋂ k^{N}. An invariant given in this presentation (i.e. given as a rational solution F of Y′ = S^{m}(A)Y) will be called a dual first integral^{3} .

Lemma 2.4 Let P be a polynomial invariant and let C be the vector of its coefficients. Then, Sym^{m}(U) · C is a dual first integral. Conversely, let F be a dual first integral and let the vector C be such that F = Sym^{m}(U) · C. Then, C is the vector of coefficients of a polynomial invariant P. Moreover, the first entry of F is the value of P.

PROOF. The first two statements follow from equation (2) and the fact that λ is a G-isomorphism. For the third statement note that the first row in Sym^{m}(U) is the vector of all monomials in y_{1}, …, y_{n}. Hence the first entry of F equals P (y_{1}, …, y_{n}), i.e. P with y_{i} substituted for X_{i}. □

A different way to explain the relation between invariants and dual first integrals is given by the following proposition.

Proposition 2.5 Let K be the Picard-Vessiot extension and G the differential Galois group of L. Define the \overline{C}-algebra homomorphism

ϕ : Sym(V (L)) → K[X_{1}, …, X_{n}]

by (it suffices to define ϕ for homogeneous elements of degree 1)

ϕ(y) = ∑_{i=1}^{n} X_{i}y^{(i−1)} for y ∈ V (L).

Then ϕ is an embedding (as \overline{C}-algebra and as G-module) of Sym(V (L)) in K[X_{1}, …, X_{n}].

PROOF. ϕ gives an embedding (as \overline{C}-vector space and as G-module) of V (L) in K · X_{1} + ⋯ + K · X_{n}. Furthermore if y_{1}, …, y_{n} is a \overline{C}-basis of V (L) then their images form a K-basis of K · X_{1} + ⋯ + K · X_{n} because the Wronskian of y_{1}, …, y_{n} has non-zero determinant. Hence ϕ is an embedding of Sym\overline{C}(V (L)) in Sym_{K}(K · X_{1} + ⋯ + K · X_{n}) = K[X_{1}, …, X_{n}]. □

If we identify the K-vector space of homogeneous polynomials of degree m with K^{N} then the maps ϕ and λ from Sym^{m}(V (L)) to W ⊂ K^{N} coincide.

3 This name comes from the fact that solutions of the dual of Y′ = S^{m}(A)Y are first integrals for Y′ = AY.

2.3 Computational aspects

The operator whose solution space is spanned by all monomials of degree m in the y_{i} is noted L^{◯ s m} and is called the m-th symmetric power of L. Its solution space is a G-homomorphic image of Sym^{m}(V (L)) ([SUl93a]); V (L^{◯ s m}) = p(W) where p : K^{N} → K is the projection on the first component. The order of L^{◯ s m} is ≤ N = \binom{n+m−1}{n−1} (the number of monomials of degree m in n variables); it is < N if and only if there is a non-zero P ∈ \overline{C}[X_{1}, …, X_{n}], homogeneous of degree m, having value 0, i.e. P (y_{1}, …, y_{n}) = 0. In this case it can happen that the value of a homogeneous polynomial P of degree m is in k even though P is not an invariant. If order(L^{◯ s m}) = N then P is invariant if and only if its value is in k. So, the standard method for computing invariants is the following: replace L (if necessary) by an operator with an isomorphic solution space, in such a way that L^{◯ s m} has the correct dimension N. Then the set of values of the invariants of degree m is the space of rational solutions of L^{◯ s m}, cf. [SUl93b,SUl96]. The usual method (given in [SUl93a]) to construct L^{◯ s m} amounts to converting the system Y′ = S^{m}(A)Y to an equation by using the (putative) cyclic vector (1, 0, …, 0).

This method has three drawbacks. First, the cost of the computation of L^{◯ s m} grows very fast with m and n (because we must perform elimination on systems whose size grows exponentially). So, in practice, the computation becomes rapidly impossible. Secondly, if L^{◯ s m} does not have the right order, then one has to perform a transformation on L and re-do the whole computation (though some information can be saved, see [Wei95a]). And thirdly, for some applications, one indeed needs the invariant in the form of a dual first integral (e.g. [Com96]).

The first motivation of this paper was not to find a faster method, but a method that would work when computation of L^{◯ s m} fails, and that would avoid the above drawbacks. The approach in this paper consists in solving directly the system Y′ = S^{m}(A)Y, without converting it to an equation. For any point x_{0} ∈ P^{1}(\overline{C}), the system has a local formal fundamental solution matrix Sym^{m}(Û) where Û is a local solution matrix of Y′ = AY. The system Y′ = S^{m}(A)Y has a rational solution F (a dual first integral) if and only if there exists C ∈ \overline{C}^{N} such that Sym^{m}(Û)C = F ∈ k^{N}. We will use this in section 4 to compute F. Thanks to lemma 2.4, this will give us the invariants in all the presentations at the same time.

3 Bounds on exponents using generalized exponents

This section addresses the question of finding the denominator and a bound on the numerator of each entry of a dual first integral F. When computing rational solutions of a differential operator L, one first computes a lower bound for the integer exponents of L at each point x_{0} ∈ P^{1}(\overline{C}). We would like to compute rational solutions of symmetric powers (and other constructions) of differential operators. In the regular singular case, [SUl94] give the bound for the integer exponents of symmetric powers L^{◯ s m} in terms of the exponents of L. In the irregular singular case, however, we can not obtain a bound for the integer exponents of L^{◯ s m} from the exponents of L. The reason is that in this case there are “too few exponents”: in the irregular singular case, there are, counted with multiplicity, less than order(L) exponents. To handle this difficulty we will use a generalization of exponents. An alternative way to get a bound (a different bound than ours) is found in lemma 3.3 in [Sin93] using a different generalization of exponents found in [BBe85].

For convenience of notation we will now assume that the point of interest is the point x = 0. Then L in C(x)[∂] is viewed as an element of the ring C((x))[δ] = C((x))[∂] where

δ = x∂.

3.1 A few preliminaries on local differential operators

In this section we list a few known facts about local differential operators that will be used in later sections.

Definition 3.1 Let L = ∑ a_{i,j}x^{i}δ^{j} ∈ C((x))[δ] be non-zero and let T be a variable. Let v be the smallest integer such that a_{v,j} ≠ 0 for some j. Then the Newton polynomial N_{0}(L) for slope 0 of an operator L is defined as ∑ a_{v,j}T^{j} ∈ C[T].

If L can be written as a product L = L_{1} · L_{2} then N_{0}(L_{2}) is a factor of N_{0}(L). The Newton polynomial is used in algorithms for computing factorizations and/or formal solutions of differential operators. One defines a Newton polygon and for each slope in the Newton polygon a Newton polynomial can be defined. Definition 3.1 gives the Newton polynomial only for slope 0 in the Newton polygon. Definitions and properties of Newton polygons and polynomials can be found in [Mal79,Tou87,Bar88,Hoe95a].

Definition 3.2 The exponents of L are those elements e ∈ \overline{C} for which there is a solution of L of the form

x^{e}s where s ∈ \overline{C((x))}[log(x)] with v(s) = 0.

Here the valuation v(s) is defined as the smallest rational number such that the coefficient of x^{v(s)} in s is non-zero.

Note: If s ∈ \overline{C((x))}[log(x)] then s ∈ \overline{C}((x^{1/r}))[log(x)] for some integer r. The smallest r with this property is called the ramification index of s. The valuation v(s) for s ≠ 0 is the largest number in Q such that sx^{−v(s)} ∈ \overline{C}[[x^{1/r}]][log(x)]. The valuation of 0 is ∞.

The following is a well-known property of exponents. It is generalized in proposition 3.7.

Lemma 3.3 An element e ∈ \overline{C} is an exponent of L if and only if e is a root of N_{0}(L).

Note: In the literature exponents are often also called indices, and the Newton polynomial N_{0}(L) is then called the indicial polynomial or indicial equation.

We denote the linear universal extension of C((x)) by V. This is a ring that contains \overline{C((x))} and a basis of solutions of all homogeneous linear differential equations over C((x)). Furthermore V is minimal with this property. A construction is given in [HvP94], lemma 2.1.1. From the way that V is constructed in [HvP94] it follows that one can define a map

Exp : \overline{C((x))} → V

with the following properties: Exp(e) is a non-zero solution of δ − e, Exp(q) = x^{q} for q ∈ Q and

Exp(e_{1} + e_{2}) = Exp(e_{1})Exp(e_{2})

for e_{1}, e_{2} ∈ \overline{C((x))}, i.e. Exp behaves like an exponential function. One can think of Exp(e) as

Exp(e) ” = ” exp(∫ \frac{e}{x}dx) ” = ” x^{e}.

We have

Exp(e) ∈ C((x)) ⟺ e ∈ ZZ + x · C[[x]]

and

Exp(e) ∈ \overline{C((x))} ⟺ e ∈ \overline{C((x))} ⋂(⋃_{r}(\frac{1}{r}ZZ + x^{1/r} · \overline{C}[[x^{1/r}]])).

Definition 3.4 Define the substitution map

S_{e} : \overline{C((x))}[δ] → \overline{C((x))}[δ]

for e ∈ \overline{C((x))} as the \overline{C((x))}-homomorphism that maps δ to δ + e.

The substitution map has the following well-known property: Exp(e)y is a solution of L if and only if y is a solution of S_{e}(L).

3.2 Definition of generalized exponents

Using the substitution map, one can rewrite the standard property of exponents (lemma 3.3) as follows:

Lemma 3.5 Let L ∈ C((x))[δ] \ \{0\}. An element e ∈ \overline{C} is an exponent of L if and only if 0 is a root of the Newton polynomial N_{0}(S_{e}(L)).

With this lemma in mind, we can generalize the exponents by replacing the set \overline{C} by a larger set of exponents E. Define

E = ⋃_{r} \overline{C}[x^{−1/r}].

In the following definition we need to generalize definition 3.1 to non-zero elements of \overline{C((x))}[δ]. Take q ∈ Q minimal such that the coefficient of x^{q} in L is non-zero. Then N_{0}(L) is this coefficient (which is in \overline{C}[δ]) with δ replaced by the variable T.

Definition 3.6 Let L ∈ \overline{C((x))}[δ] \ \{0\}. For an element e ∈ E the number ν_{e}(L) is defined as the multiplicity of the root 0 in N_{0}(S_{e}(L)). e ∈ E is called a generalized exponent^{4} of L if ν_{e}(L) > 0. The number ν_{e}(L) is called the multiplicity of the generalized exponent e in the operator L.

Alternative approaches are found in the literature (e.g. [CLe55,Put95]). The exponents are those generalized exponents that are in \overline{C}.

4 A generalized exponent was called canonical exponential part in [Hoe95b]. We changed this name to point out the use of this notion, which is to generalize methods that use exponents (for example: [SUl94]) to the irregular singular case.

The generalized exponent should not be confused with the definition of exponential part in section 3.2 in [Hoe95b]. A generalized exponent is an element of the set E, whereas an exponential part is an element of the set E/ ∼. Here the equivalence ∼ is defined by

e_{1} ∼ e_{2} ⟺ e_{1} − e_{2} ∈ \frac{1}{ram(e_{1})}ZZ

where ram(e_{1}) is the ramification index of e_{1}. From the definition in section 3.2 in [Hoe95b] it follows that the multiplicity μ_{e_{1}}(L) of an exponential part e_{1} equals the sum of the multiplicities ν_{e_{2}}(L) of the generalized exponents e_{2} taken over all e_{2} ∈ E for which e_{2} ∼ e_{1}. So by theorem 1 in [Hoe95a] it follows that

∑_{e∈E} ν_{e}(L) = order(L). (3)

Many mathematical and algorithmic difficulties with irregular singular operators are caused by the fact that for such operators there are (counted with multiplicity) “too few” exponents. Because of equation (3) these difficulties no longer exist when using generalized exponents; for our purposes the irregular singular case is not different from the regular singular case.

Computing the generalized exponents can be done using one of the several factorization algorithms. It is a subproblem of computing formal solutions, so the generalized exponents can be computed using a part of the algorithm for computing formal solutions, cf. [Tou87,Bar88]. We use “algorithm semi-regular parts” in [Hoe95a]. This algorithm is a modified version of Malgrange’s factorization algorithm [Mal79]. It uses a different type of ramifications (obtained from [Bar88]) to minimize the algebraic extensions.

The relation between generalized exponents and formal solutions is the following (this is theorem 1 in [Hoe95b]):

Proposition 3.7 Let L ∈ C((x))[δ] \ {0}. An element e ∈ E is a generalized exponent of L if and only if L has a solution of the form

Exp(e)s where s ∈ \overline{C((x))}[log(x)] and v(s) = 0.

Note that, instead of using a Newton polynomial, the generalized exponents can be defined from the formal solutions using this proposition. A different generalization of exponents by using formal solutions is found in [BBe85].

3.3 Minimal exponents

As already mentioned, our reason for introducing generalized exponents was to obtain information about the exponents of L^{◯ s m} without computing the operator L^{◯ s m}. Now a natural question arises: Given the generalized exponents of L at the point x = 0, can one determine all (generalized) exponents of L^{◯ s m}? The answer is negative, as showed by the following example.

Example 3.8 Consider the operators L_{1} = ∂^{3} + x and L_{2} = ∂^{3} + x + 1. These operators are regular at x = 0. Hence both L_{1} and L_{2} have power series solutions with valuations 0, 1 and 2 at x = 0; the exponents at x = 0 are 0, 1, 2 for both operators. Making products of these solutions, one finds solutions of L_{1}^{◯ s 2} and L_{2}^{◯ s 2} with valuations 0, 1, 2, 3, 4. Hence L_{1}^{◯ s 2} and L_{2}^{◯ s 2} will have at least the exponents 0, 1, 2, 3, 4 at x = 0. However, not all exponents of L_{1}^{◯ s 2} and L_{2}^{◯ s 2} have been determined by this. L_{2}^{◯ s 2} is regular at x = 0, so it has exponents 0, 1, 2, 3, 4, 5. But L_{1}^{◯ s 2} has exponents 0, 1, 2, 3, 4, 6 at x = 0 (x = 0 is an apparent singularity, i.e. all solutions are analytic). The conclusion of this example is that the exponents 0, 1, 2 determine the smallest exponents of the second symmetric power, but not necessarily all exponents.

Let M be a differential operator whose solution space is spanned by differential monomials in the solutions of L. If L is regular at a point x = α (where α ∈ \overline{C}), then M need not be regular at x = α. However, products, sums and derivatives of analytic functions are analytic, hence all local solutions of M at x = α are analytic. It follows that all generalized exponents of M at x = α are integers, bounded from below by 0. So in this section we only need to compute lower bounds for the exponents of M at the singularities of L and the point infinity.

These remarks and the example suggest that, instead of trying to find all generalized exponents of symmetric powers of L, we should settle for a different goal, namely to compute the minimal generalized exponents.

Definition 3.9 Let r be a positive integer. Define the following partial ordering ≤_{r} on E

e_{1} ≤_{r} e_{2} ⟺ e_{1} − e_{2} ∈ \frac{1}{r}ZZ and e_{1} − e_{2} ≤ 0.

For a set S ⊂ E define min_{r}(S) as the set of minimal elements of S with respect to the ordering ≤_{r}. For an element L ∈ C((x))[δ] \ {0} define min_{r}(L) as min_{r}(S) where S is the set of generalized exponents of L.

If L has an integer exponent e ∈ ZZ then min_{r}(L) ⋂ \frac{1}{r}ZZ contains at least one element which is ≤ e. So if we can compute min_{r} for symmetric powers of L then we find lower bounds for the integer exponents of these symmetric powers. This is done by proposition 3.11 below, using the following definition:

Definition 3.10 The symmetric product^{5} of L_{1} and L_{2}, denoted by L_{1}◯ s L_{2}, is the monic differential operator of minimal order for which

y_{1}y_{2} ∈ V (L_{1}◯ s L_{2}) for all y_{1} ∈ V (L_{1}), y_{2} ∈ V (L_{2}).

The notation L^{(i)} denotes the monic differential operator defined by:

V (L^{(i)}) = {∂^{i}y|y ∈ V (L)}.

Proposition 3.11 Let L_{1} and L_{2} be non-zero elements of C((x))[δ]. Let r be the least common multiple of all ramification indices of the generalized exponents of L_{1} and L_{2}. Define for sets S_{1}, S_{2} ⊂ E the sum S_{1} + S_{2} as {s_{1} + s_{2}|s_{1} ∈ S_{1}, s_{2} ∈ S_{2}} ⊂ E. Then

min_{r}(L_{1}◯ s L_{2}) = min_{r}(min_{r}(L_{1}) + min_{r}(L_{2})).

Corollary 3.12 Let m be a positive integer and r be the ramification index of L. Denote for S ⊂ E the set m · S as S + S + ⋯ + S (m times). Then

min_{r}(L^{◯ s m}) = min_{r}(m · min_{r}(L)).

In particular if L^{◯ s m} has an integer exponent e then min_{r}(m · min_{r}(L)) ⋂ \frac{1}{r}ZZ contains one element which is a lower bound for e. This lower bound can be computed from r, m and min_{r}(L).

Remark 3.13 The fact that such a lower bound exists is not new (lemma 3.3 in [Sin93]). However, the bound in our proposition is sharper. It gives precisely the smallest exponent of L^{◯ s m} in \frac{1}{r}ZZ. So in case all ramification indices are 1 (i.e. r = 1) our bound for the smallest integer exponent is sharp (see also example 4.2 page 21).

We postpone the proof of the proposition till after the proof of theorem 3.15 below. To prove this theorem we first need to introduce some notations.

Denote^{6} V_{e} = Exp(e)·(\overline{C} ·C((x))[e])[log(x)] as in [Hoe95a,Hoe95b]. Note that \overline{C} · C((x))[e] = \overline{C} · C((x^{1/ram(e)})) where ram(e) is the ramification index of e.

5 Strictly speaking, this name is mathematically hazy. We use it to emphasize the resemblance with the symmetric power construction L^{◯ s m}. 6 Here \overline{C} · C((x)) denotes the smallest subfield of the algebraic closure of C((x)) that contains both \overline{C} and C((x)).

We have V_{e_{1}} = V_{e_{2}} if and only if e_{1} ∼ e_{2} and (cf. theorem 3 in [Hoe95a])

V = ⨁_{e∈E/∼} V_{e}. (4)

Now define

E_{r} = \overline{C}[x^{−1/r}] ⊂ E and V_{∗,r} = ⨁_{e∈E_{r}/∼} V_{e}.

For e ∈ E_{r} define

V_{e,r} = Exp(e) · (\overline{C} · C((x^{1/r})))[log(x)].

If e_{1} − e_{2} ∈ \frac{1}{r}ZZ then V_{e_{1},r} = V_{e_{2},r} so V_{e,r} can be defined for e ∈ E_{r}/(\frac{1}{r}ZZ). V_{e,r} is the direct sum of the V_{e_{1}} taken over all e_{1} ∈ E_{r}/ ∼ for which e is e_{1} modulo \frac{1}{r}ZZ. Hence by the direct sum in equation (4) it follows that

V_{∗,r} = ⨁ V_{e,r} (5)

where the sum is taken over all e ∈ E_{r}/(\frac{1}{r}ZZ).

Definition 3.14 The ramification index ram(L) of L ∈ C((x))[∂] is defined as the least common multiple of the ramification indices of all generalized exponents of L.

From theorem 3 in [Hoe95a], it follows that V (L) ⊂ V_{∗,r} if and only if ram(L) divides r. V_{∗,r} is a differential ring extension of \overline{C((x))} consisting of all solutions of all L ∈ C((x))[δ] for which ram(L) divides r. Hence if the ramification indices of two operators L_{1} and L_{2} divide the integer r then the same holds for the operators L_{1}◯ s L_{2}, L_{1}^{(1)} (definition 3.10) and for the least common left multiple of L_{1} and L_{2}.

Theorem 3.15 Let L ∈ C((x))[δ] \ {0} be of order d and let r be a positive integer. Suppose that the ramification indices of the generalized exponents divide the integer r.

(i) There exists a basis y_{1}, …, y_{n} of V (L) which satisfies the following condition

y_{i} = Exp(e_{i})s_{i} for some s_{i} ∈ (\overline{C} · C((x^{1/r})))[log(x)], v(s_{i}) = 0 (6)

where e_{1}, …, e_{n} ∈ E.

(ii) Suppose y_{1}, …, y_{n} is a basis of the solution space V (L) which satisfies condition (6). Then

min_{r}(L) = min_{r}({e_{1}, …, e_{n}}).

PROOF. Let e ∈ min_{r}(L). Since {e_{1}, …, e_{n}} is a subset of the set of all generalized exponents of L (there are at most order(L) = d different generalized exponents) it follows that the number of elements in min_{r}({e_{1}, …, e_{n}}) can not be larger than the number of elements in min_{r}(L). So we only need to prove that e ∈ min_{r}({e_{1}, …, e_{n}}). Without loss of generality we may assume that e_{i}−e ∈ \frac{1}{r}ZZ for i ≤ t and e_{i}−e ∉ \frac{1}{r}ZZ for i > t. We need to show that t ≠ 0 and that there is one i ≤ t with e_{i} − e ≤ 0. Then the theorem is proven as follows: We may assume that e_{i} − e ∈ \frac{1}{r}ZZ is minimal, so e_{i} ∈ min_{r}(\{e_{1}, …, e_{n}\}). Because of the minimality of e we can not have e_{i} − e < 0 hence e = e_{i}.

The algorithm in section 8.4 in [Hoe95a] produces a basis ỹ_{1}, …, ỹ_{n} of formal solutions (see also the proof of theorem 1 in [Hoe95b]) where each basis element can be written in the form ỹ_{i} = Exp(ẽ_{i})s̃_{i} with s̃_{i} ∈ C((x))[ẽ_{i}, log(x)] and v(s̃_{i}) = 0 and where every generalized exponent ẽ_{i} of L occurs. Since (C((x))[ẽ_{i}])[log(x)] ⊂ (\overline{C} · C((x^{1/r})))[log(x)] this basis satisfies condition (6). Furthermore the generalized exponent e of L occurs in this basis so one of the elements of this basis is of the form y = Exp(e)s (with s ∈ C((x))[e, log(x)] and v(s) = 0). Then y ∈ V_{e,r} and y ∈ V (L).

Because of condition (6), each y_{i} is an element of V_{e_{i},r}. Since the y_{i} form a basis of V (L), it follows that y is a \overline{C} linear combination of y_{1}, …, y_{n}. Because of the direct sum in equation (5), it follows that y is a linear combination of y_{1}, …, y_{t}, since e_{i} for i > t is not equal to e modulo \frac{1}{r}ZZ and so y_{i} is in a different component than V_{e,r} for i > t. Dividing by Exp(e), we obtain that s ∈ V_{0,r} = (\overline{C} ·C((x^{1/r})))[log(x)] is a linear combination of the Exp(e_{i} −e)s_{i} ∈ V_{0,r} for i ≤ t. Hence the valuation of at least one of the Exp(e_{i} −e)s_{i} is ≤ v(s) = 0. The valuation of the s_{i} is 0 and the valuation of Exp(e_{i} − e) ∈ C((x^{1/r})) is e_{i} − e. So for at least one i ≤ t we have e_{i} − e ≤ 0 and so the theorem follows. □

Remark 3.16 Without the condition s_{i} ∈ (\overline{C}·C((x^{1/r})))[log(x)] the statement need not hold. Take for example L = δ^{2}−\frac{1}{2}δ and r = 1. Then min_{r}(L) = \{0, \frac{1}{2}\}. Now take e_{1} = e_{2} = 0, s_{1} = 1 and s_{2} = 1 + x^{1/2}. Then s_{2} does not satisfy condition (6) and min_{r}(L) ≠ min_{r}({e_{1}, e_{2}}) = {0}.

Remark 3.17 The existence result i) is also found in [CLe55] (with a different terminology, though)

[ Proof of proposition 3.11:] Let y_{i} = Exp(e_{i})s_{i}, i = 1, …, order(L_{1}) be a basis of V (L_{1}) and ỹ_{j} = Exp(ẽ_{j})s̃_{j}, j = 1, …, order(L_{2}) be a basis of V (L_{2}) which both satisfy condition (6). Then the products y_{i}ỹ_{j} span V (L_{1}◯ s L_{2}). Let S be a set of pairs (i, j) such that {y_{i}ỹ_{j}|(i, j) ∈ S} is a basis for V (L_{1}◯ s L_{2}). Now y_{i}ỹ_{j} = Exp(e_{i} + ẽ_{j})s_{i}s̃_{j} and s_{i}s̃_{j} ∈ (\overline{C} · C((x^{1/r})))[log(x)] with v(s_{i}s̃_{j}) = 0. Hence by theorem 3.15 it follows that min_{r}(L_{1}◯ s L_{2}) = min_{r}({e_{i} + ẽ_{j}|(i, j) ∈ S}).

Now {e_{i} + ẽ_{j}|(i, j) ∈ S} is a subset of the set T of all e_{i} + ẽ_{j}. So for each e ∈ min_{r}({e_{i} + ẽ_{j}|(i, j) ∈ S}) there must be precisely one ẽ ∈ min_{r}(T) such that ẽ ≤_{r} e. Furthermore T is a subset of the set of all generalized exponents of L_{1}◯ s L_{2}. Hence for each ẽ ∈ min_{r}(T) there must be precisely one e ∈ min_{r}(L_{1}◯ s L_{2}) = min_{r}({e_{i} + ẽ_{j}|(i, j) ∈ S}) such that e ≤_{r} ẽ. Then it follows that min_{r}(T) equals min_{r}(L_{1}◯ s L_{2}). □

Lemma 3.18 Let L ∈ C((x))[δ] be non-zero and let r be the ramification index of L. If 0 ∉ min_{r}(L) then

min_{r}(L^{(1)}) = {e + v(e) − 1|e ∈ min_{r}(L)}.

If 0 ∈ min_{r}(L) then

min_{r}(L^{(1)}) = {m} ⋃{e + v(e) − 1|e ∈ min_{r}(L) \ {0}}

where m ∈ ZZ, m ≥ −1, or

min_{r}(L^{(1)}) = {e + v(e) − 1|e ∈ min_{r}(L) \ {0}}.

Note that order(L^{(1)}) = dim(∂(V (L))) = dim(V (L)) − dim(V (L) ⋂ V (∂)). So order(L^{(1)}) = order(L) − 1 if 1 ∈ V (L) and order(L^{(1)}) = order(L) otherwise.

PROOF. If y = Exp(e)s where s ∈ (\overline{C} · C((x^{1/r})))[log(x)] with v(s) = 0 and e ≠ 0 then the derivative y′ is of the form Exp(e + v(e) − 1)t for some t ∈ (\overline{C} · C((x^{1/r})))[log(x)] with v(t) = 0. Now the first statement follows by applying theorem 3.15.

For the second statement we note that ν_{0}(L) > 0 means that there is a formal solution y ∈ (\overline{C} · C((x^{1/r})))[log(x)] of L with v(y) = 0. The valuation of the derivative y′ is ∞ or is an integer ≥ −1. Now distinguish the two cases: v(y′) ∈ min_{r}(L^{(1)}) (then: v(y′) is an integer m ≥ −1) or v(y′) ∉ min_{r}(L^{(1)}) (then the other case holds). □

Note that in the case 0 ∈ min_{r}(L) one can get a slightly stronger statement about min_{r}(L^{(1)}) by noting that −1 ∈ min_{r}(L^{(1)}) if and only if ν_{0}(L) > 1. We will not use this small improvement of the lemma.

Define v′ : E → Q as follows: v′(e) = v(e) for all e ∈ E \ {0} and v′(0) = 0. It follows from the lemma that for each e ∈ min_{r}(L^{(1)}) there is an ẽ ∈ min_{r}(L) such that e − (ẽ + v′(ẽ) − 1) is a non-negative integer. Repeating this, we see that for each e ∈ min_{r}(L^{(i)}) there is an ẽ ∈ min_{r}(L) such that e−(ẽ+i·v′(ẽ)−i) is a non-negative integer.

Theorem 3.19 Let L be a non-zero differential operator in C((x))[δ] and r be its ramification index. Let m_{0}, …, m_{n−1} be non-negative integers and M the symmetric product of the operators (L^{(i)})^{◯ s m_{i}}. Define B_{i} = {e+i·v′(e)−i}|e ∈ min_{r}(L)} and B = m_{0} · B_{0} + ⋯ + m_{n−1} · B_{n−1}. Suppose M has a non-zero solution y in (\overline{C} ·C((x^{1/r})))[log(x)]. Then B ⋂ \frac{1}{r}ZZ contains an element ≤ v(y).

The theorem gives a lower bound for the valuation of solutions of M in (\overline{C} · C((x^{1/r})))[log(x)]. The bound can be computed from m_{0}, …, m_{n−1}, r and min_{r}(L).

To compute the bound we need to compute the set of sums m_{0} · B_{0} + ⋯ + m_{n−1} · B_{n−1} and to take the smallest element which is in \frac{1}{r}ZZ. This means computing in a splitting field: it is not sufficient to take only one generalized exponent in each conjugacy class of generalized exponents. One can try to avoid splitting fields for computing this bound by various tricks (for example floating point computations) but we will not go into this.

In the following procedure, the notation l_{α} denotes the \overline{C}-automorphism of \overline{C}(x)[∂] given by l_{α}(x) = x + α and l_{α}(∂) = ∂; this transformation moves the point x = α to x = 0. Similarly, l_{∞} is a \overline{C}-automorphism of \overline{C}(x)[∂] given by l_{∞}(x) = 1/x and l_{∞}(∂) = −x^{2}∂; this moves the point infinity to x = 0.

Algorithm 1 (Procedure global-bounds) — Input: An operator L ∈ C(x)[∂], and non-negative integers m_{0}, …, m_{n−1} Output: A rational function Q ∈ C(x) and an integer N such that every rational solution y ∈ \overline{C}(x) of M = (L^{(0)})^{◯ s m_{0}}◯ s ⋯ ◯ s (L^{(n−1)})^{◯ s m_{n−1}} can be written as the product of Q and a polynomial in x of degree ≤ N.

(i) Q := 1. (ii) After multiplication on the left by an element of C[x], we may assume

that L = a_{n}∂^{n} + ⋯ + a_{0}∂^{0} with a_{i} ∈ C[x] and gcd(a_{0}, …, a_{n}) = 1. (iii) For each irreducible factor p of a_{n} in C[x] (not \overline{C}) do

(a) Let α ∈ \overline{C} be a root of p.

(b) Compute the generalized exponents of l_{α}(L) at the point x = 0.

(c) Let r be the ramification index of l_{α}(L) at x = 0; compute the min_{r} of the set of generalized exponents.

(d) Compute the set B from theorem 3.19.

(e) If B ⋂ \frac{1}{r}ZZ is empty then stop the algorithm and RETURN the following output: Q = 0 and N = 0.

(f) Let b_{α} ∈ ZZ be the smallest element of B ⋂ \frac{1}{r}ZZ, rounded above to an integer.

(g) Replace Q by Q · p^{b_{α}}. (iv) Perform steps 3b, 3c, 3d, 3e, 3f with α = ∞.

(v) Add 2 · (0 · m_{0} + 1 · m_{1} + ⋯ + (n − 1) · m_{n−1}) to b_{∞}. (vi) Let N be −b_{∞} plus the valuation of Q at infinity (this valuation is the

degree of the denominator of Q minus the degree of the numerator of Q). (vii) Return: Q and N.

Remark 3.20 Note that, even if Q = N = 0, there may be an invariant (whose value is zero): see the Hurwitz example in the next section for an illustration.

The fact that the algorithm works follows from the following observations:

– Because algebraic conjugation over C is an automorphism of the differential field \overline{C}(x), it follows that if α_{1}, α_{2} ∈ \overline{C} are conjugated over C then the two bounds b_{α_{1}}, b_{α_{2}} ∈ ZZ will be the same. Hence we need to take only one α in every conjugacy class of the singularities of L. In other words: We need to compute the bound for only one root of each factor of a_{n} in C[x]. Furthermore the function Q ∈ \overline{C}(x) will be an element of C(x).

– Note that for all α ∈ P^{1}(\overline{C}) the map l_{α} on \overline{C}(x)[∂] commutes with taking symmetric products and LCLM’s (least common left multiples) because the map l_{α} on \overline{C}(x) commutes with multiplication and addition. However, l_{α} does not commute with derivation if α = ∞. So l_{α} only commutes with the construction L ↦ L^{(1)} on \overline{C}(x)[∂] if α ∈ \overline{C}. The solution space of l_{∞}(L^{(1)}) equals x^{2} times the solution space of l_{∞}(L)^{(1)}, so the valuations are 2 higher than in lemma 3.18. For the point x = ∞ there is a lemma similar to lemma 3.18 with the following differences: e + v′(e) − 1 is replaced by e + v′(e)+1 and m ≥ −1 is replaced by m ≥ 1. We need a different theorem 3.19 specifically for the point x = ∞, i.e. for operators L ∈ C((\frac{1}{x}))[∂] instead of L ∈ C((x))[∂]. The only difference will be that e + i · v′(e) − i needs to be replaced by e + i · v′(e) + i. The algorithm computes the bound given by theorem 3.19 and then adds 2 · (0 · m_{0} + 1 · m_{1} + ⋯ + (n − 1) · m_{n−1}) to correct for this difference.

Example 3.21 (The P SL_{3} example:) The following example was adapted from N. Katz by Elie Compoint ([Com96]). The component of the identity of its Galois group is P SL_{3}(C) in its 8-dimensional representation. Let δ = x\frac{d}{dx} and

L = δ(δ − \frac{1}{2})(δ − \frac{1}{4})(δ + \frac{1}{4})(δ − \frac{1}{8})(δ − \frac{5}{8})(δ + \frac{1}{8})(δ + \frac{5}{8}) − x(δ + \frac{1}{3})(δ − \frac{1}{3}).

We want to compute the invariants of degree 2 and 3.

The generalized exponents of l_{∞}(L) are −1/3, 1/3 and all conjugates of x^{−1/6}+ 2/3. The ramification index is r = 6. Now −1/3 ≤_{r} 1/3 and all the other generalized exponents are different modulo \frac{1}{r}ZZ. Hence min_{r}(l_{∞}(L)) contains 7 elements; all generalized exponents except 1/3. Now the smallest element in (\frac{1}{r}ZZ) ⋂(2 · min_{r}(l_{∞}(L))) is −2/3. Rounded above to an integer this is 0. The smallest element in (\frac{1}{r}ZZ) ⋂(3 · min_{r}(l_{∞}(L))) is −1.

The generalized exponents (which are in fact exponents) of L at x = 0 are −5/8, −1/4, −1/8, 0, 1/8, 1/4, 1/2 and 5/8. So the ramification index is r = 1. Since all exponents are different modulo \frac{1}{r}ZZ the set min_{r}(L) equals the set of exponents. Now the smallest element in (\frac{1}{r}ZZ) ⋂(2 · min_{r}(L)) is 0 and the smallest element in (\frac{1}{r}ZZ) ⋂(3 · min_{r}(L)) is −1.

So the procedure global-bounds gives the following output for the second symmetric power of L: Q = 1 and N = 0. This means that the values of all invariants of degree 2 are constants. For the third symmetric power the output is Q = 1/x and N = 2, which means that the values of the invariants of degree 3 must be of the form \frac{1}{x} ·(c_{0}x^{0} +c_{1}x^{1} +c_{2}x^{2}) for some constants c_{0}, c_{1}, c_{2} which will be computed in the next section.

4 The algorithm for computing invariants

We now have all ingredients for an algorithm. There exists an invariant of degree m if and only if there is a rational solution F to Y′ = S^{m}(A)Y. The previous section gives the denominators and bounds for the degrees of the numerators of the entries of F. Thus, the problem can be reduced to linear algebra. To obtain the numerators in F, we consider a local fundamental solution matrix Û at some (possibly singular) point x_{0} ∈ P^{1}(\overline{C}). We can assume x_{0} = 0 in our algorithm after having applied the map l_{x_{0}}. Now F = Sym^{m}(Û) · C for some constant vector C. We start with undetermined constants in C, compute sufficiently many terms of the power series in Sym^{m}(Û) and then express the numerators in F in terms of the constants in C, see below. As the evaluation of the series is usually the most costly part of the algorithm, our main goal below will be to reduce the number of columns and rows of Sym^{m}(Û) that need to be evaluated during the process.

4.1 Computing candidate invariants

First, we perform the above idea only on the first row of Sym^{m}(Û). The reason for choosing this particular row is that the bounds in theorem 3.19 (compare the generalized exponents of L and L^{(1)}) are smaller than the bounds for other rows (unless there are less than 3 singularities).

Suppose ŷ_{i}, 1 ≤ i ≤ n = order(L), is a basis of formal solutions satisfying condition (6) on page 15. Then a monomial in these ŷ_{i} (i.e. a product ∏(ŷ_{i})^{m_{i}}) is again of the form (6), where the generalized exponent equals ∑ m_{i}e_{i}. The following lemma reduces the number of columns of Sym^{m}(Û) that need to be evaluated.

Lemma 4.1 Let ŷ_{i} be a basis of local formal solutions satisfying condition (6) on page 15 and let r be the ramification index. An entry of a vector C of coefficients of an invariant can only be non-zero if the generalized exponent of the corresponding monomial is in \frac{1}{r}ZZ.

PROOF. Let N = \binom{n+m−1}{n−1} and let Û be a formal fundamental matrix of Y′ = AY such that the first row is ŷ_{1}, …, ŷ_{n}, i.e. the entries of Û are the 0, …, (n−1)-th derivatives of ŷ_{1}, …, ŷ_{n}. Let P be a polynomial invariant and C be the vector of its coefficients. Then Sym^{m}(Û) · C ∈ \overline{C}(x)^{N} ⊂ (V_{0,r})^{N}. Note that each column of Sym^{m}(Û) is an element of (V_{e,r})^{N} where e is the generalized exponent of the first element (which is a monomial in the ŷ_{i}) of this column. From the fact that the columns of Sym^{m}(Û) are linearly independent and the direct sum (5) on page 15, it follows from Sym^{m}(Û)C ∈ (V_{0,r})^{N} that C can only have a non-zero entry for those columns which are in (V_{0,r})^{N}, i.e. for those monomials whose generalized exponent is in \frac{1}{r}ZZ. □

Note that the above lemma is sharp, i.e. we must consider the generalized exponents in \frac{1}{r}ZZ. Taking only generalized exponents in ZZ is not sufficient as is shown by the following example.

Example 4.2 Let L ∈ Q((x))[∂] be the monic operator of order 4 which has the following local solutions at x = 0

y_{1} = Exp(\frac{1}{x^{1/2}} + \frac{1}{10}) · (1 + x^{1/2})

y_{2} is the conjugate (replace x^{1/2} by −x^{1/2}) of y_{1}

y_{3} = Exp(\frac{1}{x^{1/2}} + \frac{4}{10})

and y_{4} is the conjugate of y_{3}. The ramification index of L is 2. L has an invariant of degree 2, even though none of the monomials y_{i}y_{j} has a generalized exponent in ZZ. The monomials y_{1}y_{4} and y_{2}y_{3} have generalized exponent 1/2. And, in fact, y_{1}y_{4} − y_{2}y_{3} = 2x is the value of an invariant of degree 2.

Algorithm 2 (Heuristic for computing invariants) — Input: an operator L, an integer m, a point x_{0}, and a number ν. Output: a vector space of candidate invariants of degree m and their corresponding candidate values, given as a parametrized candidate vector invariant and candidate value.

(i) If x_{0} ≠ 0 then apply recursion on l_{x_{0}}(L) as follows: replace L and x_{0} by l_{x_{0}}(L) and 0, apply this algorithm and then apply the inverse of l_{x_{0}} on the candidate values of the invariants.

(ii) Use the procedure global-bounds (page 18) to find the bounds Q_{1}, N_{1} for

rational solutions of the m-th symmetric power of L. (iii) Compute a basis of formal solution ŷ_{i} at x = 0 having property (6) in

theorem 3.15. Let r be the ramification index.

Let Y denote the vector of all monomials of degree m in the ŷ_{i}. Each of these monomials has a generalized exponent in \overline{C}[x^{−1/r}].

(iv) Take a vector C of unknown constants and set to zero every entry corresponding to a monomial with a generalized exponent that is not an element of \frac{1}{r}ZZ.

(v) Compute p_{1} := \frac{1}{Q_{1}} Y · C mod x^{N_{1}+ν+1} (vi) Build a linear system on C by equating to zero every term with degree

higher than N_{1} and all terms involving a log or a non-integer power of x. (vii) Return: the solution of this system (this is a vector space consisting

of candidate vector invariants) and the corresponding (vector space of) rational functions f_{1} := p_{1}Q_{1}.

Proposition 4.3 Denote by W_{L,m,ν} the vector space of candidate vector invariants produced by the above heuristic. Denote W_{L,m,∞} = ⋂_{ν} W_{L,m,ν}. Then:

(i) For all ν ∈ N, any vector of coefficients of an invariant of degree m is in W_{L,m,ν}.

(ii) There exists ν_{0} ∈ N such that W_{L,m,∞} = W_{L,m,ν_{0}}.

PROOF. Recall from section 3 that the value of any invariant of degree m is the product of Q_{1} by a polynomial of degree at most N_{1}. According to lemma 4.1, we have only computed necessary (but in general not sufficient) linear conditions. Hence i) follows. Increasing ν adds more conditions on C so W_{L,m,i+1} ⊂ W_{L,m,i}. As W_{L,m,∞} is finite dimensional, this implies ii). □

This algorithm is called a heuristic because there is no easy way of deciding whether its output is the vector space of invariants or a larger vector space. The number ν is can be chosen arbitrarily; the strategies of choice will be discussed in the next section.

Remark 4.4 We have order(L^{◯ s m}) < N = \binom{n+m−1}{n−1} if and only if the solutions of L satisfy a homogeneous polynomial relation of degree m. In this case, the value of a non-zero invariant can be zero, and furthermore it can also happen that W_{L,m,∞} contains elements that are not invariants (see the F_{36} example page 27). Note that since we do not compute L^{◯ s m} we have no easy way of checking if this case order(L^{◯ s m}) < N occurs, so this would be a serious problem if we only had the heuristic to compute invariants. We do not have this problem if we use the algorithm Invariants below; then the case order(L^{◯ s m}) < N does not cause difficulties anymore.

Example 4.5 (The P SL_{3} example (continued)) Let L be the 8-th order operator in the P SL_{3} example in section 3. We had found the bounds for rational solutions of L^{◯ s 2} and L^{◯ s 3}. Applying the heuristic with x_{0} = 0 and ν = 10 the following (candidate) invariants are obtained:

P_{2} = \frac{352}{32805}c_{0}X_{1}X_{7} − \frac{3249799168}{215233605}c_{0}X_{4}^{2} − \frac{36064}{6561}c_{0}X_{3}X_{5}

+ \frac{20240}{6561}c_{0}X_{6}X_{2} + \frac{12397}{3645}c_{0}X_{8}^{2}

and P_{2}(ŷ) = c_{0}

P_{3} = −\frac{15167488}{405}c_{1}X_{8}X_{6}X_{2} + \frac{35500589056}{12301875}c_{1}X_{1}^{2}X_{6}

+ \frac{659456}{10125}c_{1}X_{1}X_{8}X_{7} − \frac{36929536}{54675}c_{1}X_{7}X_{4}X_{5}

− \frac{743206912}{22275}c_{1}X_{8}X_{3}X_{5} + \frac{106172416}{3267}c_{1}X_{2}X_{3}^{2}

− \frac{3479057727488}{81192375}c_{1}X_{1}X_{3}X_{4} + \frac{46450432}{3375}c_{1}X_{8}^{3}

+ \frac{12176702046208}{66430125}c_{1}X_{8}X_{4}^{2} + \frac{424689664}{18225}c_{1}X_{6}X_{5}^{2}

+ c_{1}X_{7}^{2}X_{2}

and P_{3}(ŷ) = c_{1} \frac{1}{x}(1 + \frac{9144576}{3025}x + \frac{17832200896512}{3826625}x^{2})

where c_{0}, c_{1} denote arbitrary constants. Note that L^{◯ s 2} and L^{◯ s 3} have order 36 and 120 respectively. L^{◯ s 2} can be computed in several days on a big computer and L^{◯ s 3} is out of reach of computers, whereas the above computation only takes a few minutes.

4.2 Strategies for the heuristic

In the heuristic the point x_{0} and the number ν can be chosen. The advantage of choosing a singular point x_{0} is that the number of monomials that need to be considered (lemma 4.1) in the heuristic is often smaller, and so we need to evaluate fewer columns of Sym^{m}(Û). This still holds (and is important for the efficiency) for the algorithm Invariants below.

Example 4.6 (The P SL_{3} example (continued again)) In the P SL_{3} example of section 3, if we would take a regular point x_{0} then the heuristic would need to evaluate 36 monomials for the invariants of degree 2, and 120 monomials for degree 3. However, when taking the singularity x_{0} = 0, only 5 monomials of degree 2 have an integer exponent (the algorithm only considers monomials with an exponent in \frac{1}{r}ZZ, and r = 1 in this example). And only 15 monomials of degree 3 have an integer exponent. So when using the singular point x_{0} = 0 the computation for both the heuristic and the algorithm is much quicker than, say, with the regular point x_{0} = 1.

Taking a point in which a ramification occurs can be disadvantageous, because computing modulo x^{N} in C[[x^{1/r}]] involves more coefficients in C than computing modulo x^{N} in C[[x]]. So, the point x_{0} = 0 (ramification index is 1) in the P SL_{3} example is more favorable than the point x_{0} = ∞ (ramification index is 6). A point where the generalized exponents require algebraic extensions can have both advantages and disadvantages. The disadvantage is obvious: Computing the formal solutions and evaluating monomials will be more costly. The advantage is that many monomials need not be considered, for example:

Example 4.7 Suppose that order(L) = 3 and that at the point x_{0} = 0 we have 3 generalized exponents e_{1}, e_{2}, e_{3} which are algebraic over C((x)) of degree 3. From c_{1}e_{1} + c_{2}e_{2} + c_{3}e_{3} ∈ \frac{1}{r}ZZ and c_{1}, c_{2}, c_{3} ∈ ZZ it follows that c_{1} = c_{2} = c_{3} and hence only 1 monomial needs to be considered. So, for order 3, what would appear to be the worst case (the e_{i} are algebraic of degree 3), is in fact a relatively easy case.

By reasoning as in section 1.c of [Sin93], an application of Cramer’s formulas shows that we can take the following value for the number ν_{0} in proposition 4.3: N(1 + (N − 1)d + Nd_{1}) (where N = \binom{n+m−1}{n−1}, d is the maximum degree of the a_{i}, and d_{1} bounds the degrees of the numerator and denominator of Q_{1}). Thus, the above heuristic could be turned into an algorithm (but then the kernel problem order(L^{◯ s m}) < N of remark 4.4 would need to be addressed as well). However, this value for ν_{0} is usually much larger than necessary. So it is more efficient first to use the heuristic with a small value of ν, and then to apply the full algorithm Invariants from section 4.3.

If one already has some information about the group then sometimes the heuristical algorithm is sufficient to compute the invariants. Because if we know how many linearly independent invariants of degree m exist, we can simply use the heuristic by just increasing the value of ν. If at a certain point the space of candidate vector invariants has the correct dimension then it is certain (even in the problem case order(L^{◯ s m}) < N) that all invariants have been determined because the invariants form a subspace of the candidate invariants. In practice, the required number ν is usually much smaller than the theoretical bound ν_{0} above.

Example 4.8 (The Hurwitz example) The following operator has Galois group G_{168} ([SUl94]). Let ∂ = \frac{d}{dx} and

L = ∂^{3} + \frac{7x − 4}{x(x − 1)}∂^{2} + \frac{2592x^{2} − 2963x + 560}{252x^{2}(x − 1)^{2}}∂ + \frac{−40805 + 57024x}{24696x^{2}(x − 1)^{2}}.

The Galois group has invariants of degree 4, 6, 14, 21. The heuristic with m = 4, x_{0} = ∞, ν = 10 yields (in 0.75 seconds) a one-dimensional space generated by P = 1728X_{1}X_{2}^{3} + X_{1}^{3}X_{3} − 1728X_{2}X_{3}^{3} together with the value 0. The fact that the space of invariants of degree 4 has dimension exactly 1 proves that P is indeed an invariant. Similarly, the heuristic yields the other invariants quickly (see also [Wei95b]): for the invariant of degree 21, we need to compute 37 monomials at infinity (using a regular point it would have been 253 monomials).

4.3 Finding and proving which candidates are invariants

Let the monomial μ be a product of y^{(i)} to the power m_{i}, i = 0, …, n − 1. If y is a solution of L then μ is a solution of the symmetric product of the operators (L^{(i)})^{◯ s m_{i}}. By “applying procedure global-bounds on μ”, we mean applying the procedure global-bounds on these numbers m_{0}, …, m_{n−1}.

Algorithm 3 (Algorithm Invariants) — Input: L, m, x_{0} (optional: ν) Output: the space of invariants in vector and dual forms

(i) Like in the heuristic, if x_{0} ≠ 0 then apply the transformation l_{x_{0}}, use recursion, and transform back.

(ii) Now we may assume x_{0} = 0. Compute a basis of formal solutions of L at the point x = 0 having property (6) in theorem 3.15. Construct Û, the fundamental solution matrix of Y′ = AY from this.

(iii) Obtain F_{1} = f and C from the heuristic (page 22). Note that f and C contain parameters.

(iv) for i from 2 to N do:

– Let μ_{i} be the i-th monomial of degree m in y, y′, …, y^{(n−1)}.

Obtain Q_{i} and N_{i} from procedure global-bounds (page 18) applied to L and μ_{i}.

– Let p_{i} := \frac{1}{Q_{i}} Sym^{m}(Û)_{i}C mod x^{N_{i}+1} and F_{i} := p_{i} · Q_{i}.

Sym^{m}(Û)_{i} denotes the i-th row of Sym^{m}(Û).

– equate all terms involving logarithms or fractional powers of x to 0 (this gives a set of linear equations in the parameters. If the equations are non-trivial we use them to reduce the number of parameters).

(v) Now F is a vector of rational functions and C is a vector of constants. F and C contain parameters. The relation F′ − S^{m}(A)F = 0 yields a system of linear equations in the parameters. Solve this system.

(vi) Output: a basis of solutions C of this system and the corresponding dual first integrals F ∈ k^{N}.

Theorem 4.9 The output of this algorithm is exactly the space of invariants of degree m and the corresponding dual first integrals.

PROOF. That any vector of coefficients of an invariant is an element of the vector space produced by the algorithm follows from the fact that this was true for our heuristic, and from the fact that we only added necessary linear conditions in this algorithm. Hence also every dual first integral F is an element of the vector space produced by the algorithm. Conversely, as the F produced by the algorithm are rational vectors satisfying F′ = S^{m}(A)F, they are dual first integrals. So, by lemma 2.4, the corresponding C are indeed vectors of coefficients of invariants. □

4.4 Improvements and variants

Lemma 4.1 provided a speedup in the algorithm; it reduces the number of rows of Sym^{m}(Û) that need to be evaluated. It turns out that the number of columns that need to be evaluated can be reduced as well, using the following observation: F′ = S^{m}(A)F is not a random system of differential equations; there are recurrence relations so that, once some entries of F are known, the other entries can be deduced straightforwardly. Hence these latter entries, and the corresponding rows in Sym^{m}(Û), need not be considered in step 4 in algorithm 3. This provides a significant improvement of our algorithm, see below.

The recurrence relations are found as follows: Let y be a solution of L(y) = 0 and let μ_{[m_{1},...,m_{n}]} = y^{m_{1}} ⋯ (y^{(n−1)})^{m_{n}} with m_{1} + … + m_{n} = m denote a differential monomial of degree m. Then:

μ′_{[m_{1},...,m_{n}]} = m_{1}μ_{[m_{1}−1,m_{2}+1...,m_{n}]} + … + m_{n−1}μ_{[m_{1},...,m_{n−1}−1,m_{n}+1]}

+m_{n} (− ∑_{j=1}^{n} \frac{a_{j−1}}{a_{n}} μ_{[...,m_{j}+1,...,m_{n}−1]})

(with the convention that μ_{[..,−1,..]} = μ_{[..,m+1,..]} = 0). So, by replacing m_{n−1} by m_{n−1} + 1 in the above, we have:

μ_{[m_{1},...,m_{n−1},m_{n}+1]} = \frac{1}{m_{n−1}+1}(μ_{[m_{1},...,m_{n−1}+1,m_{n}]} − m_{1}μ_{[m_{1}−1,m_{2}+1,...,m_{n−1}+1,m_{n}]}

− ⋯ + m_{n} ∑_{j=1}^{n} \frac{a_{j−1}}{a_{n}} μ_{[...,m_{j}+1,...,m_{n−1}+1,m_{n}−1]}) (∗)

These relations must be understood as relations among the components of a solution F of F′ = S^{m}(A)F. So, if we know all the entries of F corresponding to monomials μ_{[...,i]} (and μ_{[...,i−1]} if i > 0), then the relations (∗) provide the entries corresponding to the μ[…, i + 1]. Thus, in the algorithm, we only need to evaluate the rows corresponding to the μ[…, 0]. A detailed study of the recursion shows that in the same way the rows corresponding to some of the μ_{[...,0]} can be skipped as well. This way, instead of evaluating N = \binom{n+m−1}{n−1}, we only need to evaluate^{7} \binom{n+m−3}{n−2} rows of Sym^{m}(Û) in step 3 of algorithm 3. Also note that, unless the number of singularities is < 3, the rows corresponding to μ_{[...,i]} have larger bounds as i increases (theorem 3.19 or the example below). So, the above relations allow us to skip the evaluation of most of the rows, but what is even more important for the efficiency is that we can skip the rows that have the worst bounds.

Step 5 of the algorithm can also be improved along the same lines. In the above process, we used the derivatives of the rows corresponding to the μ_{[m_{1},...,m_{n−1},m_{n}]} with m_{n−1} ≠ 0 to construct the entries of F. Thus, the corresponding relations are automatically satisfied in step 5 of algorithm 3. Hence in step 5 we only need to consider the rows corresponding to monomials μ_{[...,m_{n−2},0,m_{n}]}. This eliminates a lot of redundant equations.

Example 4.10 (The F_{36} example) Let

L = ∂^{3} + \frac{5 (9x^{2} + 14x + 9) ∂}{48(x + 1)^{2}x^{2}} − \frac{5 (81x^{3} + 185x^{2} + 229x + 81)}{432(x + 1)^{3}x^{3}}.

This example, taken from [GUl96], has Galois group F_{36} (in their notations). We search for invariants of degree m = 6.

7 This is for n ≥ 4. For n = 2 we need one row, and for n = 3, we need m + 1 rows.

The number of monomials to be evaluated is 8 at x = 0 and x = ∞, it is 16 at x = −1 and 28 at a regular point, so we work at x = 0. The heuristic at x = 0 with ν = 10 yields a 4-dimensional space of candidate invariants of degree 6 in 0.5 seconds. The complete algorithm (with the use of the above recurrence relations) then gives (as expected) a 2-dimensional space of invariants of degree 6 generated by

\frac{45}{4}X_{1}^{6} + 135X_{3}X_{2}y_{1,1}^{4} + 15X_{2}^{4}X_{1}^{2} + \frac{675}{2}y_{1,3}^{2}X_{2}^{2}X_{1}^{2}

− \frac{3645}{16}X_{3}^{4}X_{1}^{2} − 90X_{3}^{3}X_{2}^{3} − \frac{3645}{2}X_{3}^{5}X_{2} = −\frac{32805}{16}(x + 1)^{2}x^{5}

and

−\frac{9}{2}X_{1}^{6} + \frac{135}{2}X_{3}X_{2}y_{1,1}^{4} − 135X_{3}^{2}X_{2}^{2}X_{1}^{2} + \frac{3645}{4}X_{3}^{4}X_{1}^{2}

+ 6X_{3}X_{2}^{5} − 45X_{3}^{3}X_{2}^{3} \frac{729}{8}X_{3}^{5}X_{2} = \frac{6561}{8}(x + 1)^{2}x^{5}

The corresponding dual first integrals are

F_{1} = [−\frac{32805}{16}(x + 1)^{2}x^{5}, −\frac{10935}{32}x^{4}(x + 1)(7x + 5),
− \frac{1215}{128}x^{3} (−23 + 49x^{2} + 74x) , −\frac{243}{128}x^{3} (743 + 1463x^{2} + 2086x) ,
…, \frac{5}{339738624(x+1)^{9}x^{6}} (258162545x^{3}+8763967375x^{8}+3454336210x^{6}+4799496375x^{9}
+1442053125x^{10}+184528125x^{11}+8908159010x^{7}+430643385x^{2}−2112440530x^{5}
−2098573250x^{4}−44227701x−23914845),
\frac{5}{2717908992(x+1)^{10}x^{7}} (2491291800x^{3}+553584375x^{12}+30189672025x^{8}−9719448300x^{6}
+30901623000x^{9}+16603359750x^{10}+4702779000x^{11}+9196060400x^{7}−804729978x^{2}
−7242904720x^{5}+2863772665x^{4}−156676680x+71744535)]

and

F_{2} = [\frac{6561}{8}(x + 1)^{2}x^{5}, \frac{2187}{16}x^{4}(x + 1)(7x + 5),
\frac{243}{128}x^{3} (−55 + 89x^{2} + 130x) , \frac{243}{128}x^{3} (299 + 587x^{2} + 838x) ,
\frac{243x^{2}(205x^{3}+445x^{2}+83x−93)}{256(x+1)}, …,
\frac{1}{1358954496(x+1)^{9}x^{6}} (−1661573743x^{3}+3565866895x^{8}+1902296722x^{6}+3956079015x^{9}

+1442053125x^{10}+184528125x^{11}−1020492670x^{7}−1846340487x^{2}+12149631662x^{5}

+9005910718x^{4}+236373147x−23914845),

\frac{1}{10871635968(x+1)^{10}x^{7}} (−12566581020x^{3}+553584375x^{12}+41118775225x^{8}+41735111940x^{6}

+33723573540x^{9}+17335654830x^{10}+4791352500x^{11}+42844772456x^{7}+3414074670x^{2}
+13326450920x^{5}−21103738535x^{4}−68103180x+71744535]

To use the recurrence relations (∗), we have to compute the rows corresponding to monomials y^{i}(y′)^{6−i}, so that makes 7 rows (instead of 28). The bounds are more favorable for these 7 rows than for the other rows, and indeed the corresponding 7 entries of F_{1} and F_{2} are smaller expressions than most of the other entries. We printed the first 4 entries and the 2 last entries above. One sees that the last entries are significantly larger expressions. Precisely these large (hence: costly to compute) entries can be skipped in step 4 of algorithm 3, as these are the entries that are given by the recursion. This is the main reason why these relations are crucial for the efficiency. The computation time is 36.7 seconds and uses 1.5Mb of memory. We performed the same computation without using the recursion improvements, it took 263.5 seconds and used 2.5Mb. We then tried the first step of the standard method (computation of L^{◯ s 6}): this took 4587 seconds and more than 10Mb of memory.

5 Conclusion

We do not claim that our method is always better than the method via symmetric powers of operators. However, we have practical evidence that this method can handle much larger examples (and generally faster) than the previous one at our disposal.

To compute all the invariants of a given equation, we now face the following open problem: given L, can one bound the degrees of the generators of the invariant ring (when G is reductive)? As shown in [Com96], a solution to this problem would yield an algorithm for computing reductive unimodular Galois groups.

The method extends readily to systems: we then need formal solutions of systems (e.g. via cyclic vectors); but we lose the recurrence relations that enhance the algorithm, so the best there seems (surprisingly) to convert the system to an equation, apply the above algorithms, and then perform the correct change of variables to obtain the invariants of the original system.

We believe that the philosophy heuristic-checking is very suitable for computation. Information on the invariants can be obtained quickly by the heuristic and by modular arithmetic. If desired, this information can then be checked by algorithm 3. Furthermore algorithm 3 provides additional useful information, namely the dual first integrals corresponding to the invariants.

Applications of this algorithm are the computation of first integrals ([Wei95a]), the computation of differential relations satisfied by the solutions ([Com96]), the computation of algebraic and liouvillian solutions ([SUl93b,UWe94,SUl96]) and, more generally, to compute information on the Galois group. Extensions of the above techniques to other constructions on V (L) (and several applications) will be described in a subsequent paper.

References

[ABP95] Abramov S.A, Bronstein M & Petkošek M On polynomial solutions of linear operators Proceedings of ISSAC’95, ACM Press, 1995.

[Bar88] Barkatou A. Rational Newton Algorithm for computing formal solutions of linear differential equations Proceedings of ISSAC’88, ACM Press, 1988.

[Ber86] Bertrand D. Théorie de Galois différentielle Cours de DEA, Notes rédigées par R. Lardon, Université de Paris VI, 1986.

[BBe85] Bertrand D & Beukers F Équations différentielles linéaires et majorations de multiplicités, Ann. scient. Éc. Norm. Sup. 4ème série, t. 18 (1985), p. 181-192.

[Beu92] Beukers F Differential Galois theory In: From Number Theory to Physics (Ed: Waldschmidt, Moussa, Luck, Itzykson), Springer 1992.

[CLe55] Coddington E & Levinson NTheory of ordinary differential equations MacGraw-Hill, 1955.

[Com96] Compoint E Équations différentielles et relations algébriques Preprint (1995), Université de Paris 6.

[GUl96] Geiselmann W & Ulmer F Constructing a third order differential equation, preprint (1996), proceedings of the ”4-th Rhine Workshop on Computer Algebra”.

[HvP94] Hendriks P & Van der Put M Galois action on solutions of a differential equation, J. Symb. Comp., 1995.

[Hoe95a] van Hoeij M Formal solutions and factorization of differential operators with power series coefficients, University of Nijmegen Report nr. 9528.

[Hoe95b] van Hoeij M Factorization of differential operators with rational coefficients, University of Nijmegen Report nr. 9552.

[Lan92] Lang S. Algebra Third edition, Addison-Wesley, 1992.

[Lev89] Levelt A.H.M Differential Galois theory and tensor products Indagationes mathematicae (1989) + erratas.

[Mal79] Malgrange B Sur la réduction formelle des équations différentielles à singularités irrégulières, Manuscript (1979).

[MRa89] Martinet J & Ramis J.P Généralités sur la théorie de Galois différentielle In Computer Algebra and Differential Equations, Ed. E. Tournier, New York: Academic Press. (1990)

[Put95] van der Put M Singular complex differential equations: an introduction Nieuw Achief voor Wiskunde, 4de serie 13 No. 3 (1995) pp. 451-470.

[Sin89] Singer M.F An outline of differential Galois theory In Computer Algebra and Differential Equations, Ed. E. Tournier, New York: Academic Press. (1990).

[Sin93] Singer M.F Moduli of linear differential equations on the Rieman sphere Pac J. Math. Vol. 160, No. 2, (1993).

[Sin94] Singer M.F Testing reducibility of linear differential operators: a group theoretic perspective J. of Appl. Alg. in Eng. Comm. and Comp., 1995.

[SUl93a] Singer M.F & Ulmer F Galois groups for second and third order linear differential equations J.Symb.Comp 16 No. 1 (1993), pp 1-36.

[SUl93b] Singer M.F & Ulmer F Liouvillian and algebraic solutions of second and third order linear differential equations J.Symb.Comp 16 No. 1 (1993), pp 37-73.

[SUl94] Singer M.F & Ulmer F Necessary conditions for liouvillian solutions of (third order) linear differential equations J. of Appl. Alg. in Eng. Comm. and Comp. vol 6 No. 1 (1995) pp 1-22.

[SUl96] Singer M.F & Ulmer F Linear differential equations and products of linear forms. Preprint, presented at the MEGA’96 conference (Eindhoven, 6-8 june 1996).

[Tou87] Tournier E. Solutions formelles d’équations différentielles., Thèse d’Etat, Faculté des Sciences de Grenoble, 1987.

[UWe94] Ulmer F & Weil J.A Note on Kovacic’s algorithm Prepublication IRMAR 94-13, Rennes Juillet 94 (to appear in J. Symb. Comp.).

[Wei95a] Weil J.A First integrals and Darboux polynomials of homogeneous linear differential systems, Proceedings of AAECC 11 (Ed. M. Giusti & T. Mora), Lect. Notes in Comp. Sci. 948, Springer, 1995.

[Wei95b] Weil J.A Constantes et polynômes de Darboux en algèbre différentielle : application aux systèmes différentiels linéaires, PhD dissertation, École Polytechnique, 1995.