A Modular Algorithm to Compute the Exponential Solutions of a Linear Differential Operator

Thomas Cluzeau^{∗} LACO, Université de Limoges 123 avenue Albert Thomas 87060 Limoges, France

       thomas.cluzeau@unilim.fr

Mark van Hoeij^{†} Department of mathematics Florida State University Tallahassee, FL 32306, USA

                                               hoeij@math.fsu.edu

February 18, 2004

Abstract We present a new algorithm to compute exponential solutions of differential operators with rational function coefficients. We use a combination of local and modular computations, which allows to reduce the number of possibilities in the combinatorial part of the algorithm. We also show how unnecessarily large algebraic extensions of the constants can be avoided in the algorithm.

Introduction

An ordinary differential equation

y^{(n)} + a_{n−1}y^{(n−1)} + ⋯ + a_{0}y = 0  (1)

corresponds to a differential operator

∂^{n} + a_{n−1}∂^{n−1} + ⋯ + a_{0}∂^{0}  (2)

acting on y. Let C ⊂ \overline{ℚ} be some number field (we will allow more general C in Section 7.1). The topic of this paper is finding exponential solutions of (1) when the a_{i} are in C(x). By this we mean: finding all r ∈ \overline{ℚ}(x) for which y = exp(∫ r) is a solution of (1). This is equivalent to finding first order right hand factors ∂ − r with r ∈ \overline{ℚ}(x) of (2). Consequently, it can be viewed as part of the more general problem of factoring differential operators. Beke gave two factoring algorithms in [Bek94] (see also [PS03, 4.1 and 4.2.1]), one algorithm for first order factors (which we will refer to as Beke’s algorithm in this paper) and one algorithm (which uses first order factors) for higher order factors. Other factoring algorithms are the local to global method of [Hoe97a, Hoe97b] and the eigenring method ([Sin96], [PS03, 4.2.2]). Computing exponential solutions (computing first order factors) with these algorithms can be slow. It is not difficult to find examples where Maple’s expsols takes a long time, even though it uses an almost rational version of Beke’s algorithm, the local to global method, as well as the eigenring method. The goal in this paper is to give a faster algorithm. We combine modular methods with an almost rational version of Beke’s algorithm to obtain Algorithm ExpSolsInC which computes the exponential solutions defined over a given field C. Then we study the field problem, which is necessary for giving an efficient algorithm ExpSols for computing all exponential solutions.

∗Supported by a grant from the Région Limousin

†Supported by NFS grant 0098034

We sketch Beke’s method using the vocabulary of this paper. First, at each singularity, one computes a finite number (at most n) of objects called generalized exponents. Then, for each combination of generalized exponents (choose one generalized exponent at each singularity), the problem reduces to computing polynomial solutions (e.g., with [ABP95]) of some other differential operator. There are two main problems:

• a combinatorial problem: one has to choose a generalized exponent at each singularity which could lead to a large number (at most n^{m} where m is the number of singularities) of possible combinations.

• a field problem: singularities and generalized exponents are often defined over algebraic extensions of C, so some combinations can be defined over large algebraic extensions of C.

If neither problem occurs then Beke’s algorithm usually works well in practice. Thus we are mainly interested in situations where one or both of these two problems do occur. If we can discard combinations in advance then that would address the combinatorial problem, and in some cases even the field problem if the discarded combinations include those that are defined over large extensions of C. Our first main goal is to do this with computations modulo a prime number p.

The algebraic theory of ordinary differential equations in characteristic p is studied in [Hon81, Kat82, CC85, Put95]. For the problem of factoring differential operators with coefficients in \overline{𝔽}_{p}(x), van der Put showed in [Put95] (this is also recalled in [Put96, 2.1] and [PS03, 13.1]) that the relevant object is the p-curvature since it determines the differential module mod p and all factorizations of the operator. More recently, an algorithm to factor Ore polynomials with coefficients in 𝔽_{p}(x) is included in [GZ03]. This algorithm is based on the eigenring method adapted to characteristic p and generalizes the work of [Gie98] that tackles the problem of factoring differential operators over 𝔽_{q} where q = p^{l} for p prime and l ∈ ℕ^{∗}. An algorithm to factor differential systems over \overline{𝔽}_{p}(x), giving links between the p-curvature and the eigenring, is found in [Clu03].

The natural idea to lift modular factorizations to characteristic zero has been studied in [Put96, 4.3] and [PS03, 13.2.3]. In this paper we propose a different approach because there are many problems that prevent lifting from becoming an efficient and complete algorithm for characteristic zero. For example, there could be many factorizations in characteristic p, and it is not clear which of these should be lifted. Also, it is not clear how to do efficient Hensel lifting. Using more than one prime leads to the combinatorial problem of combining factors of different characteristics, and using large prime(s) makes p-curvature computations or factoring operators mod p expensive. We propose not to use factors of the operator mod p for problems in characteristic zero, but to use only the characteristic polynomial of the p-curvature.

To reduce the number of combinations to check in Beke’s algorithm, we prove links between eigenvalues in \overline{𝔽}_{p}(x^{p}) of the p-curvature and the reduction mod p of the generalized exponents. Algorithm CombMatchRoot determines which generalized exponents match a given eigenvalue. Algorithm ExpSolsInC uses this information to compute all exponential solutions defined over a field C (this C must be given in the input). Unfortunately, there are cases (see Section 5) where our mod p computations yield no improvements for the combinatorial problem.

The field problem refers to the problem that the usual Beke’s algorithm computes with a generalized exponent at every singularity simultaneously. So the degree of the field extension of C in Beke’s algorithm depends exponentially on m, the number of singularities (an extension of degree ≤ m! for the singularities followed by an extension of degree ≤ n^{m} for the generalized exponents). Of course this can easily choke the computation. Algorithm ExpSolsInC does not have this large field problem because the highest degree extensions of C it uses are those given by just one singularity. In fact, on complicated examples usually most of the work in ExpSolsInC consists of rational (i.e., over C) computations and not of computations in field extensions. However, for ExpSolsInC to find an exponential solution y, it needs in its input a field over which y can be defined. From now on, the field problem refers to the problem of finding such a field without constructing very large extensions that choke the computation. Our second main goal is solving this field problem.

Our final algorithm ExpSols first calls ExpSolsInC, then uses modular information to check if more exponential solutions could exist, and if so, it will first reduce the order as much as possible before tackling the field problem. After that, it calls the recursive FindASol algorithm which solves the field problem and returns an exponential solution defined over an extension, if such solution exists. We give several bounds (one of which is based on the p-curvature) for the degree of a field extension needed to find an exponential solution. Algorithm FindASol uses these bounds and makes an algebraic extension to decrease these bounds. It does this recursively and each time it calls ExpSolsInC to search for a solution over the new field. This continues until either a solution is found or until the bound becomes 1.

We summarize the two main new results in this paper:

• We use information mod p in several ways (in Section 5 for the combinatorial problem and Section 6 for improving bounds that are used in the field problem) to speed up computation of exponential solutions in characteristic 0 if possible.

To accomplish this there is a subtle but important difficulty to be overcome. There are many results known that hold for almost all p, but if such p can not easily be identified (an example is given in Section 7.3) then the algorithm can not use the result. The algorithm needs to know which p it can use. This restriction has numerous consequences, and the algorithm needs to be designed with this issue in mind. We will give a notion of good primes in Definition 3.3, and then only use results that hold for all good p. Another important issue is that the definition of good primes should discard as few primes as possible (for large p the p-curvature computation becomes slow, so only small p are useful for speeding up ExpSols).

• Our second main new result is Algorithm FindASol in Section 6.2, a recursive algorithm that constructs a field extension over which an exponential solution (if one exists) can be defined.

An efficiency advantage of Beke’s method over the local to global method in [Hoe97b] or the eigenring method in [Sin96] is that it directly computes first order factors. In contrast, the local to global method often needs to compute higher order factors before it can reach the first order factors, and the eigenring method first has to compute the eigenring. The greatest disadvantage of Beke’s method in comparison to the eigenring and local to global approaches is its use of a splitting field, and this is the issue that Algorithm FindASol resolves.

The paper is organized as follows. Section 1 contains the vocabulary needed to read this paper, in particular the notions of generalized exponents, reduction mod p, and p-curvature. Section 2 illustrates in a few examples how mod p computations can reduce the combinatorial problem. In Sections 3 and 4 we give links between the p-curvature and generalized exponents. Section 5 contains Algorithm CombMatchModp that reduces the number of combinations to check and Algorithm ExpSolsInC that computes all exponential solutions defined over a given field C. Section 6 addresses the field problem. Combining these algorithms in Section 7 we give Algorithm ExpSols that computes a basis of all exponential solutions, up to conjugation over C.

1 Preliminaries

In this section we briefly recall the concepts that are needed to read the next sections. The first subsection contains the necessary material in characteristic zero, the second subsection deals with characteristic p, and the third deals with reduction from characteristic 0 to characteristic p.

1.1 Differential operators in characteristic zero

Let \overline{ℚ} denote the algebraic closure of ℚ. We endow the fraction field \overline{ℚ}(x) of the ring \overline{ℚ}[x] of polynomials in the indeterminate x with the usual differentiation ′ := \frac{d}{dx}. The field of constants of the differential field (\overline{ℚ}(x),′) is then \overline{ℚ}. Consider the non-commutative ring \overline{ℚ}(x)[∂] of differential operators over the differential field (\overline{ℚ}(x),′). In this ring the multiplication is given by ∀ u ∈ \overline{ℚ}(x), ∂u = u∂ + u′. After multiplying if necessary by some element of \overline{ℚ}(x), a non-zero element L ∈ \overline{ℚ}(x)[∂] can be written as

L = a_{n}∂^{n} + a_{n−1}∂^{n−1} + … + a_{0}∂^{0}, a_{i} ∈ \overline{ℚ}[x], a_{n} ≠ 0 and gcd(a_{0}, …, a_{n}) = 1.  (3)

As usual, n will be called the order of the differential operator L.

1.1.1 Singularities and generalized exponents

We start with some classical definitions that can be found for example in [Inc26].

Definition 1.1. Let L ∈ \overline{ℚ}(x)[∂] be written as in (3). A point x_{0} ∈ \overline{ℚ} is said to be a singular point (or a singularity) of L if it is a root of a_{n}. Otherwise, it is said to be a regular (or ordinary) point of L.

In this paper we always consider the point at infinity as a singular point, whether it is actually singular or not. To each differential operator L ∈ \overline{ℚ}(x)[∂], we associate the subset S of \overline{ℚ} containing all finite singularities. The set of all singularities of L is then \overline{S} := S ∪ \{∞\}.

Let δ := x∂ and remark that \overline{ℚ}(x)[∂] = \overline{ℚ}(x)[δ] ⊂ \overline{ℚ}((x))[δ].

Definition 1.2. Let L be a non-zero differential operator in \overline{ℚ}((x))[δ]. L can be written ∑_{i=v}^{∞} x^{i}p_{i}(δ) for some polynomials p_{i} over \overline{ℚ} and p_{v} ≠ 0. Then the indicial equation of L is the polynomial p_{v}. Its roots in \overline{ℚ} are called the exponents of L at x = 0 (see [Inc26, 7.21] or [HW97, Lemma 11]).

The notion of exponents can be generalized in the sense of [HW97, 3.2]: there are several distinct ways to define generalized exponents. We choose here the algebraic definition of [Hoe97a, Hoe97b] that is easy to extend to characteristic p. We first give the definition at the point x = 0 and then extend it to all points. At a point x_{i} ∈ \overline{ℚ} ∪ \{∞\}, we define the local parameter t_{i} by

t_{i} := \{ \begin{matrix}x − x_{i} if x_{i} ∈ \overline{ℚ}, \\ 1/x if x_{i} = ∞. \end{matrix}

Definition 1.3. Let e ∈ \overline{ℚ}[x^{− \frac{1}{k}}] with k minimal in ℕ^{∗} and let L_{δ→δ+e} denote the differential operator L ∈ \overline{ℚ}(x)[δ] in which we have replaced δ by δ + e. Then e is said to be a generalized exponent of L at x = 0 if 0 is a root of the indicial equation of L_{δ→δ+e}. The number k is called the ramification (index) of e. The multiplicity of the root 0 in the indicial equation is called the multiplicity of e. If k = 1 then e is called an unramified generalized exponent. If k = 1 and L has no generalized exponent ẽ for which ẽ − e is a negative integer then e is called a ℤ-minimal generalized exponent.

Definition 1.4. Let e ∈ \overline{ℚ}[t_{i}^{− \frac{1}{k}}]. Then e is said to be a generalized exponent of L at x = x_{i} if e_{t_{i}→x} (e in which we have replaced t_{i} by x) is a generalized exponent of L_{t_{i}→x} at x = 0. Here L_{t_{i}→x} refers to the automorphism of \overline{ℚ}(x)[∂] given by x ↦ x+x_{i} at a finite point and x ↦ 1/x at infinity. This leads to ∂ ↦ ∂ at a finite point and ∂ ↦ −x^{2}∂ at infinity.

Generalized exponents are useful for computing exponential solutions, see Lemma 1.8 and Remark 1.10. At each point x_{i}, there are exactly n generalized exponents counted with multiplicity. The point x_{i} is called a regular singular point of L when all generalized exponents are constants, in which case they are just the usual exponents. At a regular point x_{i} the (generalized) exponents are 0, 1, …, n − 1. For further details on generalized exponents see [Hoe97a, Hoe97b, HW97].

1.1.2 Exponential solutions and first order right hand factors

Definition 1.5. Let L ∈ \overline{ℚ}(x)[∂] and let y ≠ 0 belong to a differential field extension of \overline{ℚ}(x) with \overline{ℚ} as field of constants. Assume that L(y) = 0. Then, y is said to be an exponential solution of L if y′/y ∈ \overline{ℚ}(x). It is said to be a radical solution if there exists a positive integer m such that y^{m} ∈ \overline{ℚ}(x) and finally, it is said to be a rational (respectively polynomial) solution if y ∈ \overline{ℚ}(x) (respectively y ∈ \overline{ℚ}[x]).

An exponential solution y of L corresponds to the first order right hand factor ∂ − y′/y of L.

Remark 1.6. In this paper it will not be necessary to distinguish an exponential function y from a scalar multiple of y. Hence we may represent y by its logarithmic derivative y′/y ∈ \overline{ℚ}(x) or by its minimal annihilating operator ∂ − y′/y.

Definition 1.7. Let L ∈ \overline{ℚ}(x)[∂] be of order n > 1. R ∈ \overline{ℚ}(x)[∂] of order 0 < m < n is said to be a right hand factor of L if there exists l ∈ \overline{ℚ}(x)[∂] such that L = lR.

Lemma 1.8 below (analogous to [Bou00, Proposition 8]) shows that the notions of generalized exponent, exponential solution and first order right hand factor are closely related.

Lemma 1.8. Let r ∈ \overline{ℚ}(x). If ∂ − r is a first order right hand factor of L ∈ \overline{ℚ}(x)[∂], then for all x_{i} in \overline{S}, there exists a generalized exponent e_{x_{i}} ∈ \overline{ℚ}[t_{i}^{−1}] of L at x_{i} such that

r = S + \frac{Q′}{Q} where S = ∑_{x_{i}∈S} \frac{e_{x_{i}}}{t_{i}} − t_{∞}e_{∞}^{∗},

and (Fuchs’ relation)

deg(Q) + ∑_{x_{i}∈\overline{S}} Const(e_{x_{i}}) = 0, in particular − ∑_{x_{i}∈\overline{S}} Const(e_{x_{i}}) ∈ ℕ  (4)

where Q ∈ \overline{ℚ}[x], the Const(e_{x_{i}}) are the constant terms of the e_{x_{i}} and e_{∞}^{∗} := e_{∞} − Const(e_{∞}).

Proof. If L = ∑ a_{i}∂^{i} then denote L_{∂→∂+S} = ∑ a_{i}(∂ + S)^{i}. Let e_{x_{i}} be the generalized exponent of ∂ − r at x_{i}. Since ∂ − r is a right hand factor of L, this is also a generalized exponent of L (see [Hoe97b, Lemma 3.2]). Now L_{∂→∂+S} has right hand factor ∂ − r̃ where r̃ = r − S. Let S′ be the set of all finite singularities of ∂ − r and let Δ = S′ \ S. So x_{i} ∈ Δ if and only if x_{i} is a regular point of L (hence e_{x_{i}} ∈ ℕ) but a singular point (hence e_{x_{i}} ≠ 0) of ∂ − r. Now r = ∑_{x_{i}} \frac{e_{x_{i}}}{t_{i}} − t_{∞}e_{∞}^{∗} if we take the sum over all finite singularities of ∂ − r. Hence r̃ = ∑_{x_{i}∈Δ} \frac{e_{x_{i}}}{t_{i}}. Then ∂ − r̃ has a polynomial solution Q = ∏_{x_{i}∈Δ} t_{i}^{e_{x_{i}}}. Fuchs’ relation (which applied to ∂ − r is just a reformulation of the fact that the sum of the residues of r is zero) says that if we sum Const(e_{x_{i}}) over all singularities of ∂ −r then the result is zero, ∑_{x_{i}∈\overline{S}∪Δ} Const(e_{x_{i}}) = 0. Then Equation (4) follows from the observation that deg(Q) = ∑_{x_{i}∈Δ} e_{x_{i}}.

Remark 1.9. In our algorithm we only need to compute unramified generalized exponents of L, those are the only ones that can be relevant for our algorithm.

If at a singularity x_{i} ∈ S there are two generalized exponents e_{x_{i},1} and e_{x_{i},2} whose difference d = e_{x_{i},1} −e_{x_{i},2} is a positive integer, then e_{x_{i},1} is not relevant in the sense that the lemma stays true if we never use e_{x_{i},1} (because replacing e_{x_{i},1} and Q by e_{x_{i},2} and t_{i}^{d}Q does not change S + Q′/Q). So, in our algorithm, we will only use ℤ-minimal generalized exponents to construct r.

The r from the previous lemma matches the exponential solution

exp(∫ S) Q.

If all e_{x_{i}} are rational numbers then this exponential solution can be written as

Q ∏_{x_{i}∈S} t_{i}^{e_{x_{i}}}  (5)

which is a radical solution. So the radical solutions are those exponential solutions for which all associated generalized exponents are rational numbers.

Remark 1.10. If we only know the generalized exponents defining some exponential solution modulo rational numbers (respectively modulo integers), then finding such exponential solution reduces to finding radical solutions (respectively rational solutions) of L_{∂→∂+S}.

1.2 Differential operators in characteristic p

Let \overline{𝔽}_{p} denote the algebraic closure of the finite field 𝔽_{p}. We define the ring \overline{𝔽}_{p}(x)[∂] of differential operators in the same way as in characteristic zero. The main difference is that the field of constants of (\overline{𝔽}_{p}(x),′) is \overline{𝔽}_{p}(x^{p}) (see [Put01] or [GZ03, Lemma 3.3]).

1.2.1 The p-curvature and the map y ↦ y^{(p−1)} + y^{p}

The central tool for studying differential operators in characteristic p is the p-curvature; we recall a definition and we refer to [Kat82, CC85] and [PS03, 13.2.2] for more information. The differential field \overline{𝔽}_{p}(x) is a finite dimensional vector space over its field of constants \overline{𝔽}_{p}(x^{p}):

\overline{𝔽}_{p}(x) = ⨁_{i=0}^{p−1} \overline{𝔽}_{p}(x^{p})x^{i}.

To any differential operator L ∈ \overline{𝔽}_{p}(x)[∂], one can associate the differential module

M_{L} := \overline{𝔽}_{p}(x)[∂]/\overline{𝔽}_{p}(x)[∂]L.

This module is equipped with a \overline{𝔽}_{p}(x^{p})-linear map ˜∂ (coming from the left multiplication by ∂ on \overline{𝔽}_{p}(x)[∂]) satisfying

˜∂fm = f′m + f ˜∂m for all m ∈ M_{L}, f ∈ \overline{𝔽}_{p}(x).

Definition 1.11. Let L ∈ \overline{𝔽}_{p}(x)[∂]. The p-curvature of L is the \overline{𝔽}_{p}(x)-linear map ˜∂^{p} acting on the differential module M_{L} associated to L.

An algorithm to compute the p-curvature matrix is given in [PS03, 13.2.2]. The map

\begin{aligned}τ : & \overline{𝔽}_{p}(x) & → & \overline{𝔽}_{p}(x^{p}), \\ & y & ↦ & y^{(p−1)} + y^{p} \end{aligned}

is very useful for order one right hand factors in characteristic p. Details can be found in [[PS03], 13.2.1] and [Put95]; we only recall a lemma that we will often use:

Lemma 1.12 ([Put95], Lemma 1.4.2). The map τ : \overline{𝔽}_{p}(x) → \overline{𝔽}_{p}(x^{p}) is a surjective additive map with kernel \{ \frac{z′}{z} | z ∈ \overline{𝔽}_{p}(x)^{∗} \}.

In [Put95], van der Put proves that for the problem of factoring a differential operator in characteristic p (or in an equivalent way its associated differential module M_{L}), the p-curvature is the relevant object since its characteristic properties give all possible factorizations of the operator.

We will note χ_{p}(L) the characteristic polynomial of the p-curvature of L ∈ \overline{𝔽}_{p}(x)[∂]. As we are only interested in first order factors, the only modular information that we will use are the roots in \overline{𝔽}_{p}(x^{p}) (with multiplicities) of χ_{p}(L). We denote R(χ_{p}(L)) as the set of roots of χ_{p}(L) in \overline{𝔽}_{p}(x^{p}).

Lemma 1.13. Let L ∈ \overline{𝔽}_{p}(x)[∂]. (a). If L = L_{1}L_{2} then χ_{p}(L) = χ_{p}(L_{1})χ_{p}(L_{2}) with deg(χ_{p}(L_{i})) = order(L_{i}). (b). The map μ from the set of all first order monic right hand factors of L in \overline{𝔽}_{p}(x)[∂] to R(χ_{p}(L)) given by μ(∂ − r) := τ(r) is well defined and surjective.

Proof. The claim (a) follows from the classification (the equivalence of categories) of [Put95] and [PS03, 13.1]. That the map μ is well defined follows from the first claim and the fact that τ(r) is the p-curvature of ∂ − r (see [Put96, Lemma 2.2]). Surjectivity follows from the classification as well.

1.3 Reduction modulo p and factorization

Let v_{p} : ℚ → ℤ ∪ {∞} denote the standard p-adic valuation on ℚ. There are infinitely many ways to extend this valuation to a valuation v_{p} : \overline{ℚ} → ℚ ∪ \{∞\}. We assume that one such choice is made. Let R_{p} = \{a ∈ \overline{ℚ} | v_{p}(a) ≥ 0\}, and consider the maximal ideal I = \{a ∈ R_{p} | v_{p}(a) > 0\}. Now choose an isomorphism R_{p}/I ≅ \overline{𝔽}_{p}.

Definition 1.14. Having fixed the above choices, we denote the image of a ∈ R_{p} in \overline{𝔽}_{p} by a[p]. This is called the reduction of a mod p.

Note that a[p], the reduction of a mod p, depends not only on a and p but depends on the choice of v_{p} as well. However, for compactness of notation we will only mention p even when we refer to both p and v_{p}. If a ∈ \overline{ℚ} is not in R_{p} then a[p] is not defined. In our algorithm, we will apply the reduction map [p] : R_{p} → \overline{𝔽}_{p} to only finitely many algebraic numbers that are known in advance. Thus, it is not difficult to choose a prime p and valuation v_{p} in such a way that a[p] is always defined for all numbers a for which the reduction map [p] is used. The a’s we reduce mod p are finite in number and thus they are elements of a finite extension K of ℚ. There are only finitely many possible choices for v_{p} on K. We will assume that such p and v_{p} are chosen, so that all instances of the notation a[p] are defined. For more details on this choice see Section 5.3.

In practice, reducing an algebraic number a mod p is done as follows: Let a_{1}, …, a_{n} be those algebraic numbers that have previously been reduced mod p, but take only those whose reduction mod p may appear in the same expression as a[p] in the algorithm (if the reductions of a_{1}, …, a_{n} and some other numbers b_{1}, …, b_{m} are used completely independently, then the algorithm will still be valid if it computes the reductions of these two sets of numbers independently, using one valuation on ℚ(a_{1}, …, a_{n}) and an independently chosen valuation on ℚ(b_{1}, …, b_{m}). This way the larger field ℚ(a_{1}, …, a_{n}, b_{1}, …, b_{m}) does not need to be constructed). Then determine the minimal polynomial of a over K := ℚ(a_{1}, …, a_{n}). Reduce this minimal polynomial mod p, then the coefficients are in a field we will denote by K[p], which is defined as 𝔽_{p}(a_{1}[p], …, a_{n}[p]). Factor this polynomial over K[p], and choose one of the irreducible factors. Then determine a field extension of K[p] in which this irreducible factor has a root, and let a[p] be that root.

We will apply the reduction map [p] to the singularities of L, which are in \overline{ℚ} ∪ \{∞\}. If a = ∞ then we define a[p] as ∞. Two distinct singularities of L can have the same image under [p], in which case we say that p is not a good prime (see condition (C3) in Definition 3.3 in Section 3). Note that this definition of good prime depends not only on p, but depends on v_{p} as well (recall that when we refer to both p and v_{p} we will only mention p for brevity).

If a is some object defined over \overline{ℚ} (like L, a local parameter t_{i}, or a generalized exponent e_{x_{i}}) we can define a[p] by applying the reduction map [p] to all coefficients.

We show now that without too many assumptions on p a factorization of a differential operator in characteristic zero can be reduced mod p. This result can be found in [Put96].

Proposition 1.15. Let L ∈ \overline{ℚ}(x)[∂] and let p be such that L can be reduced mod p. If L = L_{1}L_{2} then, after possibly replacing L_{1} and L_{2} by cL_{1} and \frac{1}{c}L_{2} for some constant c ∈ \overline{ℚ}, we have L[p] = L_{1}[p] L_{2}[p].

We allow to replace L_{1} and L_{2} by cL_{1} and \frac{1}{c}L_{2} for some constant c to avoid pathologies of the form L_{1} = \frac{1}{p}∂ and L_{2} = p∂.

Proof. The only thing to prove is that cL_{1} and \frac{1}{c}L_{2} can be reduced mod p for some constant c. Any p-adic valuation v_{p} on \overline{ℚ} can be extended to \overline{ℚ}[x] (this statement is a form of Gauss’ lemma. To define v_{p}(f) for a non-zero polynomial one takes the smallest valuation of the coefficients). The valuation can then be extended to the fraction field \overline{ℚ}(x), and then be extended to \overline{ℚ}(x)[∂] (again, take the minimal valuation of the coefficients). Then L[p] is defined if and only if v_{p}(L) ≥ 0. Let c be a constant with valuation −v_{p}(L_{1}), so v_{p}(cL_{1}) = 0, hence cL_{1} can be reduced mod p. Then v_{p}(\frac{1}{c}L_{2}) = v_{p}(cL_{1}) + v_{p}(\frac{1}{c}L_{2}) = v_{p}(cL_{1} \frac{1}{c}L_{2}) = v_{p}(L_{1}L_{2}) = v_{p}(L) ≥ 0, so \frac{1}{c}L_{2} can be reduced mod p as well.

One obtains the following two criteria (these can also be easily deduced from [Put95]):

Corollary 1.16. Let L ∈ \overline{ℚ}(x)[∂] and let p be a prime such that L can be reduced mod p. Assume further that the order n of L does not drop after the reduction. If χ_{p}(L) is irreducible over \overline{𝔽}_{p}(x^{p}) then L is irreducible over \overline{ℚ}(x). If χ_{p}(L) has no roots in \overline{𝔽}_{p}(x^{p}) then L has no exponential solutions.

Proof. It follows directly from the previous proposition, [Clu03, Theorem 4.1-Step 2] and Lemma 1.13(b).

As an example, let L = (x^{2} + x + 8)∂^{2} + (−x^{8} + x + 6)∂ + 1 and p = 3. Then χ_{p}(L) has no roots and this is a very fast way to show that L has no exponential solutions. Without modular methods this would have taken much longer because one of the degree bounds (one ends up searching for polynomial solutions of some other operators) turns out to be very high in this example.

In the corollary, the hypothesis that the order does not drop can not be omitted. For example, the operator (p∂ + 1)∂ is reducible in characteristic 0 but its reduction L[p] = ∂ is irreducible.

2 Examples

We will illustrate our modular improvements for the combinatorial problem (see the introduction) with a few examples. These improvements rely on Proposition 1.15 and Lemma 1.13(b) in the following way. A first order factor in characteristic zero leads to a first order factor in characteristic p by Proposition 1.15, which in turn corresponds to a root of χ_{p}(L) by Lemma 1.13(b).

2.1 Example 1

Consider the following differential operator:

L = ∂^{3} − \frac{2x^{2} − x + 4}{2x^{2}}∂^{2} − \frac{3x^{3} − 4x^{2} − 3x − 2}{2x^{4}}∂ + \frac{2x^{3} − 3x − 2}{2x^{4}}.

The singularities are 0 and ∞. Generalized exponents are 0, \frac{5}{2} + \frac{1}{x}, 2+ \frac{1}{x} at 0 and −t_{∞}^{−1}, −t_{∞}^{− \frac{1}{2}} , t_{∞}^{− \frac{1}{2}} at ∞ (recall that t_{∞} = \frac{1}{x}). The usual Beke’s algorithm will have three combinations to check since we only have one unramified generalized exponent at ∞, and three at 0. We will use this example to illustrate our modular method, even though our method generally does not speed up trivial examples (L is trivial for Beke’s algorithm because the number of combinations is very small and there is no field problem).

Let p = 3 and note c = x^{p}. The characteristic polynomial χ_{p}(L) of the p-curvature matrix is (λ^{2} + λ/c^{2} + 2/c + 1/c^{4})(λ + 2) and has only the root λ = 1. We will see in Section 5.2 that this root excludes the choices \frac{5}{2} + \frac{1}{x} and 2 + \frac{1}{x} at 0; consequently, there remains only one combination to check.

For each of the 3 combinations, Beke’s algorithm will check Fuchs’ relation (Equation (4) in Lemma 1.8). This check will be very fast because since there is no field problem in this example, we do not have to construct a complicated number field in order to do this check. One of the three combinations passes this check, and it yields a solution. In this example, the modular method would not provide a speedup because it discarded only two combinations (which otherwise would have been discarded anyway by Fuchs’ relation). In the next example, our modular method discards more combinations, and moreover, the discarded combinations include all “hard” combinations, i.e., those not invariant under the Galois group of \overline{ℚ} over ℚ.

2.2 Example 2

Let

L = 9(x^{3} − 2)^{5}∂^{3} + (x^{3} − 2)(2x^{10} − 12x^{7} + 108x^{5} + 24x^{4} − 216x^{2} − 16x − 9)∂
−2x(190x^{6} − 274x^{3} − 27x − 212).

There are three finite singularities; the roots of x^{3} − 2, with generalized exponents 2, \frac{1}{36} \frac{α^{2}}{(x−α)} − \frac{1}{18}α

and −\frac{1}{36} \frac{α^{2}}{(x−α)} + 4 + \frac{1}{18}α where α := RootOf(x^{3} − 2). At infinity, the generalized exponents are 0, −\frac{5}{3}, −\frac{4}{3}. So the usual Beke’s algorithm will have 81 combinations to try. Let p = 5 and look at the reduction mod p of the singularities and generalized exponents. First, we factor x^{3} − 2 = (x + 2)(x^{2} + 3x + 4) mod p, so α reduces to 3 as well as to β := RootOf(x^{2} + 3x + 4) in characteristic p. After reduction mod p the situation is then as follows: we have the singularity 3 with generalized exponents 2, \frac{4}{x+2} − 1, −\frac{4}{x+2} and the singularities β with generalized exponents 2, \frac{β^{2}}{(x−β)} +3β, −\frac{β^{2}}{(x−β)} +4+2β. At infinity the generalized exponents become 0 (with multiplicity 2)

and 2. Now, we compute χ_{p}(L) and see that its only root is λ = s where s := \frac{4c^{4}+2c+1}{c^{6}+c^{3}+4} and c = x^{p}. Our method consists then in finding which generalized exponents reduced mod p could match such a root s. Here, the result is that s can only correspond to the choice \frac{4}{x+2} − 1 at 3, \frac{β^{2}}{(x−β)} + 3β at β and any choice at infinity. This means that in characteristic zero, a possible exponential solution must involve the generalized exponent \frac{1}{36} \frac{α^{2}}{(x−α)} − \frac{1}{18}α at α. This reduces the number of possible combinations from 81 to 3. Moreover, we have eliminated all “hard” combinations, i.e., those for which the sum S in Lemma 1.8 would not be in ℚ(x).

If two generalized exponents differ by a rational number, then our mod p computation can not detect which one to take, which is an important limitation. This is why three combinations remained after the mod p computation even though there was only one root of χ_{p}(L). Of the three remaining combinations only one satisfies Fuchs’ relation, and this combination yields an exponential solution.

2.3 Example 3

Let p be a prime number and consider the operator L = x^{4}∂^{2} − px∂ − 2x − 1. Then p is a good prime in the sense of Definition 3.3 below. However, despite the fact that L[p] exists and can even be factored, neither of the two generalized exponents 3 + 1/p + p/x^{2} and −1/p of L at x = 0 can be reduced mod p. So both can be discarded by Lemma 3.4 and hence L has no exponential solutions.

3 Roots of χ_{p}(L) and τ-reduced terms

Definition 3.1. If r is an object in characteristic zero we set τ(r) := τ(r [p]). Let e_{x_{i}} ∈ \overline{ℚ}[t_{i}^{−1}] be an unramified generalized exponent of L ∈ \overline{ℚ}(x)[∂] at a point x_{i} ∈ \overline{ℚ} ∪ \{∞\}. Let p be such that x_{i} and e_{x_{i}} can both be reduced mod p. The τ-reduced term corresponding to e_{x_{i}} is then defined as

τ_{red}(e_{x_{i}}) := \{ \begin{matrix}τ(\frac{e_{x_{i}}}{t_{i}}) if x_{i} ∈ \overline{ℚ}, \\ −τ(t_{∞}e_{∞}^{∗}) if x_{i} = ∞ \end{matrix}

where e_{∞}^{∗} denotes e_{∞} without its constant term.

Remark 3.2. If e ∈ ℚ then e[p] ∈ 𝔽_{p} and Lemma 1.12 implies τ_{red}(e) = 0.

Definition 3.3. A prime p is a good prime for L if the following conditions hold: (C1)- The coefficients of L can be reduced mod p and the order of L does not decrease after reduction mod p, i.e., a_{n}[p] ≠ 0. (C2)- All singularities can be reduced mod p. (C3)- Distinct singularities stay distinct after reduction mod p.

If L is written in the form (3), then Condition (C2) simply means “deg_{x}(a_{n}[p]) = deg_{x}(a_{n})”. And (C3) means that the largest square-free factor of a_{n} remains square-free mod p. Checking if a prime p is good is fairly easy in practice. All three conditions are useful in our algorithm, for example, without (C3) we would have to replace Algorithm CombMatchRoot in Section 5.2 by a more complicated algorithm. The lemma below uses (C3) as well. Nevertheless, if p is not a good prime, then it is often still possible to obtain some useful information from a computation mod p, but we will not detail this.

Lemma 3.4. Let L ∈ \overline{ℚ}(x)[∂] and p be a good prime for L. If ∂ − r with r ∈ \overline{ℚ}(x) is a monic first order right hand factor of L and if e_{x_{i}} ∈ \overline{ℚ}[t_{i}^{−1}] is the generalized exponent of ∂ − r at x_{i}, then e_{x_{i}} can be reduced mod p.

Proof. From Lemma 1.8, there exists a generalized exponent e_{x_{i}} ∈ \overline{ℚ}[t_{i}^{−1}] of L at each singularity x_{i} such that

r = ∑_{x_{i}∈S} \frac{e_{x_{i}}}{t_{i}} − t_{∞}e_{∞}^{∗} + \frac{Q′}{Q} ,  (6)

where Q ∈ \overline{ℚ}[x] and e_{∞}^{∗} is e_{∞} without its constant term. Since p satisfies (C1), r can be reduced mod p (see Proposition 1.15). Now Q is defined up to a multiplicative constant (which we can choose in such a way that Q[p] is defined and is not zero), so Q′/Q can be reduced mod p and hence the remaining part ∑_{x_{i}∈S} e_{x_{i}}/t_{i} − t_{∞}e_{∞}^{∗} can be reduced mod p as well. Conditions (C2)+(C3) in Definition 3.3 imply that we can recover e_{∞}^{∗}[p] and e_{x_{i}}[p] for x_{i} ∈ S from (∑_{x_{i}∈S} e_{x_{i}}/t_{i} − t_{∞}e_{∞}^{∗})[p] and therefore e_{x_{i}}[p] is well defined. Then the same is true for e_{∞}[p] by Fuchs’ relation.

We now give two new results linking the τ-reduced terms and the eigenvalues of the p-curvature in \overline{𝔽}_{p}(x^{p}), i.e., the elements of R(χ_{p}(L)). First another definition:

Definition 3.5. Let L ∈ \overline{ℚ}(x)[∂]. We say that s ∈ R(χ_{p}(L)) comes from characteristic zero if there exists r ∈ \overline{ℚ}(x) such that ∂ − r is a right hand factor of L in characteristic zero and s = τ(r).

Proposition 3.6. Let L ∈ \overline{ℚ}(x)[∂] and let p be a good prime for L. If s ∈ R(χ_{p}(L)) comes from characteristic zero, then ∀ x_{i} ∈ \overline{S}, there exists a generalized exponent e_{x_{i}} at x_{i} such that:

s = ∑_{x_{i}∈\overline{S}} τ_{red}(e_{x_{i}}).

Proof. By definition, s ∈ R(χ_{p}(L)) comes from characteristic zero means that s = τ(r) for some r ∈ \overline{ℚ}(x) where ∂ − r is a right hand factor of L. From Lemma 1.8, there exists then a generalized exponent e_{x_{i}} ∈ \overline{ℚ}[t_{i}^{−1}] of L at each singularity x_{i} such that r can be written in the form (6) where Q ∈ \overline{ℚ}[x] and e_{∞}^{∗} is e_{∞} without its constant term. Since p is a good prime for L and Q is defined up to a multiplicative constant (which we can choose in such a way that Q[p] is defined and is not zero), the reduction mod p of each term in Equation (6) is well defined (see also Lemma 3.4 above). Now using that s = τ(r), the result follows directly from Lemma 1.12 and the definition of the τ_{red}(e_{x_{i}}).

Proposition 3.7. For almost all primes p the condition that s comes from characteristic 0 may be omitted.

Proof. As we will not need this result for our algorithm, we will only sketch a proof. First we note that generalized exponents can be defined in characteristic p. Definitions 1.3 and 1.4 are completely algebraic: the operations required are additions and multiplications in \overline{ℚ}(x), tests to zero in \overline{ℚ} and computing the roots of a polynomial in \overline{ℚ}[x]. Obviously these operations also make sense in characteristic p so the definitions are easily translated to characteristic p provided that the ramification k is not divisible by p (this is not an issue because we only consider k = 1 in this paper). The generalized exponents are here defined mod p, which is coherent with the fact that we want t_{i}^{p} to have generalized exponent 0 at x_{i} (because t_{i}^{p} is in the field of constants) as well as generalized exponent p at x_{i} because it vanishes with multiplicity p at x_{i}. Note that the notions of singularity, indicial equation and exponents in characteristic p have already been introduced in [Hon81] and [CC85].

Lemma 1.8 can also be adapted to characteristic p (the lemma is just a statement about the relation between a rational function and its local expansions). Then one can also give the analogue of Proposition 3.6 in characteristic p, with L ∈ \overline{𝔽}_{p}(x)[∂], with the x_{i} in \overline{𝔽}_{p} ∪ \{∞\}, the unramified generalized exponents e_{x_{i}} in \overline{𝔽}_{p}[t_{i}^{−1}], but now s can be any root of χ_{p}(L).

Let L ∈ \overline{ℚ}(x)[∂], let p be a good prime, and let s be any root of χ_{p}(L), not necessarily coming from characteristic 0. Applying the mod p analogue of Proposition 3.6 to L[p] we see that even if s does not come from characteristic 0, the statement in Proposition 3.6 still holds if we replace the generalized exponents e_{x_{i}} of L by the generalized exponents of L[p]. This completes the proof because for all but finitely many good primes p, reduction mod p will induce a bijection between the (unramified) generalized exponents of L and those of L[p].

One obtains the following corollary which is related to the converse of [Hon81, Theorem 2] and [CC85, Corollary 1.4]. Indeed, their statement (attributed to Katz) is: if L[p] has “sufficiently many solutions in a weak sense”, which is equivalent to “L has nilpotent p-curvature”, for almost all p, then all exponents of L are rational numbers.

Corollary 3.8. Let L ∈ \overline{ℚ}(x)[∂]. If all generalized exponents of L are rational numbers, then for almost all primes, R(χ_{p}(L)) ⊆ {0}.

Proof. If e_{x_{i}} ∈ ℚ then τ_{red}(e_{x_{i}}) = 0 and the result follows from Proposition 3.7.

4 Exponential solutions and roots of χ_{p}(L)

4.1 Classes of exponential solutions

Definition 4.1. The class \overline{e_{x_{i}}} of a generalized exponent e_{x_{i}} is \overline{e_{x_{i}}} := e_{x_{i}} mod ℚ, the equivalence class of e_{x_{i}} modulo ℚ.

Two generalized exponents e_{x_{i}} and e′_{x_{i}} that have the same class \overline{e_{x_{i}}} and that can both be reduced mod p satisfy τ_{red}(e_{x_{i}}) = τ_{red}(e′_{x_{i}}) by Remark 3.2. So we can define τ_{red}(\overline{e_{x_{i}}}) := τ_{red}(e_{x_{i}}).

Definition 4.2. Let y be an exponential solution of L, let ∂ − r be the associated first order right hand factor of L (so r = y′/y) and let e_{x_{i}} be the associated generalized exponents in the sense of Lemma 1.8. Then the class Cl(y) of y is defined as the image of y′/y in \overline{ℚ}(x)/∼ where a ∼ b ⇔ a − b = \frac{g′}{g} for some radical function g ≠ 0. Viewing \overline{e_{x_{i}}}/t_{i} as elements of \overline{ℚ}(x)/∼ we can write

Cl(y) = \frac{y′}{y} mod ∼ = ∑_{x_{i}∈S} \frac{\overline{e_{x_{i}}}}{t_{i}} − t_{∞}e_{∞}^{∗}.

Two exponential solutions have the same class if and only if their quotient is a radical function, if and only if at every point, the associated generalized exponents have the same class. Now, to L ∈ \overline{ℚ}(x)[∂] we associate the class-set

Cl(L) := \{ Cl(y) ∈ \overline{ℚ}(x)/∼ | y exponential solution of L \}.

4.2 Links between Cl(L) and R(χ_{p}(L))

Lemma 4.3. Let L ∈ \overline{ℚ}(x)[∂] and p be a good prime for L. There is a map

Ψ : Cl(L) → R(χ_{p}(L)),
Cl(y) ↦ τ(\frac{y′}{y}) = ∑_{x_{i}∈\overline{S}} τ_{red}(\overline{e_{x_{i}}}).

Proof. Let Cl(y) ∈ Cl(L). An exponential solution y corresponds to a first order right hand factor ∂ −r of L in characteristic zero where r = y′/y. Since p is a good prime for L, this factorization can be reduced mod p (see Proposition 1.15) and we have thus a first order right hand factor ∂ − r[p] of L[p]. Now from Lemma 1.13, τ(r) ∈ R(χ_{p}(L)). The lemma is then clear using Proposition 3.6 since τ(r) is an element of R(χ_{p}(L)) that comes from characteristic zero.

Given L, there are two natural questions about Ψ: (Q1)- Is Ψ surjective for almost all good primes? (Q2)- Does there exist a good prime p such that Ψ is injective?

Proposition 4.4. The answer to questions (Q1) and (Q2) can be no.

Proof. We give one example for each question. For (Q1), if we take a differential operator that has a basis of algebraic (but no exponential) solutions, then Cl(L) = ∅. The proved part of Grothendieck’s conjecture (see [Kat82], [Cha01] or [Put01] for statements and histories of this conjecture) says that for almost all primes the p-curvature is zero so that R(χ_{p}(L)) = {0}. For (Q2), take an operator L with the following 8 exponential solutions:

x^{e_{0}}(x − 1)^{e_{1}}(x − 2)^{e_{2}}

for e_{0} ∈ \{0, \sqrt{2}\}, e_{1} ∈ \{0, \sqrt{3}\} and e_{2} ∈ \{0, \sqrt{6}\}. Modulo every prime p, at least one of the three square roots becomes an element of 𝔽_{p}. Suppose for example that \sqrt{6}[p] ∈ 𝔽_{p}. Consider then the exponential solutions y_{1} and y_{2} corresponding to choices respectively (0, 0, 0) and (0, 0, \sqrt{6}) for (e_{0}, e_{1}, e_{2}). Then Cl(y_{1}) ≠ Cl(y_{2}) but they have the same image under Ψ.

The previous proposition says that we can not always choose a prime number for which each root of χ_{p}(L) corresponds to at most one class of exponential solutions. The example illustrates the cause of this, and one can show the following (note that we do not use this in our algorithm):

Proposition 4.5. Let L ∈ \overline{ℚ}(x)[∂]. If there exists a good prime p for L for which:

(∗) ∀ x_{i} ∈ \overline{S}, e_{x_{i},2} − e_{x_{i},1} ∉ ℚ ⇒ e_{x_{i},2}[p] − e_{x_{i},1}[p] ∉ 𝔽_{p}

then, for that prime, Ψ is injective.

Proof. Suppose that Cl(y_{1}) = ∑_{x_{i}∈S} \frac{\overline{e_{x_{i},1}}}{t_{i}} − t_{∞}e_{∞,1}^{∗} and Cl(y_{2}) = ∑_{x_{i}∈S} \frac{\overline{e_{x_{i},2}}}{t_{i}} − t_{∞}e_{∞,2}^{∗} are two distinct elements of Cl(L) and suppose that Ψ(Cl(y_{1})) = Ψ(Cl(y_{2})).

Ψ(Cl(y_{1})) = Ψ(Cl(y_{2})) ⇔  ∑_{x_{i}∈\overline{S}} τ_{red}(e_{x_{i},1}) = ∑_{x_{i}∈\overline{S}} τ_{red}(e_{x_{i},2}),
⇔ τ(∑_{x_{i}∈S} \frac{e_{x_{i},1}}{t_{i}} − t_{∞}e_{∞,1}^{∗}) = τ(∑_{x_{i}∈S} \frac{e_{x_{i},2}}{t_{i}} − t_{∞}e_{∞,2}^{∗}),
⇔ ∑_{x_{i}∈S} \frac{e_{x_{i},2}−e_{x_{i},1}}{t_{i}} − t_{∞}(e_{∞,2}^{∗} − e_{∞,1}^{∗}) = \frac{g′}{g} , for g ∈ \overline{𝔽}_{p}(x)^{∗}.

This leads to: ∀_{x_{i}∈S}, e_{x_{i},2}[p] − e_{x_{i},1}[p] ∈ 𝔽_{p} and e_{∞,2}^{∗}[p] − e_{∞,1}^{∗}[p] = 0 which, by Fuchs’ relation, implies that also e_{∞,2}[p] − e_{∞,1}[p] ∈ 𝔽_{p}. On the other hand, the hypothesis that Cl(y_{1}) ≠ Cl(y_{2}) means that ∃ i ∈ \overline{S}, e_{x_{i},2} − e_{x_{i},1} ∉ ℚ which contradicts condition (∗).

5 Modular improvements of the combinatorial problem

Definition 5.1. Let L ∈ \overline{ℚ}(x)[∂]. L is said to be rad-regular if there exists r ∈ \overline{ℚ}(x) such that L_{∂→∂+r} has all its unramified generalized exponents in ℚ. Otherwise, L is said to be rad-singular; it means that there exists a singularity x_{i} with two unramified generalized exponents whose difference is not in ℚ.

As mentioned in the introduction, we want to address two problems in Beke’s algorithm: the combinatorial problem and the field problem. In this section we focus only on the combinatorial problem and not on the field problem; we suppose that a number field C is given and develop an algorithm that finds all exponential solutions defined over C. We show how a careful modular analysis can effectively reduce the combinatorial problem if L is rad-singular. Unfortunately, in the rad-regular case no reduction in the combinatorial problem is obtained (see also Section 7.2).

5.1 Some precisions on the vocabulary

Let C be a number field. In light of Remark 1.6 we use the terminology:

Definition 5.2. An exponential function y is defined over C if y′/y ∈ C(x).

When C is not algebraically closed, we will not consider all roots of the leading coefficient a_{n} of L, we will only consider one root for each irreducible factor over C. So the number of singularities that our algorithm considers is the number of irreducible factors of a_{n} (this number depends on the field C) plus one (the point at infinity). The only extensions of C that we will work with in Section 5 are extensions given by a single irreducible factor of a_{n} in C[x], so Algorithm ExpSolsInC in Section 5.4 does not use computationally expensive nested extensions.

Definition 5.3. A place is either a monic irreducible polynomial P_{i} in C[x] or ∞. At each place P_{i} ≠ ∞, we define the point x_{i} as the class of x in C[x]/(P_{i}) noted RootOf(P_{i}). At the place ∞, the point is also ∞.

If the degree of P_{i} is one, then we can view x_{i} as an element of C. Recall that we have the local parameter t_{i} = x − x_{i} at finite points and t_{i} = 1/x at ∞. We adapt the notion of singularity and generalized exponents to the purpose of this section:

Definition 5.4. The C-singularities of L ∈ C(x)[∂] written as in (3) are the points at the places P_{i} for all monic irreducible factors P_{i} of a_{n} in C[x] and the point ∞. We keep the notations S for the set of finite C-singularities and \overline{S} = S ∪ \{∞\}. Furthermore, a C-generalized exponent of L at x_{i} is an unramified generalized exponent of L at the C-singularity x_{i} that further belongs to C(x_{i})[t_{i}^{−1}]. The degree of a C-singularity x_{i} is the degree of P_{i} if x_{i} ≠ ∞, and it is 1 if x_{i} = ∞.

A first order operator ∂ − r ∈ C(x)[∂] has a C-generalized exponent at every C-singularity. Hence, if at a C-singularity of L there are no C-generalized exponents, then L can not have exponential solutions defined over C.

Definition 5.5. A C-combination is a list composed of one C-generalized exponent e_{x_{i}} at each C-singularity x_{i} ∈ \overline{S}.

The exponential solutions that we are looking for (those defined over C) match first order right hand factors ∂ − r with r ∈ C(x) as in Lemma 1.8 but here we write:

r = ∑_{x_{i}∈S} Tr_{i}(\frac{e_{x_{i}}}{t_{i}}) − t_{∞}e_{∞}^{∗} + \frac{Q′}{Q}  (7)

where Tr_{i} is the trace over the field extension C(x_{i}, x) ⊃ C(x). The term Tr_{i}(e_{x_{i}}/t_{i}) ∈ C(x) is the sum of all conjugates of e_{x_{i}}/t_{i} over C(x). So the above r is as in Lemma 1.8 with the additional property r ∈ C(x), which means that for any σ in the Galois group of \overline{ℚ}(x) over C(x), if we choose e_{x_{i}} at x_{i} in Lemma 1.8 then we must choose σ(e_{x_{i}}) at σ(x_{i}) (σ maps the unramified generalized exponents of L ∈ C(x)[∂] at x_{i} to those at σ(x_{i})). We extend Definition 3.1 to finite C-singularities x_{i} as follows:

τ_{red}(e_{x_{i}}) := τ(Tr_{i}(\frac{e_{x_{i}}}{t_{i}})).

Definition 5.6. Let P be a problem defined over a field C of characteristic 0 for which the following holds: if s is a solution of P defined over an algebraic extension of C, then all conjugates of s over C are solutions as well. Then we say that s_{1}, …, s_{k} are the solutions of P up to conjugation over C when:

- for each s_{i} a finite extension C_{i} of C is given such that s_{i} is defined over C_{i} and,

- for any solution s of P defined over any algebraic extension C′ of C there is precisely one i such that C_{i} can be embedded in C′ over C and s_{i} corresponds to s under this embedding.

We define the degree (or algebraic degree) of s_{i} over C as [C_{i} : C]. The total number of solutions defined over \overline{C} equals the sum of the degrees over C of the s_{i}.

Each such problem P that we need to solve (such as: find the singularities, find the generalized exponents at a singularity x_{i}, find the exponential solutions) will be solved up to conjugation, this in order to prevent computationally expensive splitting fields. In this terminology, “computing the C-singularities” and “computing the singularities up to conjugation over C” have identical meaning: factor a_{n} over C and construct a field C_{i} = C(x_{i}) = C[x]/(P_{i}) for each factor P_{i} (the singularity at infinity has degree 1 hence its field is C). Solving a problem up to conjugation over C does not involve constructing fields that contain more than one C_{i}, which means that we are not able to perform arithmetic between objects defined over distinct C_{i}’s (unless of course when the degrees are 1, i.e., when these C_{i}’s are just C). That is why the e_{x_{i}}/t_{i} in equation (7) must first be “traced down” to an object defined over C before we can add them.

When we compute generalized exponents at x_{i} (recall that we only need the unramified generalized exponents), we compute them up to conjugation over C(x_{i}). The C-generalized exponents are those of (algebraic) degree 1 over C(x_{i}). Generalized exponents of degree > 1 will not be used in Section 5, they will be used in Algorithm FindASol in Section 6.2.

5.2 Finding C-combinations matching the modular information

Let L ∈ C(x)[∂] where C is a number field. Suppose that L has m finite C-singularities x_{1}, …, x_{m} and consider the following set:

Comb := \{(\overline{e_{x_{1}}}, …, \overline{e_{x_{m}}}, \overline{e_{∞}}) | e_{x_{i}} is a C-generalized exponent of L at x_{i} \}.

The elements of Comb are called C-combinations mod ℚ. Note that the set Comb depends on C in two ways. First, the x_{i} correspond to irreducible factors of a_{n} over C. Second, the e_{x_{i}} must be defined over C(x_{i}). Consider now the map

\begin{aligned}Φ : & Comb & → & \overline{𝔽}_{p}(x^{p}), \\ & (\overline{e_{x_{1}}}, …, \overline{e_{x_{m}}}, \overline{e_{∞}}) & ↦ & ∑_{x_{i}∈\overline{S}} τ_{red}(\overline{e_{x_{i}}}). \end{aligned}

Definition 5.7. Let s ∈ R(χ_{p}(L)) and h ∈ Comb. We say that h matches s if Φ(h) = s. We say that h satisfies Fuchs’ relation if ∑_{x_{i}∈S} Tr_{i}(Const(\overline{e_{x_{i}}})) + Const(\overline{e_{∞}}) = 0 mod ℚ.

Given an element s of R(χ_{p}(L)) for a good prime p, we can effectively find the finite set of elements h of Comb that match s. A naive way to do that consists in taking the image under Φ of all possible C-combinations mod ℚ and keeping those with image s. A reason for including Condition (C3) in the definition of a good prime is so that we can give a more practical method (Algorithm CombMatchRoot below).

In the sequel, 𝔽 := C[p] denotes the image of the number field C after reduction mod p. In practice C is given by some (possibly nested) RootOf’s so that this reduction makes sense (recall that choosing a reduction [p] means choosing factors mod p of the polynomials defining the RootOf’s).

Since 𝔽 is finite, the map U ↦ U^{p} is an isomorphism from 𝔽[x] to 𝔽[x^{p}], and we denote the inverse isomorphism by V ↦ V^{1/p}. If r = \frac{a}{b} +c with a, b, c ∈ 𝔽[x] and deg(a) < deg(b) then τ(r) can be written as \frac{A}{B} + C with deg(A) < deg(B), B = b^{p} (here gcd(A, B) need not be 1) and C = τ(c). So if V ∈ 𝔽[x^{p}] divides the denominator of τ(r) then V^{1/p} ∈ 𝔽[x] divides the denominator of r.

Algorithm CombMatchRoot

Input: a number field C, a good prime p for L, the reduction [p], s ∈ R(χ_{p}(L)) ∩ 𝔽(x^{p})

and the data structure of C-singularities and C-generalized exponents of L. Output: the elements h of Comb matching s and satisfying Fuchs’ relation. 1 - For each x_{i} ∈ \overline{S} and each e_{x_{i}} at x_{i} compute τ_{red}(e_{x_{i}}). 2 - Compute the partial fraction decomposition ∑_{i=1}^{N_{s}} \frac{U_{i}}{V_{i}^{n_{i}}} + W of s with U_{i}, W ∈ 𝔽[x^{p}],

V_{i} irreducible in 𝔽[x^{p}] and deg(U_{i}) < deg(V_{i}^{n_{i}}). 3 - Remove all \overline{e_{∞}} for which τ_{red}(\overline{e_{∞}}) ≠ W; if none remain, then return ∅ and stop. 4 - For each i = 1 … N_{s}, 4a - Find the C-singularity x_{j} = RootOf(P_{j}) such that V_{i}^{1/p} ∈ 𝔽[x] divides P_{j}[p], 4b - Remove all \overline{e_{x_{j}}} for which the partial fraction decomposition of τ_{red}(\overline{e_{x_{j}}}) does not

contain the term \frac{U_{i}}{V_{i}^{n_{i}}}; if none remain, then return ∅ and stop. 5 - For all remaining \overline{e_{x_{j}}}, compute all combinations h = (\overline{e_{x_{1}}}, …, \overline{e_{x_{m}}}, \overline{e_{∞}}) that satisfy

Fuchs’ relation and return those as output.

In step 1, if the reduction mod p fails then e_{x_{i}} can be discarded by the same argument as in Lemma 3.4. In steps 3 and 4b, if nothing remains, then there are no first order right hand factors ∂ − r of L with r ∈ C(x) and s = τ(r). Condition (C3) in the definition of good primes is used in step 4a to make sure that the P_{j} for which V_{i}^{1/p} divides P_{j}[p] is unique (the existence of such a P_{j} follows from the fact that a pole of a root of χ_{p}(L) is a singularity of L[p]).

Remark 5.8. It follows from Lemma 3.4 that if e_{x_{i}} can not be reduced mod p then it can not be relevant for the algorithm, it can not correspond to any first order right hand factor. It follows from Algorithm CombMatchRoot above that a generalized exponent that can be reduced mod p but that does not match any element of R(χ_{p}(L)) can not be relevant either. If, adding these notions of relevant to those of Remark 1.9 ( i.e., unramified, ℤ-minimal), no relevant generalized exponents remain at a certain singularity x_{i} then there are no exponential solutions and we can stop the computation. If at a certain x_{i} only one relevant generalized exponent remains then we call x_{i} a semi-apparent singularity because it does not contribute to the combinatorial problem. In Algorithm FindASol in Section 6 it is important that we use only relevant generalized exponents to prevent making unnecessary field extensions.

Algorithm CombMatchModp:

Input: an operator L ∈ C(x)[∂], a number field C, the data structure E of C-singularities and C-generalized exponents of L. Output: the set of C-combinations mod ℚ that match a root of χ_{p}(L) for a good prime p and that satisfy Fuchs’ relation. 1 - If at some C-singularity there are no C-generalized exponents, return ∅ and stop. 2 - Choose a good prime p for L (choose a reduction [p]). 3 - Compute R(χ_{p}(L)) and (to be used in Section 7) Nroots := the total number of roots

of χ_{p}(L) in \overline{𝔽}_{p}(x^{p}) counted with multiplicity. 4 - Return the union of CombMatchRoot(C,p,[p],s,E) for s ∈ R(χ_{p}(L)) ∩ 𝔽(x^{p}).

One can often (use Condition (∗) in Proposition 4.5) find a prime such that each s matches at most one C-combination mod ℚ (i.e., the set CombMatchRoot(C,p,[p],s,E) has ≤ 1 elements for each s). But Proposition 4.4 shows that such prime does not always exist so that, in practice, we do not try to find one.

5.3 Some remarks on the choice of p

We can choose p as follows. Start with p = 2, and as long as a reduction mod p of a coefficient of L fails, or if one of the other conditions of a good prime in Definition 3.3 is not met, we take the next prime. So we compute the smallest good prime. There are, however, certain cases where another choice could be better. In hard cases (when we need to use Algorithm FindASol, or when the number of combinations is very high) we may want to try a few more primes and select the best one. In very easy examples such as in Section 2.1 it can be best to use no primes at all because the combinatorial problem was already practically trivial without modular methods. Another situation where there is no improvement in the combinatorial problem is when L is rad-regular (see Definition 5.1 and also Section 7.2) and the p-curvature computation does not exclude existence of exponential solutions. In most cases we only use the smallest good p, and in the unusual event that this p is large then we skip the modular improvements given in Section 5 because p-curvature computations become expensive for large p. Exceptions are very easy cases (use no modular methods) and very hard cases (check more than one prime and select the best one).

5.4 Finding all exponential solutions defined over a given field

Algorithm CombMatchModp computes a set of possible combinations for C-generalized exponents mod ℚ that could appear in an exponential solution defined over C. We will use this to find all exponential solutions defined over C.

Recall that with the notations of Equation (7) and Lemma 1.8, an exponential solution can be written as exp(∫ S) Q where S = ∑_{x_{i}∈S} Tr_{i}(\frac{e_{x_{i}}}{t_{i}}) − t_{∞}e_{∞}^{∗} and Q is a polynomial. After CombMatchModp, it remains to check which C-combinations mod ℚ lead to exponential solutions and to compute these exponential solutions. A natural idea is: at each C-singularity, given a \overline{e_{x_{i}}} we take the finitely many generalized exponents e_{x_{j}} satisfying e_{x_{j}} mod ℚ = \overline{e_{x_{i}}} and we have thus a finite number of C-combinations (e_{x_{1}}, …, e_{x_{m}}, e_{∞}), hence a finite number of possibilities for S after checking Fuchs’ relation. Then for each S, we construct the operator L_{∂→∂+S} having as solutions those of L divided by exp(∫ S), and we search for polynomial solutions Q of that operator.

We define an E-combination as a C-combination that uses only generalized exponents that appear in our data structure E. As already mentioned in Remark 1.9, replacing e_{x_{i}} by e_{x_{i}} − d and Q by t_{i}^{d}Q for some d ∈ ℤ leaves the expression exp(∫ S) Q invariant. So we only need one generalized exponent in each equivalence class mod ℤ. To assure that Q is a polynomial we only put ℤ-minimal generalized exponents in our data structure E, i.e., unramified generalized exponents that are minimal in their equivalence class mod ℤ. This is also true for the point at infinity; although the constant term of e_{∞} does not contribute to S in Lemma 1.8, it does contribute in Equation (4) to the degree bound for Q (the number N in step 4 is used as a degree bound for the Q_{i} in step 4b), and using only ℤ-minimal generalized exponents at infinity ensures that this degree bound is not too low.

Algorithm ExpSolsInC:

Input: an operator L ∈ C(x)[∂], a number field C.

Output: Sol, a basis of exponential solutions of L defined over C.

1 - Compute the data structure E of C-singularities and ℤ-minimal C-generalized exponents. 2 - M := CombMatchModp(L,C,E). 3 - Sol := ∅. 4 - For each E-combination (e_{x_{1}}, …, e_{x_{m}}, e_{∞}) whose class mod ℚ is in M, and that satisfies:

N := −Const(e_{∞}) − ∑_{x_{i}∈S} Tr_{i}(Const(e_{x_{i}})) ∈ ℕ, do: 4a - Let S := ∑_{x_{i}∈S} Tr_{i}(\frac{e_{x_{i}}}{t_{i}}) − t_{∞}e_{∞}^{∗} and L̃ := L_{∂→∂+S}, 4b - Compute a basis Q_{1}, …, Q_{w} of polynomial solutions in C[x] for L̃, 4c - If w > 0 then add the exp(∫ S) Q_{j} to Sol. 5 - Return Sol.

We can pre-compute the Tr_{i}(Const(e_{x_{i}})) and Tr_{i}(\frac{e_{x_{i}}}{t_{i}}), from then on all the work in the loop in step 4 (which may dominate the computation time, see also Section 7.1) takes place over C.

To use this procedure in the recursive algorithm FindASol in Section 6.2, we need to allow two extra inputs. First, an option called “just one”; if this option is given then we only need to return just one (if it exists) exponential solution defined over C but not a totally arbitrary one. We return a solution exp(∫ S) Q with S as in the algorithm, but with Q of minimal possible degree. Second, an optional set F may be given in the input. Each element of F is a couple (x_{i}, e_{x_{i}}) where x_{i} is a C-singularity of degree 1 (see Definition 5.4) and e_{x_{i}} is a C-generalized exponent at x_{i}. This couple encodes the command: “in the data structure E in step 1 of ExpSolsInC, throw away all C-generalized exponents at x_{i} except e_{x_{i}}”. This reduces the number of E-combinations to be checked. Note that usually not all C-singularities will appear in F, so in general there is still a combinatorial problem. Furthermore, if F ≠ ∅, then it is of course possible that no exponential solutions over C are found even if such solutions exist.

6 A way to handle the field problem

Let L ∈ C(x)[∂] with C ⊂ \overline{ℚ}. We want to find a basis of all exponential solutions, so we may have to search for exponential solutions defined over algebraic extensions of C. Before constructing such extensions, we first want to bound their degrees.

6.1 Two ways to obtain bounds

We start with a few definitions, in particular we need the notion of type. Many useful properties of the type of a differential operator can be found in [Ore32] (where type is called Art-Begriff). We will only use the type for first order operators, or equivalently, for exponential solutions:

Definition 6.1. Two exponential functions y_{1} and y_{2} are said to be of the same type if y_{1}/y_{2} is a rational function. This defines an equivalence relation, and type(y_{1}), the type of y_{1}, is defined to be the equivalence class of y_{1}. We say that y_{1} and y_{2} have the same local type at x_{i} if y_{1}/y_{2} is meromorphic at x_{i}, which is equivalent to saying that the generalized exponents of ∂ − y′_{1}/y_{1} and ∂ − y′_{2}/y_{2} at x_{i} differ by an integer.

We will identify an exponential solution y with its minimal operator ∂ − y′/y ∈ \overline{ℚ}(x)[∂], and so we shall not distinguish y from a constant multiple of y. This makes it easy to describe what we mean by a conjugate of y; if σ is an element of the Galois group of \overline{ℚ} over some number field C, then by σ(y) we simply mean a non-zero solution of ∂ − σ(y′/y).

Two exponential functions have the same type if and only if they have the same local type at every point x_{i}. The type of the exponential solution y = exp(∫ S) Q with S and Q as in Lemma 1.8 can be uniquely represented as ∑_{x_{i}∈S} \frac{e_{x_{i}} mod ℤ}{t_{i}} − t_{∞}e_{∞}^{∗} where e_{∞}^{∗} is e_{∞} without its constant term (recall that this constant term depends on the constant terms of the other e_{x_{i}} by Fuchs’ relation). The Galois group of \overline{ℚ} over C acts on the set of all types.

Definition 6.2. Let C be a number field, y be an exponential function. Let K be the field of definition of y over C, which is defined as the smallest number field K that contains C for which y′/y ∈ K(x). Then the algebraic degree of y over C is [K : C]. We say that y is of minimal algebraic degree m over C if [K : C] = m and there exist no exponential function ỹ, of the same type as y, having smaller algebraic degree over C.

Let C be a number field, let y be an exponential function, and let K be the field of definition of y over C. Let L_{1} = ∂ − y′/y ∈ K(x)[∂]. Let L_{1}, …, L_{r} ∈ \overline{ℚ}(x)[∂] be the conjugates of L_{1} over C, and let L = LCLM(L_{1}, …, L_{r}) ∈ C(x)[∂]. It is easy to calculate a basis of all exponential solutions of L that have the same type as y; compute L̃ := L_{∂→∂+y′/y} ∈ K(x)[∂], compute a basis b_{1}, …, b_{k} ∈ K(x) of rational solutions of L̃, then take b_{1}y, …, b_{k}y. Let R ∈ K(x)[∂] be the monic operator with b_{1}y, …, b_{k}y as solutions. Now R is the unique monic operator whose solution space is the set of all solutions of L with the same type as y, so R is uniquely determined by two things: L ∈ C(x) and type(y). Hence, any element σ of the Galois group of \overline{ℚ} over C that leaves type(y) invariant leaves R invariant as well. Conversely, if σ leaves R invariant, it leaves type(y) invariant as well because R has precisely one type of exponential solutions. Thus, one can define a field of definition of type(y) over C, namely as the smallest field C′ that contains C for which R ∈ C′(x)[∂].

Note that since R ∈ K(x)[∂], one has C′ ⊆ K, hence: the field of definition of y over C contains the field of definition of type(y) over C.

We will now prove constructively that y is of minimal algebraic degree over C if and only if the fields of definition of y and type(y) are the same. To do this, we will show that R has an exponential solution defined over C′, in fact, we will show that the algorithm ExpSolsInC with input R, C′ will find such solutions. Since R has only one type, there is exactly one local type at each singularity x_{i}, so generalized exponents at x_{i} are unique modulo integers. So if σ is in the Galois group of \overline{ℚ} over C′(x_{i}), then it can only move a generalized exponent e_{x_{i}} by an integer, but this implies that σ leaves e_{x_{i}} invariant, and thus e_{x_{i}} ∈ C′(x_{i})[t_{i}^{−1}], i.e., e_{x_{i}} is a C′-generalized exponent. Hence there is exactly one ℤ-minimal C′-generalized exponent at every singularity x_{i}, which implies that ExpSolsInC(R,C′) will try exactly one C′-combination. This combination has the same type as y, and since all solutions of R have this type, ExpSolsInC will find a basis y_{1}, …, y_{k} of solutions of R, defined over C′, where k is the order of R. We conclude the following:

Lemma 6.3. Let C be a number field and L ∈ C(x)[∂]. If y is an exponential solution of L defined over some algebraic extension K of C, then its conjugates over C are also solutions of L. Given L and y, we can compute y_{1}, …, y_{k} of minimal algebraic degree over C that form a basis of all exponential solutions of L of the same type as y. If y is of minimal algebraic degree over C, then the distinct conjugates of y over C are of distinct type and hence linearly independent.

Proof. The fact that for L ∈ C(x)[∂], the conjugates over C of solutions are again solutions is clear. From the foregoing, given L and y, we can compute an operator R ∈ K(x)[∂] whose solution space is precisely the set of all exponential solutions of L of the same type as y. We can calculate the smallest field C′ containing C for which R ∈ C′(x)[∂], and can find a basis y_{1}, …, y_{k} of solutions of R with ExpSolsInC(R,C′). Then y_{1}, …, y_{k} are defined over C′, which is the field of definition of their type, so they have minimal algebraic degree over C. And they form a basis of all exponential solutions of L of the same type as y.

If σ is in the Galois group of \overline{ℚ} over C, and if two objects have the same field of definition over C, then σ leaves the first object invariant if and only if σ acts as the identity on the field of definition, if and only if σ leaves the second object invariant. Now if y is of minimal algebraic degree over C, then it has the same field of definition as type(y), so σ leaves type(y) invariant if and only if σ leaves y (which we identify with ∂ − y′/y because we do not want to distinguish y from scalar multiples of y) invariant. Thus, distinct conjugates of y must have distinct types.

Therefore, if we have an exponential solution of minimal algebraic degree m, then this exponential solution gives in fact m linearly independent solutions. In the output of our algorithm ExpSols we will only return exponential solutions up to conjugation over C. In order to make it easier to count how many linearly independent solutions such an output represents, we will only allow exponential solutions in the output of ExpSols that have minimal algebraic degree over C.

Definition 6.4. If C is a number field, x_{i} a singularity of L and e_{x_{i}} a generalized exponent at x_{i}, then the field C(x_{i}; e_{x_{i}}) is the extension of C given by x_{i} and the coefficients of e_{x_{i}}.

Proposition 6.5. Let C be a number field, L ∈ C(x)[∂], x_{i} a singularity, t_{i} the local parameter, and let e ∈ \overline{ℚ}[t_{i}^{−1}]. Let b_{1} be the number of distinct generalized exponents e_{x_{i}} of L at x_{i} with e_{x_{i}} − e ∈ ℤ. Suppose that y is an exponential solution of minimal algebraic degree m over C(x_{i}; e) with generalized exponent e at x_{i}. Then m ≤ b_{1}.

Proof. By the previous Lemma, the conjugates of y over C(x_{i}; e) form m linearly independent exponential solutions, all of which have generalized exponent in e + ℤ at x_{i}. These exponential solutions (viewed as formal solutions at x_{i}) can be written as ϕ_{j}(t_{i}) exp(∫ e/t_{i} dt_{i}) where ϕ_{j}(t_{i}) ∈ \overline{ℚ}((t_{i})) for all j ∈ \{1, …, m\}. After Gaussian elimination on the vectors of coefficients of the ϕ_{j}(t_{i}), we obtain m formal solutions at x_{i} with distinct generalized exponents in e + ℤ. Then m ≤ b_{1} because there are only b_{1} such generalized exponents at x_{i}.

The modular approach yields a different bound:

Proposition 6.6. Let C be a number field, let L ∈ C(x)[∂] and let p be a good prime for L. Let b_{2} be the number of roots of χ_{p}(L) (counted with multiplicity) that match (in the sense of Algorithm CombMatchRoot) the choice \overline{e_{x_{i}}} at x_{i}. Let y be an exponential solution of minimal algebraic degree m over C(x_{i}; e_{x_{i}}), and generalized exponent e_{x_{i}} at x_{i}. Then m ≤ b_{2}.

Proof. Let R ∈ C(x_{i}; e_{x_{i}})(x)[∂] be the right hand factor whose solutions are spanned by all exponential solutions y of L that have generalized exponent e_{x_{i}} at x_{i}. Since there are at least m independent such solutions, the order d of R is at least m. All solutions of R have the same local type at x_{i} (namely the local type e_{x_{i}} mod ℤ, see Definition 6.1). Now R factors as a product of first order factors R = R_{1}R_{2} ⋯ R_{d} because R has a basis of exponential solutions. All R_{i} have the same local type at x_{i}, so for each R_{i}, the generalized exponent at x_{i} must be congruent to e_{x_{i}} mod ℤ. Hence the root of χ_{p}(R_{i}) matches \overline{e_{x_{i}}} at x_{i}. The proposition now follows from the fact that R and R_{1}, …, R_{d} can be reduced mod p (see Proposition 1.15) and Lemma 1.13(a).

We now have two bounds b_{1} and b_{2} for the degree of the extension needed over C(x_{i}; e_{x_{i}}). We do not know in advance which is smaller, so we compute both and take the minimum.

Definition 6.7. To each couple (x_{i}, e_{x_{i}}) we associate the bound b_{(x_{i},e_{x_{i}})} := min(b_{1}, b_{2}).

6.2 An algorithm to find an exponential solution over an algebraic extension

Let L ∈ C(x)[∂] with C ⊂ \overline{ℚ}. We suppose that L has no exponential solutions defined over C and we want to find (if it exists) an exponential solution defined over an algebraic extension of C. We will use Algorithm ExpSolsInC as well as the bounds from the previous section. To find an exponential solution over an algebraic extension of C, we increase the field step by step by adding extensions coming from singularities or generalized exponents, and at each step we use ExpSolsInC to search for an exponential solution over the new field. The bounds tell us when we can stop making field extensions.

For efficiency reasons we first want to reduce the order of the problem as much as possible before introducing field extensions, which is precisely what Algorithm ExpSols in Section 7 does before it calls Algorithm FindASol below. And since the order can be reduced whenever a solution is found, we want FindASol to stop computing as soon as it finds a solution.

In the algorithm we denote e_{i,j}, j = 1, …, n_{i} as all, up to conjugation over C(x_{i}), see Definition 5.6, relevant (see Remarks 1.9 and 5.8) generalized exponents at x_{i}. If n_{i} = 0 then we can stop the algorithm, the output is ∅. Let F be a data structure as in the comments after Algorithm ExpSolsInC. If x_{i} appears in F, which means that a generalized exponent at x_{i} has already been chosen, then n_{i} := 1. We denote d_{i,j} := [C(x_{i}; e_{i,j}) : C(x_{i})], which is the degree as in Definition 5.6 of e_{i,j} over C(x_{i}). Now we define b_{(x_{i},e_{i,j},F)} as in Definition 6.7 with as difference that in Proposition 6.6 we only count those roots of χ_{p}(L) that match all generalized exponents specified in F. Of course, before we enter step 5 we discard all e_{i,j} for which b_{(x_{i},e_{i,j},F)} = 0 because we want n_{i} to be as small as possible, and if any n_{i} vanishes then we stop.

Algorithm FindASol:

Input: an operator L ∈ C(x)[∂], a field C, a data structure F and a positive integer b.

Output: an exponential solution of L or ∅.

1 - If b = 0, then return ∅ and stop.

2 - Run ExpSolsInC(L,C,options={“just one”, F }). 3 - If this gives a solution, then return it and stop. 4 - If b = 1, then return ∅ and stop. 5 - If we can choose x_{i} in step 6 with degree 1 over C, or if b′ ≤ b/2 for all b′ in step 6a,

then go to step 6, else go to step 8. 6 - Choose one C-singularity x_{i}, not appearing in F. Then, for j = 1, …, n_{i} do 6a and 6b: 6a - Let b′ be the minimum of b_{(x_{i},e_{i,j},F)} and ⌊\frac{b}{d_{i,j}} ⌋. 6b - If FindASol(L, C(x_{i}; e_{i,j}), F ∪ {(x_{i}, e_{i,j})}, b′) ≠ ∅, then return it and stop. 7 - Return ∅ and stop. 8 - Do step 8a for all x_{i} not appearing in F for which d_{i} := [C(x_{i}) : C] is not 1 and

there are at least two j’s with d_{i,j} = 1. 8a - If FindASol(L, C(x_{i}), F, ⌊b\frac{⌊d_{i}/2⌋}{d_{i}}⌋) ≠ ∅, then return it and stop. 9 - For all i, j with x_{i} not appearing in F and 1 < d_{i,j} ≤ b do step 9a: 9a - If FindASol(L, C(x_{i}; e_{i,j}), F ∪ \{(x_{i}, e_{i,j})\}, ⌊\frac{b}{d_{i,j}} ⌋) ≠ ∅, then return it and stop. 10 - Return ∅.

For the modular computations in recursive calls of ExpSolsInC, we always use the same prime number so that the modular information that we need is only computed once. One can also reuse information on singularities and generalized exponents from previous computations, but this information needs to be updated to the new field which involves factoring polynomials. For the singularities this means: factoring a_{n} over the new field. The generalized exponents need to be updated as well, for example, if we extend C then the degree over C(x_{i}) of a generalized exponent can decrease, which means that the polynomials defining the field extension C(x_{i}; e_{i,j}) over C(x_{i}) have become reducible, so they need to be factored. This way one conjugacy class of generalized exponents can split up into several conjugacy classes of generalized exponents.

We are going to verify (we must prove this recursively, i.e., by induction) that FindASol satisfies the following:

Algorithm specification: if an exponential solution y exists that satisfies the fixed choices in F and has degree ≤ b over C, then FindASol will return some exponential solution of L (but not necessarily of degree ≤ b over C).

Proof. We first show that step 6 is correct. Suppose such y exists. Then the conjugates of y over C are also solutions of L, and satisfy the same fixed choices made in F (because in F we only have singularities and generalized exponents defined over C) as well as the degree bound b. Now the e_{i,j} constitute, up to conjugation, a complete list of relevant ℤ-minimal generalized exponents at x_{i}, so there must exist j such that a conjugate of y has generalized exponent in e_{i,j} + ℤ at x_{i}. We may assume that conjugate is just y. So for this j, there is a solution y defined over a field C′ with [C′ : C] ≤ b having generalized exponent in e_{i,j} + ℤ at x_{i}. We now replace C by C(x_{i}; e_{i,j}) and F by F ∪ {(x_{i}, e_{i,j})} so we may use the degree bounds from the previous section. We included a new bound in the algorithm as well, namely ⌊\frac{b}{d_{i,j}} ⌋, which holds because:

e_{i,j} is a C′-generalized exponent since it corresponds (mod ℤ) to the solution y defined over C′. So C′(x_{i}) = C′(x_{i}; e_{i,j}). Let b̃ = [C′(x_{i}) : C(x_{i})]. Now C′(x_{i}) = C′(x_{i}; e_{i,j}) contains C(x_{i}; e_{i,j}), which has degree d_{i,j} over C(x_{i}), and so [C′(x_{i}) : C(x_{i}; e_{i,j})] = b̃/d_{i,j}. Now y is defined over C′, hence over C′(x_{i}), which is an extension of C(x_{i}; e_{i,j}) of degree b̃/d_{i,j} and this is ≤ ⌊\frac{b}{d_{i,j}} ⌋ because b̃ ≤ [C′ : C] ≤ b.

Now steps 8 and 9. Again suppose that an exponential solution y exists that satisfies the fixed choices in F and has degree ≤ b over C. Let again C′ be the field of definition of y over C. If every C′-combination was also a C-combination then ExpSolsInC(L, C) and ExpSolsInC(L, C′) would have to find an equal number of solutions; so FindASol will return a solution in step 3. Thus we may assume that not every C′-combination is a C-combination, so there are more C’-combinations than C-combinations. This can only happen when at least one of the following is true: (I). There are more C’-singularities than C-singularities, or: (II). There exists a C-singularity x_{i} at which there are more C′-generalized exponents than Cgeneralized exponents. Case (I) will be handled by step 8 and case (II) by step 9. Suppose that (I) holds. Then the minimal polynomial of some C-singularity x_{i} must be reducible over C′. This minimal polynomial has degree d_{i}. One of the factors over C′ must have degree d′ ≤ ⌊d_{i}/2⌋. By abuse of notation we denote x_{i} as a root of that factor as well. So then [C′(x_{i}) : C′] = d′, [C′ : C] ≤ b, so [C′(x_{i}) : C] ≤ d′b. Now C(x_{i}) is a field between C′(x_{i}) and C having degree d_{i} over C, so [C′(x_{i}) : C(x_{i})] = [C′(x_{i}) : C]/d_{i} ≤ d′b/d_{i} ≤ bϵ where ϵ := ⌊d_{i}/2⌋/d_{i}. Note that ϵ = 1/2 when d_{i} is even, and ϵ < 1/2 when d_{i} is odd. Now y is defined over C′, hence also over C′(x_{i}), which is an extension of degree ≤ bϵ over our new field C(x_{i}). Hence the degree bound used in step 8a is correct.

Note that if there exists an i for which there is no j with d_{i,j} = 1 then a solution y necessarily involves a generalized exponent at x_{i} that is not defined over C(x_{i}), the case that is treated by step 9, and in this case the algorithm stays correct (but becomes more efficient) if we skip 8 and use only this i in step 9. And if there is only one j with d_{i,j} = 1, then either we must choose this j (at x_{i} and its conjugates over C), or the problem is handled by step 9 (in which case it is enough to use only this i in 9).

Increasing the number of singularities does not increase the number of combinations if there is only one choice to be made at these singularities, thus skipping i in step 8 when there is only one j with d_{i,j} = 1 is correct.

Now assume that (I) does not hold (then every C′-singularity is a C-singularity) and assume that (II) holds. Then at some C-singularity x_{i} we have a C’-generalized exponent e_{x_{i}} that is not a C-generalized exponent. This generalized exponent is defined over C′(x_{i}) but has degree d_{i,j} over C(x_{i}). The degree bound ⌊\frac{b}{d_{i,j}} ⌋ in step 9a can then be proved with the same arguments as before, which completes the proof that FindASol satisfies its specification.

The degree bound drops at least a factor 2 each time we use recursion with an extension of C. So if the initial degree bound b is less than 4 then the extensions we construct are not greater than those given by a single x_{i} (in step 8) or a single couple (x_{i}, e_{x_{i}}) in step 9. If b < 8 then we will stack at most 2 such extensions (each given by an x_{i} or by a couple (x_{i}, e_{x_{i}})) on top of each other, which is generally feasible in practice. So unless the degree bound is very high we do not expect the algorithm to choke on these algebraic extensions. And even if the degree bound is very high (order ≥ Nroots ≥ 8), the alternatives currently used by Maple’s expsols need not be better.

We can further improve the algorithm with modular methods. For example, certain x_{i} in step 8 may be omitted on the basis of modular information, namely when none of the roots of χ_{p}(L) allows a combination that uses distinct generalized exponents at two singularities conjugated to x_{i}.

If a root of χ_{p}(L) can only match a solution defined over C (this can happen when it must use an e_{i,j} for which C(x_{i}; e_{i,j}) = C and the bound in Definition 6.7 is 1) then we can discard this root after step 2. Reducing the number of roots of χ_{p}(L) in this way reduces the bounds and the n_{i}’s which may help prevent making some unnecessary field extensions.

Step 8a can be improved in the following way. If the recursive call for x_{i} returns ∅, then later recursive calls in step 8a should be prevented from factoring the minimal polynomial of x_{i} over the new fields, and prevented from making that same extension x_{i} again.

A similar improvement is possible in step 6b. If the recursive call for e_{i,j} returns ∅, then we could discard e_{i,j}. This will prevent later recursive calls from trying e_{i,j} again. But it also reduces the number of generalized exponents, which may improve the bound b (see the number N_{i} in comment [1] after Algorithm ExpSols in Section 7).

7 An algorithm to find all exponential solutions

In this section, we give a complete algorithm to find a basis, up to conjugation, for all exponential solutions of a differential operator L ∈ C(x)[∂] where C ⊂ \overline{ℚ}. The general idea in this recursive algorithm is the following. Using our procedure ExpSolsInC, we first compute the “easy” first order factors, that is, those defined over C. Then by looking at the value of Nroots (from Algorithm CombMatchModp) we check if there may exist more exponential solutions, which are then necessarily not definable over C. If such exponential solutions could exist then we first recursively remove all “easy” factors, both on the right as well as on the left using the adjoint L^{∗} of L (see [Inc26, 5.3] or [PS03, 2.1]) until there are no more easy factors on either side. We do this because we want to make the order of the operator as small as possible before entering the “hard” case (i.e., searching for exponential solutions defined over algebraic extensions using Algorithm FindASol).

The specifications of the following ExpSols algorithm are the following: first, it should not return the “same” solution more than once, more precisely: the output gives only one exponential solution in each conjugacy class. This is because we want the output to be a basis up to conjugation. Second, the solutions given in the output must all be of minimal algebraic degree over C, because this guarantees (see Lemma 6.3) that their conjugates over C are linearly independent. To count the number of elements in the basis of exponential solutions, we must count a solution y given in the output d(y) times, where d(y) is the degree over C of the field of definition of y. In the following, Card(Sol) denotes the number of elements of the set Sol. The number of solutions represented by Sol is the sum of d(y) taken over y ∈ Sol. Finally, we mark with a [.] the steps for which explanations or comments are given after the algorithm.

Algorithm ExpSols: Input: a linear differential operator L ∈ C(x)[∂] (C ⊂ \overline{ℚ}) and the field C. Output: Sol, a basis of exponential solutions of L up to conjugation over C.

– Sol := ExpSolsInC(L,C).

– If Card(Sol) ≥ Nroots −1, then return Sol and stop ([1]).

– If Sol ≠ ∅ then,

- Write L = L̃ LCLM(∂ − y′/y | y ∈ Sol),

- ExpSols(L̃,C) ([2]).

- Keep only the types not defined over C of the solutions found,

- For each remaining type represented by t ([3]), . Compute a basis R_{1}, …, R_{s} of rational solutions of L_{∂→∂+t}, . Add the exp(∫ t) R_{i} to Sol.

- Return Sol and stop.

– Else,

- If the order is 2 then go to the next Else below,

- Make L monic,

- ExpSolsInC(L^{∗},C) where L^{∗} is the adjoint of L ([4]).

- If it finds solutions y_{1}, …, y_{r} where r > 0 then, . Write L = LCLM(∂ − y′_{1}/y_{1}, …, ∂ − y′_{r}/y_{r})^{∗} L̃, . Return ExpSols(L̃,C) and stop.

- Else ([5]), . FindASol(L,C,∅,Nroots) ([6]). If ∅ then return ∅ and stop. . Optimize the solution y found ([7]), . Compute a basis R_{1}, …, R_{s} of rational solutions of L_{∂→∂+y′/y}, . Add the y R_{i} to Sol, . Write L = L̃ LCLM(∂ − r_{1}, …, ∂ − r_{s}, “and conjugates over C”), ([8]) . Remove recursively from L̃ the solutions of the same type as y ([9]), . ExpSols(L̃,C). . Keep only the types not defined over C of the solutions found, . For each remaining type represented by t ([3]),

Compute a basis R_{1}, …, R_{s} of rational solutions of L_{∂→∂+t}, Add the exp(∫ t) R_{i} to Sol.

. Return Sol and stop.

We explain or comment the points marked with a [.]:

[1] - The number Nroots, which was calculated in CombMatchModp, is used as upper bound for the number of linearly independent exponential solutions. We have already found Card(Sol) independent solutions and so there are at most d := Nroots − CardSol linearly independent exponential solutions left. This d is an upper bound for the degree of the field extension of remaining exponential solutions of minimal algebraic degree. If d ≤ 1 then no extensions are necessary so we can stop the computation. At a singularity x_{i}, let N_{i} be the number of distinct unramified generalized exponents in \overline{ℚ}[t_{i}^{−1}]: this is another upper bound for the total number of linearly independent exponential solutions. If N_{i} < Nroots, then replace Nroots by N_{i}.

[2] - At this point the roots of the p-curvature have already been computed, and need not be computed again; we can simply take the roots (with multiplicity) for χ_{p}(L) and reduce the multiplicities according to the solutions found.

[3] - We have applied recursion on a left factor of L so the exponential solutions found are in general not exponential solutions of L. We will only use the types of the exponential solutions of this left factor, not the exponential solutions themselves. We skip types defined over C because all exponential solutions over C have already been found. From [Ore32] or [Hoe97b, Lemma 7.1], the types of exponential solutions of L that are not yet found must be among the types in the left factor L̃ of L.

[4] - We will apply FindASol only when all easy factors have been removed both on the left and on the right. Removal of easy factors, left or right, does not cause solutions to be lost because we only do this after all easy solutions have already been computed and stored in the set Sol, in the first step of the algorithm. To find the easy left factors, we apply ExpSolsInC to the adjoint operator L^{∗} to find the easy exponential solutions of L^{∗} which correspond to easy left factors of L (see [PS03, 2.1]). Then we apply recursion on the remaining right hand factor.

[5] - If the order of the operator is two, and if the bounds in FindASol do not immediately rule out a solution defined over an extension, then the formulas given in [UW96] would be a good alternative to FindASol. If the order is three, then an alternative would be the eigenring method (see [PS03]) but we do not expect that to be better than FindASol, especially if one implements some shortcuts for order 3 in FindASol (such as: we only need to consider extensions in FindASol that have a degree 3 subfield. And: if at x_{i} there are two non-conjugated generalized exponents that do not differ by an integer, then C(x_{i}) is the only extension that needs to be considered).

[6] - Here, we reach the point where there is nothing left to do than entering the “hard case”, trying to find a solution over an extension of C with FindASol. An a priori first bound is Nroots, the number of roots in \overline{𝔽}_{p}(x^{p}) counted with multiplicity of χ_{p}(L) for a good prime p. Note that a more detailed analysis could lead to a sharper bound. For example, if two roots of χ_{p}(L) can not correspond to conjugated exponential solutions, or if one can conclude using the bounds in Section 6.1 that some root of χ_{p}(L) can not correspond to an exponential solution of minimal algebraic degree Nroots, then we could give FindASol a better bound.

[7] - Let y be an exponential solution. We can write y′/y = P_{1}/P_{2} where P_{1} and P_{2} are polynomials with gcd 1 and P_{2} is monic. The field of definition of y is then the field generated over C by the coefficients of P_{1} and P_{2}. By “optimizing the solution” we mean two things: 1). Using this solution y to find a solution of minimal algebraic degree over C. And 2). Making sure that the field that the algorithm gives for y (this field contains the field of definition of y) is actually equal to the field of definition of y. Both 1) and 2) are important. We want y to be of minimal algebraic degree so that we know that its conjugates are linearly independent. But the way we count the number of conjugates of y, i.e., the way we determine the number d(y), is not by looking at y, but by looking at the field given for y. The field provided for y by Algorithm FindASol contains, but need not be not equal to, the field of definition of y. So to optimize y, we take the field given for y, and then determine the subfield generated over C by the coefficients of P_{1} and P_{2}. Then we find defining equations (i.e., new RootOf’s) for this subfield, and use them to rewrite the coefficients of P_{1} and P_{2}. This then takes care of 2). To do 1), we could use the approach in Lemma 6.3, however, this is not necessary because the special choice that Algorithm ExpSolsInC makes when the option “just one” is given causes 1) to be automatically satisfied. The polynomial Q is of minimal degree with this option given, which leads to uniqueness of y′/y which in turn causes y to already be of minimal algebraic degree over C.

[8] - Here r_{i} = R′_{i}/R_{i} + y′/y. One can compute this LCLM using the method of undetermined coefficients and solving linear equations over C(x), this can be done without actually constructing the conjugates of y′/y.

[9] - What is meant here is to also remove solutions whose type is conjugated over C to the type of y. We will remove such solutions recursively because if we do not then it would be non-trivial to ensure the specification of the algorithm that the “same” solution is not returned more than once. So before we call ExpSols on the remaining left factor L̃ of L, we first remove from L̃ all exponential solutions of the same type as (conjugates of) y, and we keep repeating this until L̃ no longer has solutions of this type. A recursive procedure to achieve that is the following: compute the rational solutions R_{i} of L̃_{∂→∂+y′/y}; each one gives a first order right hand factor ∂ − r_{i} of L where r_{i} = R′_{i}/R_{i} +y′/y. If there are none then return L̃ else write L̃ = L̃_{1} LCLM(∂ −r_{1}, …, ∂ −r_{s}, “and conjugates over C”) and apply recursion on L̃_{1}.

7.1 Some remarks on the algorithm

Remark on ExpSolsInC: Even though we use Algorithm CombMatchModp, there could still be combinations in step 4 that do not lead to an exponential solution (if for example the number of combinations is greater than Nroots then we already know in advance that this will be the case). If the degree bound (the number N in step 4) for polynomial solutions is high, or if C is a complicated field, or if there are many combinations to be checked, then to make the algorithm efficient we need a quick way to discard non-solution-combinations before doing all the work in steps 4a and 4b, because each of these two steps can dominate the computation time.

In step 4, we can of course pre-compute the traces, so that all computations in step 4 are rational (i.e., over C). Computing S in step 4a then only involves additions in C(x). However, that does not imply that computing S is cheap, because normalizing S (writing S as P_{1}/P_{2} where P_{1}, P_{2} are polynomials with no common factors) can be an expensive operation if C is a complicated field. Computation of L̃ = L_{∂→∂+S} can be even more expensive because we need to multiply, differentiate, and add in C(x). And since step 4a is applied to each combination, one can easily spend more time in step 4a than in steps 1 and 2 combined. And if N is large then 4b can take even more time.

Often the N’s are very small, and one can improve the running time of ExpSolsInC substantially by implementing some cases for step 4b separately, such as the case N = 0 (in that case we should use a fast zero test for L̃(1) without computing L̃ itself) and the case N small non-zero: here, one would be tempted to take an ansatz ∑_{i=0}^{N} c_{i}x^{i}, and to compute the coefficients of L̃(ansatz) which gives linear equations for the unknowns c_{i}. But for small N and large L̃ this results in much more equations than we need because there are only N + 1 unknowns. In this case, to make the ansatz approach efficient we should not fully evaluate L̃(ansatz) because then we computed much more equations than we need. Instead we should evaluate just enough coefficients of L̃(ansatz) (and do this without fully computing L̃ itself) so that we have enough equations to determine the c_{i}, and once the c_{i} ∈ C are found then we finish with a quick zero test for L̃(∑_{i=0}^{N} c_{i}x^{i}). For large N we propose to compute L̃ and to use [ABP95].

If we find w > 0 polynomial solutions in step 4b then the computation time in 4a and 4b was well spent even if those steps took a lot of time. But when w = 0, and if C is a complicated field or N is a large number, then we want to avoid 4a and 4b with high probability. This is done as follows: take a sufficiently large prime number p (in general this p is not the same prime as we used for the p-curvature because for the p-curvature we only use small primes) such that C[p] = 𝔽_{p} (such primes have density at least 1/[C : ℚ] so they are easy to find). Then compute L̃[p] ∈ 𝔽_{p}(x)[∂] and check if it allows a polynomial solution of degree ≤ N. If not, skip 4a and 4b.

More general fields: It is not difficult to generalize our algorithm to arbitrary fields of characteristic 0. We may assume that C is the field generated by the coefficients of L, so then C is a finitely generated extension of ℚ, say C = ℚ(t_{1}, …, t_{k})[α_{1}, …, α_{l}] where t_{1}, …, t_{k} are algebraically independent and the α_{i} are algebraic over ℚ(t_{1}, …, t_{k}). Then the reduction mod p works as follows: choose random values for t_{i} in ℚ. This gives a map f_{1} from a subring R_{1} of ℚ(t_{1}, …, t_{k}) to ℚ. If not all coefficients of the minimal polynomial of α_{1} are in R_{1}, then choose different random values for t_{i}. After that, the minimal polynomial of α_{1} is mapped to a polynomial over ℚ, we can choose an irreducible factor, a root α′_{1} of that factor, and map α_{1} to α′_{1}. Do the same for α_{2}, …, α_{l}. This way we have a map f_{2} from a subring R_{2} of C to the field ℚ(α′_{1}, …, α′_{l}). We can then choose a reduction mod p in the same way as before, so we have a map g from a subring of ℚ(α′_{1}, …, α′_{l}) to a finite field. Then let the map [p] be the composition g ∘ f_{2} of these two maps, which will be defined on some subring R_{3} of R_{2}. One can always find evaluations in ℚ for the t_{i} and a prime number p such that all elements of C that we want to reduce mod p are in R_{3}.

Logarithms: Generalized exponents correspond to formal solutions. Since we are only interested in exponential solutions, the formal solutions that contain a logarithm are not of interest to us. Because of this, we can disregard those generalized exponents that do not correspond to a formal solution without a logarithm. Note that in any given equivalence class mod ℤ, this observation could eliminate some but never all unramified generalized exponents, so it does not help much for the combinatorial problem. However, this observation may reduce the number of distinct unramified generalized exponents in some equivalence classes mod ℤ, which can reduce the bound in Proposition 6.5 as well as the number N_{i} in comment [1] after Algorithm ExpSols. This could speed up Algorithm FindASol. Note that if L has an exponential solution of minimal algebraic degree n − 1 where n is the order of L, then the adjoint of L has an exponential solution defined over C. Algorithm ExpSols removes this solution, hence, if ExpSols calls FindASol with degree bound (the number b in FindASol) N_{i} = n − 1 then we may use n − 2 instead.

Deciding which generalized exponents match formal solutions without logarithms can be done if we compute d + 1 terms of the formal solutions, where d is the exponent-difference. We propose to do this computation only when ExpSols has to call FindASol. To speed up this computation, one could compute the formal solutions modulo a suitably large prime number, if a logarithm occurs then a logarithm must occur in characteristic zero as well.

A related implementation: Although we do not have an implementation for the algorithm ExpSols in this paper, there is an implementation in Maple 9 for the difference case, i.e., a procedure that computes the hypergeometric solutions of a difference equation. This implementation follows some of the ideas in this paper and performs very well in practice (a paper concerning the difference case is planned). We expect our ExpSols to perform equally well.

7.2 Computing radical solutions

Computing radical solutions is almost the same problem as computing exponential solutions, with only a few differences:

1. One only uses generalized exponents that are in ℚ, see Equation (5) and Remark 1.10.

2. τ_{red} maps rational numbers to zero, which implies that the only root of χ_{p}(L) that matters is 0, but also implies that there is no reduction in the combinatorial problem. This means that for computing radical solutions, the only useful modular information is the multiplicity of the root 0 in χ_{p}(L), which we can use as a bound for the number of linearly independent radical solutions, hence as a bound for the degree of the field extension over which such solutions (of minimal algebraic degree) are defined.

To save computation time, we can compute this multiplicity without computing χ_{p}(L), in the following way: Take a good prime p. Then compute a basis z_{1}, …, z_{s} of rational solutions of L[p]. If s = 0 then return 0, otherwise write L[p] = L_{1}L_{2} where z_{1}, …, z_{s} is a basis of solutions of L_{2}, apply recursion to L_{1}, then add s to the result for L_{1}, and return the answer. We expect this to be faster than computing χ_{p}(L) because rational solutions of a differential operator in characteristic p can be computed very quickly, O(max(l, p)^{2}p) field operations where l is a bound on the degree of the coefficients of L, see [Clu03].

7.3 Almost all primes

There are many results known that hold for all but finitely many p. For example, it is easy to show that

Lemma 7.1. For all but finitely many p, the number of linearly independent radical solutions of L is bounded by the dimension of the space of rational solutions of L[p].

Thus, for almost all p, we could have used the number s in the previous section as a bound. The reader may wonder why we did not do so. The problem here is that of the many results that are known for all but finitely many p, the only ones we can use for the algorithm are those for which a small prime p can be exhibited with little computation. This problem is also the reason that the only information we use from the p-curvature is its characteristic polynomial.

References

[ABP95] S. A. Abramov, M. Bronstein, M. Petkovšek. On Polynomial Solutions of Linear Operator Equations. In Proceedings of ISSAC’95, ACM, New York, 1995.

[Bek94] E. Beke. Die Irreduzibilität der homogenen linearen differentialgleichungen. In Math. Ann., 45, 278-294, 1894.

[Bou00] D. Boucher. Sur les Équations Différentielles Linéaires Paramétrées, une Application aux Systèmes Hamiltoniens. Thèse, Université de Limoges, 2000.

[Cha01] A. Chambert-Loir. Théorèmes d’Algébricité en Géométrie Diophantienne. Séminaire Bourbaki, exposé No. 886, Mars 2001.

[CC85] D. V. Chudnovsky, G. V. Chudnovsky. Applications of Pade Approximations to the Grothendieck Conjecture on Linear Equations. In Lecture Notes in Mathematics, Number Theory, 1135: 52-100, 1985.

[Clu03] T. Cluzeau. Factorization of Differential Systems in Characteristic p. In Proceedings of ISSAC’03, August 3-6, Philadelphia, Pennsylvania: 58-65, 2003.

[Gie98] M. Giesbrecht. Factoring in Skew-Polynomial Rings over Finite Fields. In Journal of Symbolic Computation, 24(5): 463-486, 1998.

[GZ03] M. Giesbrecht, Y. Zhang. Factoring and Decomposing Ore Polynomials over 𝔽_{p}(t). In Proceedings of ISSAC’03, August 3-6, Philadelphia, Pennsylvania: 127-134, 2003.

[Hoe97a] M. van Hoeij. Formal Solutions and Factorization of Differential Operators with Power Series Coefficients. In Journal Of Symbolic Computation, 24: 1-30, 1997.

[Hoe97b] M. van Hoeij. Factorization of Differential Operators with Rational Functions Coefficients. In Journal Of Symbolic Computation, 24: 537-561, 1997.

[HW97] M. van Hoeij, J.-A. Weil. An Algorithm for Computing Invariants of Differential Galois Groups. In Journal of Pure and Applied Algebra, 117-118: 353-379, 1997.

[Hon81] T. Honda. Algebraic Differential Equations. In Symp. Math., 24: 169-204, 1981.

[Inc26] E. L. Ince. Ordinary Differential Equations. Dover Publications, INC. New York, 1926.

[Kat82] N. Katz. A Conjecture in the Arithmetic Theory of Differential Equations. In Bulletin de la Société Mathématique de France, 110: 515-534, 1982.

[Ore32] O. Ore. Formale Theorie der linearen Differentialgleichungen (Zweiter Teil), J. für d. Reine u. angew. Math., 168: 233-252, 1932.

[Pfl97] E. Pflügel. An Algorithm for Computing Exponential Solutions of First Order Linear Differential Systems. In Proceedings of ISSAC’97, ACM, New York, 1997.

[Put95] M. van der Put. Differential Equations in Characteristic p. In Compositio Mathematica, 97: 227-251, 1995.

[Put96] M. van der Put. Reduction Modulo p of Differential Equations. In Indag. Math.N.S, 7(3): 367-387, 1996.

[Put01] M. van der Put. Grothendieck’s Conjecture for the Risch Equation y′ = ay+b. In Indag. Mathem., 12(1): 113-124, 2001.

[PS03] M. van der Put, M. F. Singer. Galois Theory of Linear Differential Equations. In Grundlehren der mathematischen Wissenschaften, Vol. 328, Springer, 2003.

[Sin96] M. F. Singer. Testing Reducibility of Linear Differential Operators: A Group Theoretic Perspective. In Journal of Appl. Alg. in Eng. Comm. and Comp. (AAECC), 7(2): 77-104, 1996.

[UW96] F. Ulmer, J.-A. Weil. Note on Kovacic’s algorithm. In Journal Of Symbolic Computation, 22: 179-200, 1996.