On Symmetric Powers of Differential Operators

Manuel Bronstein Thom Mulders Institute for Scientific Computation

ETH Zentrum IFW CH-8092 Zürich

{bronstein,mulders}@inf.ethz.ch
http://www.inf.ethz.ch/personal/{bronstein,mulders}/

Jacques-Arthur Weil Département de Mathématiques

Faculté des Sciences, 123 Av. Albert Thomas F-87060 Limoges

                    weil@unilim.fr
    http://Medicis.Polytechnique.fr/gage/weil.html

Abstract

We present alternative algorithms for computing symmetric powers of linear ordinary differential operators. Our algorithms are applicable to operators with coefficients in arbitrary integral domains and become faster than the traditional methods for symmetric powers of sufficiently large order, or over sufficiently complicated coefficient domains. The basic ideas are also applicable to other computations involving cyclic vector techniques, such as exterior powers of differential or difference operators.

Introduction

Let R be an integral domain of characteristic 0, D be a derivation on R, K be the quotient field of R and R[∂; D] be the corresponding ring of linear differential operators with coefficients in R. Let L = ∑_{i=0}^{n} a_{i}∂^{i} ∈ R[∂; D] with n > 0 and a_{n} ≠ 0, Y_{0}, …, Y_{n−1} be indeterminates, and consider the extension of D to K[Y_{0}, …, Y_{n−1}] given by DY_{i} = Y_{i+1} for 0 ≤ i < n−1 and DY_{n−1} = −a_{n}^{−1} ∑_{i=0}^{n−1} a_{i}Y_{i}. For any in-

teger m > 0, the m^{th} symmetric power of L, denoted L^{◯ s ^{m}}, is defined to be a nonzero element of R[∂; D] of minimal degree such that L^{◯ s ^{m}}(Y_{0}^{m}) = 0. Let H_{m} ⊂ K[Y_{0}, …, Y_{n−1}] be the set of homogeneous polynomials of degree m. Since H_{m} is a vector space of dimension N = \binom{n+m−1}{n−1} over K and is closed under D, it follows that the elements Y_{0}^{m}, D (Y_{0}^{m}) , …, D^{N} (Y_{0}^{m}) of H_{m} are linearly dependent over K, hence that L^{◯ s ^{m}} always exists and has degree at most N. Symmetric powers have many applications within algorithms for solving linear differential equations:

• the factorization patterns of L^{◯ s ^{m}} determine the structure of the Galois groups of second and third order equations, and provide necessary and sufficient conditions for the existence of Liouvillian solutions [13];

• the radical solutions of L^{◯ s ^{m}} provide explicitely the Liouvillian solutions of second and third order equations [14]; • the rational solutions of L^{◯ s ^{m}} provide explicitely the

Liouvillian solutions of completely reducible second order equations [16], and are needed as a first step in computing the Liouvillian solutions of higher-order equations [15];

For computing invariants of operators, the algorithms presented in this paper provide a rational alternative to the local method of [18].

The basic algorithm [11, 13] for computing L^{◯ s ^{m}} is to look for a linear dependence over K between Y_{0}^{m}, D (Y_{0}^{m}) , …, D^{i} (Y_{0}^{m}) for increasing i’s until i = N, the coefficients of the first dependence found being the coefficients of L^{◯ s ^{m}}. Expressing each D^{j} (Y_{0}^{m}) in the monomial basis of H_{m} converts the problem of looking for a linear dependence over K to finding the kernel of an N × (i + 1) matrix with coefficients in K. This kernel computation becomes however rapidly expensive as n and m increase, so we present in this paper algorithms for computing L^{◯ s ^{m}} that either avoid this linear algebra step (when L has order 2), or reduce its cost through preprocessing of the corresponding matrix (when L has arbitrary order). In addition, our algorithms can be carried out in a fraction–free fashion over an arbitrary integral domain rather than over its quotient field, avoiding gcd computations and allowing the use of efficient fraction–free linear algebra algorithms. Notation: for any x ∈ R and integer m ≥ 0, x\underline{m} is the descending factorial of x of degree m, i.e. x\underline{m} = ∏_{i=0}^{m−1}(x−i). Note that m\underline{m} = m\underline{m−1} = m! and m\underline{m+1} = 0. For any rational number q, ⌊q⌋ is the largest integer m such that m ≤ q.

1 An iteration for second-order operators

In the case where deg(L) = 2, it is known that L^{◯ s ^{m}} has degree m + 1. Furthermore, in the appropriate basis of H_{m} the matrix of Y_{0}^{m}, D (Y_{0}^{m}) , …, D^{m} (Y_{0}^{m}) is upper triangular, so the kernel computation can be avoided altogether, resulting in the following iteration for computing L^{◯ s ^{m}}.

Theorem 1 Let L = ∂^{2} + a∂ + b ∈ R[∂; D], m > 0 be an integer and consider the sequence given by L_{0} = 1, L_{1} = ∂ and

L_{i+1} = (∂ + ia)L_{i} + i(m − (i − 1))bL_{i−1}

for 0 < i ≤ m. Then, L_{m+1} = L^{◯ s ^{m}}.

Proof. Note that deg(L_{i}) = i for 0 ≤ i ≤ m + 1 and that each L_{i} is monic. Let Y_{0}, Y_{1} be indeterminates and consider the extension of D to K[Y_{0}, Y_{1}] given by DY_{0} = Y_{1} and DY_{1} = −aY_{1} − bY_{0}. We show by induction on i that

L_{i}(Y_{0}^{m}) = m\underline{i}Y_{0}^{m−i}Y_{1}^{i}  (1)

for 0 ≤ i ≤ m + 1. We have L_{0}(Y_{0}^{m}) = Y_{0}^{m} and L_{1}(Y_{0}^{m}) = ∂(Y_{0}^{m}) = mY_{0}^{m−1}Y_{1}, so let 1 ≤ j ≤ m and suppose that (1) holds for 0 ≤ i ≤ j. Then,

L_{j+1}(Y_{0}^{m}) = (∂ + ja)L_{j}(Y_{0}^{m}) + j(m − (j − 1))bL_{j−1}(Y_{0}^{m})
= m\underline{j}∂(Y_{0}^{m−j}Y_{1}^{j}) + jam\underline{j}Y_{0}^{m−j}Y_{1}^{j}
+j(m − (j − 1))bm\underline{j−1}Y_{0}^{m−(j−1)}Y_{1}^{j−1}
= m\underline{j}(m − j)Y_{0}^{m−(j+1)}Y_{1}^{j+1}
+jm\underline{j}Y_{0}^{m−j}Y_{1}^{j−1} (−aY_{1} − bY_{0})
+jam\underline{j}Y_{0}^{m−j}Y_{1}^{j} + jbm\underline{j}Y_{0}^{m−(j−1)}Y_{1}^{j−1}
= m\underline{j+1}Y_{0}^{m−(j+1)}Y_{1}^{j+1}

which implies that (1) holds for j + 1. Therefore,

L_{m+1}(Y_{0}^{m}) = m\underline{m+1}Y_{0}^{−1}Y_{1}^{m} = 0

which implies that L_{m+1} = L^{◯ s ^{m}} since it has degree m + 1. □

When L = p∂^{2} + a∂ + b and p is not a unit in R, we want to avoid going to the quotient field in order to avoid computing gcd’s. While we can replace L_{i} by p^{i−1}L_{i} in the iteration of Theorem 1 and obtain a sequence of operators in R[∂; D], those operators may have increasing contents, in a similar fashion to what happens in polynomial remainder sequences. Large parts of those contents can however be predicted, yielding a fraction–free iteration over R, which we present in the rest of this section.

Lemma 1 If R is a unique factorization domain then

r | q | gcd(p, Dp) ⟹ qr | pDr

for any p, q, r ∈ R.

Proof. We can assume that p ≠ 0, the result being trivial otherwise. Suppose first that r is irreducible. If r | Dr, then qr | pDr since q | p, so suppose that r |̸ Dr and write p = r^{n}h where n > 0 and r |̸ h. Then, Dp = nr^{n−1}hDr + r^{n}Dh, so r |̸ hDr implies that r^{n} |̸ Dp, hence that r^{n} |̸ q since q | Dp. Therefore, r | (p/q), which implies that qr | p, hence that qr | pDr. Let now r_{1}, r_{2} be such that r_{i} | q and qr_{i} | pDr_{i} for i = 1, 2. Then,

\frac{p}{q} \frac{D(r_{1}r_{2})}{r_{1}r_{2}} = \frac{p}{q} (\frac{Dr_{1}}{r_{1}} + \frac{Dr_{2}}{r_{2}} ) = \frac{pDr_{1}}{qr_{1}} + \frac{pDr_{2}}{qr_{2}} ∈ R

so qr_{1}r_{2} | pD(r_{1}r_{2}) and the lemma follows. □

Theorem 2 Suppose that R is a gcd domain. Let L = p∂^{2}+ a∂ + b ∈ R[∂; D] where p ≠ 0, p_{e} = gcd(p, Dp, a − Dp), p_{o} = gcd(p_{e}, bp/p_{e}), m > 1 be an integer and consider the sequence given by L_{0} = 1, L_{1} = ∂, L_{2} = p∂^{2} + a∂ + mb,

L_{3} = (\frac{p}{p_{e}} ∂ + 2 \frac{a − Dp}{p_{e}} + \frac{Dp}{p_{e}}) L_{2} + 2(m − 1)b \frac{p}{p_{e}} L_{1} ,(2)

L_{2i} = (\frac{p}{p_{o}} ∂ + (2i − 1) \frac{a − Dp}{p_{o}} + \frac{Dp}{p_{o}} + (i − 1) \frac{pDp_{e}}{p_{o}p_{e}}

+(i − 2) \frac{pDp_{o}}{p_{o}^{2}} ) L_{2i−1}

+(2i − 1)(m − (2i − 2)) \frac{bp}{p_{o}p_{e}} L_{2i−2}  (3)

and

L_{2i+1} = (\frac{p}{p_{e}} ∂ + 2i \frac{a − Dp}{p_{e}} + \frac{Dp}{p_{e}} + (i − 1) \frac{pDp_{e}}{p_{e}^{2}}

+(i − 1) \frac{pDp_{o}}{p_{e}p_{o}} ) L_{2i}

+2i(m − (2i − 1)) \frac{bp}{p_{o}p_{e}} L_{2i−1}  (4)

for 2 ≤ i ≤ ⌊(m + 1)/2⌋. Then, L_{m+1} = L^{◯ s ^{m}}.

Proof. Define Q_{0} = L_{0} = 1, Q_{1} = L_{1} = ∂, and

Q_{j} = \frac{p_{e}^{⌊(j−1)/2⌋}p_{o}^{⌊j/2⌋−1}}{p^{j−1}} L_{j}  (5)

for 2 ≤ j ≤ m + 1. We first prove that

Q_{j+1} = (∂ + j \frac{a}{p}) Q_{j} + j(m − (j − 1)) \frac{b}{p} Q_{j−1} (6)

for 0 < j ≤ m. We have

Q_{2} = \frac{1}{p} L_{2} = ∂^{2} + \frac{a}{p} ∂ + m \frac{b}{p} = (∂ + \frac{a}{p}) Q_{1} + m \frac{b}{p} Q_{0}

so (6) holds for j = 1. Using (2) we get

Q_{3} = \frac{p_{e}}{p^{2}} L_{3} = (\frac{1}{p} ∂ + 2 \frac{a − Dp}{p^{2}} + \frac{Dp}{p^{2}}) L_{2} + 2(m − 1) \frac{b}{p} L_{1}

= \frac{1}{p} ∂(pQ_{2}) + (2 \frac{a − Dp}{p} + \frac{Dp}{p}) Q_{2} + 2(m − 1) \frac{b}{p} Q_{1}

= (∂ + 2 \frac{a}{p}) Q_{2} + 2(m − 1) \frac{b}{p} Q_{1}

so (6) holds for j = 2. Let 3 ≤ j ≤ m. If j is odd, then j + 1 = 2i for some i ≥ 2, so using (3) and (5) we get

Q_{j+1} = \frac{p_{e}^{i−1}p_{o}^{i−1}}{p^{j}} L_{j+1} = \frac{p_{e}^{i−1}p_{o}^{i−1}}{p^{j}} L_{2i} = \frac{p_{e}^{i−1}p_{o}^{i−1}}{p^{j}} \{

(\frac{p}{p_{o}} ∂ + j \frac{a − Dp}{p_{o}} + \frac{Dp}{p_{o}} + (i − 1) \frac{pDp_{e}}{p_{o}p_{e}}

+(i − 2) \frac{pDp_{o}}{p_{o}^{2}} ) \frac{p^{j−1}}{p_{e}^{i−1}p_{o}^{i−2}} Q_{j}

+j(m − (j − 1)) \frac{bp}{p_{o}p_{e}} \frac{p^{j−2}}{p_{e}^{i−2}p_{o}^{i−2}} Q_{j−1}\}

= \frac{p_{e}^{i−1}p_{o}^{i−2}}{p^{j−1}} ∂ (\frac{p^{j−1}}{p_{e}^{i−1}p_{o}^{i−2}} Q_{j}) + (j \frac{a − Dp}{p} + \frac{Dp}{p}

+(i − 1) \frac{Dp_{e}}{p_{e}} + (i − 2) \frac{Dp_{o}}{p_{o}} ) Q_{j} + j(m − (j − 1)) \frac{b}{p} Q_{j−1}

= (∂ + j \frac{a}{p}) Q_{j} + j(m − (j − 1)) \frac{b}{p} Q_{j−1}

so (6) holds for j odd. A similar calculation using (4) and (5) shows that (6) holds also for j even. Since Q_{0} = 1, Q_{1} = ∂ and (6) holds for 0 < j ≤ m, Theo-

rem 1 applied to K implies that Q_{m+1} = (L/p)^{◯ s ^{m}} hence that L_{m+1} = L^{◯ s ^{m}}. □

Note that since p_{o} | p_{e} | gcd(p, Dp), Lemma 1 implies that if R is a unique factorization domain, then the coefficients of the iteration of Theorem 2 are in R, i.e. that all the quotients are exact.

Furthermore, if follows from (6) that each Q_{j} has leading coefficient 1, so (5) implies that the leading coefficient of L^{◯ s ^{m}} is

c_{m} = \frac{p^{m}}{p_{e}^{⌊m/2⌋}p_{o}^{⌊(m+1)/2⌋−1}} .

If R is a unique factorization domain, then c_{m} ∈ R and every irreducible factor of c_{m} must divide p, which implies in particular that if R is of the form k[x] with Dx = 1, then all the singularities of L^{◯ s ^{m}} are singularities of L. This last point follows also from the iteration of Theorem 1.

2 A general fraction–free method

Let L = ∑_{i=0}^{n} a_{i}∂^{i} ∈ R[∂; D] with n > 0 and a_{n} ≠ 0, Y_{0}, …, Y_{n−1} be indeterminates and consider the derivation Δ = a_{n}D on K[Y_{0}, …, Y_{n−1}]. Then, ΔY_{i} = a_{n}Y_{i+1} for 0 ≤ i < n − 1 and ΔY_{n−1} = − ∑_{i=0}^{n−1} a_{i}Y_{i}, which implies that R[Y_{0}, …, Y_{n−1}] is closed under Δ.

Lemma 2 Let m > 0 be an integer,

w_{0} = Y_{0}^{m} and w_{i+1} = Δw_{i} − iD(a_{n})w_{i} for i ≥ 0 . (7)

If ∑_{i=0}^{M} c_{i}w_{i} = 0 for some integer M ≥ 0 and c_{0}, …, c_{M} ∈ R, then

(∑_{i=0}^{M} c_{i}a_{n}^{i}∂^{i}) (Y_{0}^{m}) = 0 .

Proof. We first prove by induction on i that

w_{i} = a_{n}^{i}D^{i} (Y_{0}^{m})

for i ≥ 0. We have w_{0} = Y_{0}^{m} = a_{n}^{0}D^{0} (Y_{0}^{m}), so suppose that w_{i} = a_{n}^{i}D^{i} (Y_{0}^{m}) for some i ≥ 0. Then,

w_{i+1} = Δw_{i} − iD(a_{n})w_{i}
= a_{n}D (a_{n}^{i}D^{i} (Y_{0}^{m})) − iD(a_{n})a_{n}^{i}D^{i} (Y_{0}^{m})
= a_{n}^{i+1}D^{i+1} (Y_{0}^{m}) + ia_{n}^{i}D(a_{n})D^{i} (Y_{0}^{m})
−iD(a_{n})a_{n}^{i}D^{i} (Y_{0}^{m}) = a_{n}^{i+1}D^{i+1} (Y_{0}^{m}) .

It follows that

(∑_{i=0}^{M} c_{i}a_{n}^{i}∂^{i}) (Y_{0}^{m}) = ∑_{i=0}^{M} c_{i}a_{n}^{i}D^{i} (Y_{0}^{m}) = ∑_{i=0}^{M} c_{i}w_{i} = 0 .

□

It follows that L^{◯ s ^{m}} can be obtained from the first linear dependence found over R between w_{0}, w_{1}, …, w_{i} for increasing i until i = N, which is equivalent to finding the kernel of an N × (i + 1) matrix with coefficients in R. However, as in the case of second-order operators, the iteration (7) produces increasing contents in the rows of the resulting matrix.

Lemma 3 Let (w_{i})_{i≥0} be given by (7) and Y ∈ K[Y_{0}, …, Y_{n−1}] be a primitive monomial of total degree m. Then for any i ≥ 0, the coefficient of Y in w_{i} is divisible by a_{n}^{P (Y,i)} where

P (Y, i) = \{ \begin{matrix}F (Y) & if i < F (Y); \\ i & if F (Y) ≤ i ≤ F (Y) + G(Y); \\ F (Y) + G(Y) & if i > F (Y) + G(Y), \end{matrix}

F (Y) = ∑_{k=1}^{n−1} k deg_{Y_{k}}(Y)

and

G(Y) = n − 1 − max{k | deg_{Y_{k}}(Y) > 0} .

Proof. By induction on i. We have w_{0} = Y_{0}^{m} and F (Y_{0}^{m}) = 0, so the lemma holds for i = 0. Suppose now that i ≥ 0 and that the coefficient of Y in w_{i} is divisible by a_{n}^{P (Y,i)} for any primitive monomial in Y ∈ H_{m}. Since w_{i} ∈ H_{m}, we can write w_{i} = ∑_{Y} b_{Y} Y, where b_{Y} ∈ R and the sum is taken over all primitive monomials of total degree m. We have

w_{i+1} = Δw_{i} − iD(a_{n})w_{i}
= ∑_{Y}(a_{n}D(b_{Y})Y − iD(a_{n})b_{Y} Y + a_{n}b_{Y} D(Y)) .

It suffices to show that for all primitive monomials Y and Ỹ of total degree m, a_{n}^{P (Ỹ ,i+1)} divides the coefficient of Ỹ in a_{n}D(b_{Y})Y −iD(a_{n})b_{Y} Y +a_{n}b_{Y} D(Y). Writing b_{Y} = ba_{n}^{P (Y,i)} we get

a_{n}D(b_{Y}) − iD(a_{n})b_{Y} =
D(b)a_{n}^{P (Y,i)+1} + (P (Y, i) − i)ba_{n}^{P (Y,i)}D(a_{n}) . (8)

When P (Y, i) ≠ i we have P (Y, i + 1) = P (Y, i). From this and (8) we see that a_{n}^{P (Y,i+1)} divides a_{n}D(b_{Y}) − iD(a_{n})b_{Y}. For a primitive monomial Y ∈ H_{m} and 0 ≤ j, k < n such that deg_{Y_{j}} (Y) > 0, we define Y^{[j,k]} to be the primitive monomial Y Y_{k}/Y_{j}. Whenever 0 ≤ j < n−1 and deg_{Y_{j}} (Y) > 0, we have

F (Y^{[j,j+1]}) = F (Y) + 1  (9)

and when 0 ≤ j < n and deg_{Y_{n−1}}(Y) > 0, then

F (Y^{[n−1,j]}) = F (Y) − (n − 1) + j.  (10)

Writing Y = ∏_{j=0}^{n−1} Y_{j}^{e_{j}}, we have

a_{n}b_{Y} D(Y) = a_{n}b_{Y} ∑_{j=0}^{n−2} e_{j}Y^{[j,j+1]} − e_{n−1}b_{Y} ∑_{j=0}^{n−1} a_{j}Y^{[n−1,j]}.

For Ỹ = Y^{[j,j+1]} we have G(Ỹ) ≤ G(Y). From this and (9) it follows easily that P (Ỹ , i+1) ≤ P (Y, i)+1 and so a_{n}^{P (Ỹ ,i+1)} divides a_{n}b_{Y}.

For Ỹ = Y^{[n−1,j]} we have G(Ỹ) ≤ n − 1 − j. From this and (10) it follows that P (Ỹ , i + 1) ≤ F (Ỹ) + G(Ỹ) ≤ F (Y) ≤ P (Y, i), so a_{n}^{P (Ỹ ,i+1)} divides b_{Y}. □ From the previous lemma we see that a_{n}^{F (Y )} divides the content of the row corresponding to Y. In fact various entries of the matrix are divisible by even higher powers of a_{n}.

The following matrix shows the powers of a_{n} in the general case of the 2^{nd} symmetric power of an order 3 operator. The numbers on the left are the corresponding F (Y). ∞ means that the corresponding entry is 0.

0 1 2 2 3 4 ( \begin{matrix}0 & ∞ & ∞ & 2 & 2 & 2 & 2 \\ ∞ & 1 & ∞ & 2 & 2 & 2 & 2 \\ ∞ & ∞ & 2 & 2 & 2 & 2 & 2 \\ ∞ & ∞ & 2 & ∞ & 3 & 3 & 3 \\ ∞ & ∞ & ∞ & 3 & 3 & 3 & 3 \\ ∞ & ∞ & ∞ & ∞ & 4 & 4 & 4 \end{matrix} )

Depending on the coefficients of the operator it can happen that even higher powers of a_{n} than predicted by Lemma 3, divide some entry of the matrix as the following example shows.

Example 1 For the 2^{nd} symmetrix power of

(x^{2} + 5)∂^{3} + (5x^{2} − 2x + 25)∂^{2} + (x^{4} − 1)∂ + x^{4} − x^{3} − x − 1

we get the following matrix of powers (of x^{2} + 5) dividing the corresponding entries.

( \begin{matrix}0 & ∞ & ∞ & 2 & 3 & 3 & 3 \\ ∞ & 1 & ∞ & 2 & 4 & 3 & 3 \\ ∞ & ∞ & 2 & 2 & 3 & 4 & 4 \\ ∞ & ∞ & 2 & ∞ & 3 & 4 & 4 \\ ∞ & ∞ & ∞ & 3 & 3 & 4 & 4 \\ ∞ & ∞ & ∞ & ∞ & 4 & 4 & 4 \end{matrix} )

In this section, we have used the operator Δ = a_{n}D to get a matrix with all entries in R. In fact it might be the case that (depending on the coefficients of the operator) it is sufficient to use an operator Δ = cD with c dividing a_{n}, in order to avoid fractions. For example, for third–order operators with a_{3} = p^{3} it suffices to use c = p^{2}. However, notice that in this case R[Y_{0}, …, Y_{n−1}] is not closed under Δ. As will be seen in the next section, the extra contents brought in by using a_{n}D instead of cD can be removed through an adequate preprocessing of the matrix.

In order to avoid fractions one might also use the original operator D and multiply a column by a_{n} only when fractions do appear.

3 Preprocessing the matrix

When applying the fraction–free method described in the previous section, we have to compute a linear dependence between the w_{j}. For this we can use any method from linear algebra and in particular fraction–free Gaussian elimination ([5] §9.3) or modular methods [8] are suitable.

Choosing some order of the monomials of total degree m in Y_{0}, …, Y_{n−1} we can write the coefficients of w_{j} in a matrix A = (a_{ij}) (i.e. a_{ij} is the coefficient of the ith monomial in w_{j}). Computing a linear dependence between the w_{j} and between the columns of A are now equivalent problems.

As shown in the previous section, the entries in A can be divisible by high powers of a_{n}. In this section we will show how we can preprocess the matrix A in order to make the computation of a linear dependence between its columns more efficient.

It is clear that the kernel of a matrix does not change by multiplying the rows of that matrix by scalars. By multiplying the columns of a matrix by scalars the kernel does change, but there is an easy relationship between the kernels. Our preprocessing strategy consists of multiplying the rows and columns of the matrix A by powers of a_{n} in order to decrease the powers of a_{n} dividing the entries of the matrix. Let B = (b_{ij}) be a matrix, having entries in ℕ ∪ {∞}, such that b_{ij} = ∞ when a_{ij} = 0 and a_{n}^{b_{ij}} divides a_{ij} when a_{ij} ≠ 0. For this we can use for example the powers predicted by Lemma 3 or, when R is a unique factorization domain, we can use the true exponents of a_{n} in the entries of A. Multiplication of the ith row (resp. jth column) of A by a_{n}^{k} can now be translated to adding k to each entry of the ith row (resp. jth column) of B. We can state the preprocessing problem now as follows:

Let B be an n × m matrix with coefficients in ℕ∪{∞}. Find r_{i} ∈ ℤ (1 ≤ i ≤ n) and c_{j} ∈ ℤ (1 ≤ j ≤ m) such that the following holds:

1. b_{ij} + r_{i} + c_{j} ≥ 0 for all i, j.

2. ∑(b_{ij}+r_{i}+c_{j}) should be minimal, where the sum is taken over all i, j such that b_{ij} ≠ ∞ (i.e. the remaining powers of a_{n} should be minimal).

When we omit the condition that the variables r_{i} and c_{j} should have integer values, this is an example of a linear programming problem. For details on linear programming we refer to [9]. It is clear that this problem has a solution.

We can write the constraints as Cx ≤ b, where C is a matrix whose columns are labeled by r_{i} and c_{j} and whose rows are labeled by b_{ij} ≠ ∞, and b is a vector whose entries are labeled by (and are equal to) b_{ij} ≠ ∞. Then all entries in the b_{ij}’s row of C are 0, except for the entries in columns r_{i} and c_{j}, which are equal to −1. We see that C is the negative of the incidence matrix of the bipartite graph G, whose vertices are labeled by r_{i} and c_{j} and which has an edge from r_{i} to c_{j} when b_{ij} ≠ ∞.

An important class of matrices in linear programming are the totally unimodular matrices. One of the characterizations of these matrices is that all minors equal −1, 0 or 1. It is well known that the incidence matrix of a bipartite graph is totally unimodular, so C is totally unimodular. From this it follows that a solution to our problem has integer components.

In our problem the variables r_{i} and c_{j} are free variables (i.e. their value may be positive or negative) and the constraints are inequalities. If we want to use the standard simplex method to solve this problem we first have to transform it to a problem having only non-negative variables (i.e. their values are non-negative), which doubles the number of variables. If N is the number of constraints in the original problem then solving the new problem using the simplex method would need pivoting in an (N+1)×(N+2(n+m)+1) matrix.

The dual problem of the original problem only has N non-negative variables and n + m equalities. So solving the problem via the dual problem would need pivoting in an (n + m + 1) × (N + 1 + n + m) matrix (we not only have to solve the dual problem but also the dual of the dual which equals the original problem).

Since N = O(nm), the solution via the dual problem seems to be the most efficient one.

Note that the dual problem is essentially the Hitchcock problem (or transportation problem) for which efficient algorithms do exist ([6]).

As before, we can see that the matrix corresponding to the dual problem is totally unimodular. Next we will show that during the process of pivoting the matrix will stay totally unimodular.

Lemma 4 Let S = (s_{ij}) be a totally unimodular n × m matrix, 1 ≤ i_{0} ≤ n, 1 ≤ j_{0} ≤ m such that s_{i_{0}j_{0}} ≠ 0. Let T = (t_{ij}) be the matrix we get by pivoting matrix S using s_{i_{0}j_{0}} as pivot (i.e. clearing the j_{0}^{th} column of S). Then T is totally unimodular.

Proof. Since s_{i_{0}j_{0}} = ±1 and by permuting the rows and columns of S we see that it suffices to prove the lemma for i_{0} = j_{0} = 1 and s_{11} = 1. Notice that

t_{ij} = \{ \begin{matrix}s_{ij} & if i = 1 \\ s_{11}s_{ij} − s_{i1}s_{1j} & if i > 1 \end{matrix}

When 1 ≤ i_{1} < ⋯ < i_{k} ≤ n and 1 ≤ j_{1} < ⋯ < j_{k} ≤ m, we denote by ([i_{1}, …, i_{k}], [j_{1}, …, j_{k}])_{U} the submatrix of matrix U consisting of rows i_{1}, ⋯, i_{k} and columns j_{1}, …, j_{k}.

Since t_{i1} = 0 for i > 1, we have when i_{1} > 1

det (([i_{1}, …, i_{k}], [1, j_{2}, …, j_{k}])_{T}) = 0 .

When i_{1} = 1 then

det (([i_{1}, …, i_{k}], [j_{1}, …, j_{k}])_{T}) =
det (([i_{1}, …, i_{k}], [j_{1}, …, j_{k}])_{S})

since ([i_{1}, …, i_{k}], [j_{1}, …, j_{k}])_{T} is obtained from ([i_{1}, …, i_{k}], [j_{1}, …, j_{k}])_{S} by adding or subtracting the first row from some other rows.

When i_{1} > 1 and j_{1} > 1 then Sylvester’s identity ([5], §9.3) shows that

det (([i_{1}, …, i_{k}], [j_{1}, …, j_{k}])_{T})

= s_{11}^{k−1} det (([1, i_{1}, …, i_{k}], [1, j_{1}, …, j_{k}])_{S})

= det (([1, i_{1}, …, i_{k}], [1, j_{1}, …, j_{k}])_{S}) .
□

This lemma shows that during the process of pivoting the matrix is always totally unimodular, which implies that the pivot (which is a 1-minor) always equals ±1. This means that only integers appear during the pivoting and this fact can be used to produce a faster implementation of the algorithm.

Example 2 The solution to the linear program corresponding to the matrix B of Example 1 is

(r_{1}, …, r_{6}, c_{1}, …, c_{7}) =

(−3, −3, −4, −4, −4, −4, 3, 2, 2, 2, 1, 0, 0) .

Adding r_{i} + c_{j} to b_{ij} we get:

( \begin{matrix}0 & ∞ & ∞ & 1 & 1 & 0 & 0 \\ ∞ & 0 & ∞ & 1 & 2 & 0 & 0 \\ ∞ & ∞ & 0 & 0 & 0 & 0 & 0 \\ ∞ & ∞ & 0 & ∞ & 0 & 0 & 0 \\ ∞ & ∞ & ∞ & 1 & 0 & 0 & 0 \\ ∞ & ∞ & ∞ & ∞ & 1 & 0 & 0 \end{matrix} )

The sum of all entries in B is 78, while the sum of all entries in the optimized matrix is only 7.

If one has some kind of factorization of a_{n} one can also use this optimization technique to all factors seperately. So depending on the algorithms available for R (factorization, gcd computations or none) one can use a complete factorization, balanced factorization or no factorization of a_{n}. Using this preprocessing at each factor compensates for using a_{n}D as derivation (instead of cD for some factor c of a_{n}) when generating the matrix.

4 Implementation issues and timings

We have implemented the above algorithms in Σ^{IT} 0.1.8b [2] and in MAPLE V.3 [4]. In order to benefit fully from the advantages of the fraction–free methods, the following points must be considered:

• for second order operators, all the exact quotients appearing in Theorem 2 should be precomputed before starting the iteration;

• because of the high powers (Lemma 3) of the leading coefficient a_{n} of L appearing in the iteration (7), the elements of R should be manipulated in the form a_{n}^{d}b where d ∈ ℕ and b ∈ R, so that the powers a_{n}^{d} are not expanded;

• since the fraction–free algorithms work over arbitrary integral domains, when R = k[x] and k is the quotient field of an integral domain S, then cL ∈ S[x][∂; D] for some nonzero c ∈ S, so we can take advantage of the faster arithmetic in S[x] and compute (cL)^{◯ s ^{m}} rather than L^{◯ s ^{m}}. An important case is R = ℚ[x], since multiplication in ℤ[x] has better asymptotic complexity than in ℚ[x]^{1}.

• when R is a univariate polynomial ring (with arbitrary derivation), then the preprocessing of the previous section should be done for each balanced factor of a_{n}.

• if a ring-homomorphism (or even ℤ-module homomorphism) σ : R → 𝔽 is available for a finite field 𝔽, then the linear independence over 𝔽 of w_{0}^{σ}, w_{1}^{σ}, …, w_{i}^{σ} can be checked while the w_{i}’s are being computed. Let d be the smallest integer such that w_{0}^{σ}, …, w_{d}^{σ} are linearly dependent over 𝔽. Then, deg( L^{◯ s ^{m}}) ≥ d and for suitable R, 𝔽 and σ, that degree is d with high probability [10]. We can then apply the preprocessing of Section 3 to the matrix A whose columns consist of w_{0}, …, w_{d} and continue the iteration (7) beyond w_{d} only if ker A = 0. When R = ℤ[x], we can also attempt to directly lift the linear dependence [8]. Incidentally, this heuristic either proves that deg( L^{◯ s ^{m}}) = \binom{n+m−1}{n−1} or indicates that it is stricly smaller with high probability, an important information for the algorithm of [15].

We present in the rest of this section empirical timings of the various methods presented in this paper. In addition, since our MAPLE implementation of the fractional kernel method takes advantage of the many cancellations appearing during the linear algebra, we also give the timings obtained with version “Nov.21 1996” of the diffop package of [17]. All the timings in this section were obtained on a single CPU DEC Alpha 500/333 with 256Mb main memory running Digital Unix V4.0a. All timings are given in seconds of CPU, garbage collection included, and a ∗ indicates that the computation was aborted after one hour of CPU. In order to be able to verify the results, we use examples with known Galois groups and/or invariants from the literature.

1The Karatsuba multiplication algorithm is not effective in ℚ[x] since multiplication and addition have the same complexity in ℚ.

4.1 Timings for second-order operators

For rational coefficients, we use the following operators from [16]:

D_{2} = D^{2} − \frac{2}{2x − 1} D + \frac{3(2x − 1)^{2}(x^{4} − 2x^{3} + x + 1)}{16x^{2}(x − 1)^{2}(x^{2} − x − 1)^{2}} ,

D_{3} = D^{2}− \frac{2}{2x − 1} D+ \frac{(27x^{4} − 54x^{3} + 5x^{2} + 22x + 27)(2x − 1)^{2}}{144x^{2}(x − 1)^{2}(x^{2} − x − 1)^{2}} ,

A_{4} = D^{2} + \frac{3}{16x^{2}} + \frac{2}{9(x − 1)^{2}} − \frac{3}{16x(x − 1)}

and for an example with algebraic numbers in its coefficients, the following non-completely reducible operator from [1]:

L\sqrt{−222} = (x^{12} + 2x^{10} + \frac{151}{3} x^{8} + \frac{296}{3} x^{6} + \frac{5920}{9} x^{4}

+ \frac{10952}{9} x^{2} + \frac{5476}{9}) D^{2} − \frac{104}{25} x^{10}

− (\frac{274}{25} − \frac{22}{15} \sqrt{−222}) x^{8} − (\frac{7754}{75} − \frac{68}{15} \sqrt{−222}) x^{6}

− (\frac{11248}{75} − \frac{194}{15} \sqrt{−222}) x^{4}

− (\frac{29452}{75} − \frac{296}{5} \sqrt{−222}) x^{2} + \frac{10952}{5} + \frac{148}{3} \sqrt{−222} .

Table 1 contains the times needed for computing symmetric powers by the various methods of this paper: Iteration ℤ[x] is the fraction–free iteration of Theorem 2, Iteration ℚ(x) is the fractional iteration of Theorem 1, Kernel ℤ[x] is the fraction–free kernel method of Section 2, and Kernel ℚ(x) is the classical kernel method^{2}. We see from Table 1 that the iteration methods of Section 1 are faster than the kernel methods, and that the fraction–free iteration is faster than the fractional one. The fraction–free kernel method of Section 2 is not always faster than the fractional kernel method, which is to be expected since it reduces the cost of the linear algebra, but the matrices generated in the case of second order operators are triangular, so the linear algebra cost is negligible compared to the simplex preprocessing and overhead. As will be seen in the next section, the fraction– free kernel method becomes better when the linear algebra cost increases, i.e. for operators of sufficiently high order, or with sufficiently complicated coefficients. We remark also that the fraction–free iteration together with a fraction–free rational kernel finder makes the algorithm of [16] quite efficient: Darboux curves for D_{2}, D_{3} and A_{4} were computed in Σ^{IT} in 1.2, 1.4 and 1.1 seconds respectively, without a priori knowledge of their degrees, while the MAPLE differential equation solver was only able to solve D_{3}.

Table 2 contains the times needed for computing symmetric powers of L\sqrt{−222} by the fraction–free (ffree) and fractional (frac) iterations of Section 1. It shows as ex-

Σ^{IT}  MAPLE
ffree frac ffree frac diffop
L\sqrt{−222}^{◯ s ^{2}} 1.3  0.5  0.3 0.3  1.7
L\sqrt{−222}^{◯ s ^{3}} 2.2 71.4 0.6 0.6  4.9
L\sqrt{−222}^{◯ s ^{4}} 3.8 286.3 0.9 1.3  13.2
L\sqrt{−222}^{◯ s ^{5}} 5.9 2252 1.6 2.5  34.5
L\sqrt{−222}^{◯ s ^{6}} 8.7  ∗  2.1 4.4  69.2
L\sqrt{−222}^{◯ s ^{12}} 60  ∗  10 41.5  ∗

Table 2: second order operators over ℚ(\sqrt{−222})(x)

pected that fraction–free methods become more advantageous when the cost of computing in the quotient field of R increases. The difference between the two methods is more drastic in Σ^{IT}, which uses the classical Euclidean algorithm for computing greatest common divisors in ℚ(\sqrt{−222})[x], while MAPLE uses better algorithms, making the cost of normalizing quotients less important.

4.2 Timings for higher-order operators

For rational coefficients, we use the following operators from [18]:

F_{36} = D^{3}+ \frac{5(9x^{2} + 14x + 9)}{48x^{2}(x + 1)^{2}} D− \frac{5(81x^{3} + 185x^{2} + 229x + 81)}{432x^{3}(x + 1)^{3}} ,

G_{168} = D^{3} + \frac{7x − 4}{x(x − 1)} D^{2}

+ \frac{2592x^{2} − 2963x + 560}{252x^{2}(x − 1)^{2}} D + \frac{57024x − 40805}{24696x^{2}(x − 1)^{2}} ,

P SL_{3} = xD(xD − \frac{1}{2} )(xD − \frac{1}{4} )(xD + \frac{1}{4} )(xD − \frac{1}{8} )

(xD − \frac{5}{8} )(xD + \frac{1}{8} )(xD + \frac{5}{8} )

−x(xD + \frac{1}{3} )(xD − \frac{1}{3} ) ,

as well as the following operator with larger coefficients from [15]:

L̃ = D^{3} + \frac{128 x^{5} − 98 x^{4} + 108 x^{3} + 378 x^{2} − 486 x + 162}{x(x − 1)p(x)} D^{2}

+ \frac{a(x)}{36 x^{2}(x − 1)^{2}p(x)} D − \frac{b(x)}{72 x(x − 1)^{2}p(x)} ,

where

p(x) = 72 x^{6} − 216 x^{5} + 280 x^{4} − 126 x^{3} + 27 x^{2} + 81 x − 54 ,

a(x) = 2304 x^{8} − 13464 x^{7} + 24872 x^{6} − 24840 x^{5} + 2106 x^{4}
+29565 x^{3} − 26514 x^{2} + 9477 x − 1458 ,

and

b(x) = 4608 x^{6} − 28872 x^{5} + 67240 x^{4} − 115560 x^{3}
+131166 x^{2} − 78165 x + 27135 .

2ℚ(x) is represented as the quotient field of ℤ[x], which is a more efficient type than the quotient field of ℚ[x].

Σ^{IT}  MAPLE
Iteration Kernel Iteration Kernel diffop
ℤ[x] ℚ(x) ℤ[x] ℚ(x) ℤ[x] ℚ(x) ℤ[x] ℚ(x) ℚ(x)
D^{2}^{◯ s ^{4}} 0.04 0.2 0.3 0.3 0.08 0.2 0.4 0.3  1.2
D^{2}^{◯ s ^{12}} 0.8 3.3 6.2 4.9 0.7 1.9 5.8 3.2 37.4
D^{2}^{◯ s ^{20}} 3.5 15.5 42.1 28  3.0 8.1 28.4 14.7 398.8
D^{3}^{◯ s ^{4}} 0.06 0.2 0.3 0.3  0.1 0.1 0.4 0.2  1.3
D^{3}^{◯ s ^{12}} 0.8 3.4 7.1 6.4  0.7 2.1 6.1 3.6  38.1
D^{3}^{◯ s ^{20}} 3.7 16.2 50 31.7 3.2 9.1 30.3 16.6  490
A^{4}^{◯ s ^{6}} 0.02 0.1 0.1 0.3 0.07 0.07 0.4 0.2  0.8
A^{4}^{◯ s ^{12}} 0.1 0.7 0.4 2.1  0.2 0.4 1.6 0.8  4.9
A^{4}^{◯ s ^{20}} 0.3 2.9  3  5.8  0.5 1.4 6.5 2.7  24.4

Table 1: second order operators over ℚ(x)

Finally an example with algebraic numbers in its (small) coefficients:

L\sqrt{2} = (x^{2} + 1)D^{3} − \sqrt{2}D + x .

Table 3 contains the times needed for computing symmetric powers by the fraction–free (ℤ[x]) and fractional (ℚ(x)) kernel methods, while Table 4 contains the times needed for

Σ^{IT}  MAPLE
ℤ[x] ℚ(x) ℤ[x] ℚ(x) diffop
F^{36}^{◯ s ^{5}}  8.3 10.3 10.9  3.9  61.6
F^{36}^{◯ s ^{6}}  31.4 25.7 29.4 14.3  340.7
G^{168}^{◯ s ^{5}} 15.4 26.3 14.0  8.3  60.8
G^{168}^{◯ s ^{6}}  105 79.1 35.6 21.0  500.7
L̃^{◯ s ^{3}}  5.7  12  6.7  6.0  33.6
L̃^{◯ s ^{4}}  61.2 118.2 197.8 144.6 1201.2
L̃^{◯ s ^{5}}  895.6 1384  ∗  2635.3  ∗
P SL^{3}^{◯ s ^{2}} 163.1  ∗  218.5 197.9 200.4

Table 3: higher order operators over ℚ(x)

computing symmetric powers of L\sqrt{2} by the fraction–free (ffree) and fractional (frac) kernel methods.

Σ^{IT}  MAPLE
ffree frac ffree frac diffop
L\sqrt{2}^{◯ s ^{2}} 0.3 1.1 1.4 1.4  5.1
L\sqrt{2}^{◯ s ^{3}} 33.7 ∗  9.6 76.5 81.4
L\sqrt{2}^{◯ s ^{4}}  ∗  ∗  60.2 ∗  ∗

Table 4: higher order operators over ℚ(\sqrt{2})(x)

The Σ^{IT} columns of Table 3 illustrate the cutoff when the fraction–free method starts becoming better than the fractional one in a system with canonical expanded forms for polynomials and fractions such as Σ^{IT} or axi .om [7]: F_{36}^{◯ s ^{m}} and G_{168}^{◯ s ^{m}} have orders smaller than the generic (m+1)(m+2)/2 and their coefficients are rather small, so the fractional method remains faster. On the other hand, L̃ and its powers have quite larger coefficients, and the fraction– free method is then faster, as well as on the higher-order example P SL_{3}. This cutoff is further away on MAPLE because of two characteristics of the arithmetic in MAPLE:

• the preprocessing step of Section 3 consists of a simplex on a matrix of machine-integers. Machine integers are however not accessible in the MAPLE user language, causing the preprocessing step to take a more significant portion of the computing time than in compiled languages, thereby pushing the cutoff further away;

• the normal function in MAPLE keeps the numerators and denominators of elements of ℚ(x) factorized, and this allows an appropriately coded Gaussian elimination to take advantage of the many cancellations that happens on the matrices that arise in the symmetric powers computations, thereby giving an extra advantage to the fractional method; At the same time, the intermediate results tend to become expanded during the fraction-free elimination, giving it an additional disadvantage;

Both of the above points disappear when the constant field is a proper extension of ℚ, as illustrated by Table 4, where the fraction–free method is significantly faster than the fractional one on all systems.

Conclusions

For second order operators, the fraction–free iteration method of Theorem 2 is the most efficient way of computing symmetric powers, and it turns the algorithm of [16] into a very efficient second order linear ordinary equation solver. For higher-order operators, the fraction–free kernel method of Section 2 is more efficient than the fractional method whenever the order of the operator or its coefficients are sufficiently large, or whenever the underlying constant field is not ℚ, in which case it seems to be the only practical method. Its performance on third-order operators with smaller coefficients in ℚ(x) seems to indicate that a polyalgorithm is not needed, and that the fraction–free method can be used in general.

In addition, the techniques of applying a fraction–free iteration of a multiple of D (Lemma 2) followed by the preprocessing of Section 3 are applicable to other problems that involve computing a cyclic vector of a matrix constructed from an operator, such as the computation of exterior powers of differential or difference operators [3], or the factorization of completely reducible differential operators [12], although that last problem can also be solved without computing a cyclic vector.

Acknowledgements

We are thankful to Marco Petkov̌sek, Loic Pottier and Emo Welzl for their help and advice on the linear programming aspects of the preprocessing step.

References

[1] Bronstein, M. On radical solutions of linear ordinary differential equations. In Proceedings of the Nijmegen Cathode Workshop (1995), T. Levelt, Ed., pp. 23–24.

[2] Bronstein, M. Σ^{IT}– a strongly-typed embeddable computer algebra library. In Proceedings of DISCO’96 (1996), LNCS 1128, Springer, pp. 22–33.

[3] Bronstein, M., and Petkovšek, M. An introduction to pseudo–linear algebra. Theoretical Computer Science 157 (1996), 3–33.

[4] Char, B., Geddes, K., Gonnet, G., Benton, L., Monagan, M., and Watt, S. Maple V Library Reference Manual. Springer, New York, 1991.

[5] Geddes, K., Czapor, S., and Labahn, G. Algorithms for Computer Algebra. Kluwer Academic Publishers, Boston, 1992.

[6] Hénon, M. A mechanical model for the hitchcock problem. Comptes Rendus de l’Académie des Sciences, Paris, Série I, Mathématiques 321 (1995), 741–745.

[7] Jenks, R., and Sutor, R. Axiom – The Scientific Computation System. Springer, New York, 1992.

[8] McClellan, M. The exact solution of systems of linear equations with polynomial coefficients. Journal of the ACM 20 (1973), 563–588.

[9] Schrijver, A. Theory of Linear and Integer Programming. Wiley, 1986.

[10] Schwartz, J. Fast probabilistic algorithms for verification of polynomial identities. Journal of the ACM 27 (1980), 701–717.

[11] Singer, M. Liouvillian solutions of n^{th} order homogeneous linear differential equations. American Journal of Mathematics 103 (1980), 661–682.

[12] Singer, M. Testing reducibility of linear differential operators: a group theoretic perspective. Applicable Algebra in Engineering, Communication and Computing 7 (1996), 77–104.

[13] Singer, M., and Ulmer, F. Galois groups for second and third order linear differential equations. Journal of Symbolic Computation 16 (1993), 1–36.

[14] Singer, M., and Ulmer, F. Liouvillian and algebraic solutions of second and third order linear differential equations. Journal of Symbolic Computation 16 (1993), 37–73.

[15] Singer, M., and Ulmer, F. Linear differential equations and products of linear forms. In Proceedings of MEGA’96 (to appear).

[16] Ulmer, F., and Weil, J.-A. Note on kovacic’s algorithm. Journal of Symbolic Computation 22 (1996), 179–200.

[17] van Hoeij, M. Factorization of Linear Differential Operators. Computer science dissertation, Nijmegen, 1996. The package diffop is available from

                                                        http://www-math.sci.kun.nl/math/compalg/diffop/.

[18] van Hoeij, M., and Weil, J.-A. An algorithm for computing invariants of differential galois groups. In Proceedings of MEGA’96 (to appear).