EG-eliminations^{∗}

Sergei A. Abramov

Computer Center of the Russian Academy of Science Vavilova 40, Moscow 117967, Russia abramov@ccas.ru, sabramov@cs.msu.su

Abstract We propose an algorithm to put linear recurrent systems in a form which is convenient for using the systems to search for polynomial, power-series, Laurent-series, and other types of solutions of various linear functional systems (differential, difference and q-difference). Some algorithms to search for solutions of functional systems are described. None of the proposed algorithms requires preliminary uncoupling of linear systems.

1 Introduction

Linear recurrences with variable coefficients are of interest for combinatorics and numeric computation. Additionally they give a useful auxiliary tool for constructing solutions of linear functional equations (differential, difference and q-difference) in the form of polynomials, power and Laurent series, rational functions, and so on [4, 8, 7].

Linear recurrent systems are more general and universal instruments. But working with recurrent systems is more complicated than working with scalar recurrences. When we consider a recurrence

p_{l}(n)z_{n+l} + p_{l−1}(n)z_{n+l−1} + ⋯ + p_{t}(n)z_{n+t} = b_{n}

with an unknown sequence z = {z_{n}} then in many situations (e.g., if the coefficients of this recurrence are polynomials) the leading coefficient p_{l}(n) vanishes only for a finite set of values of n. Hence for almost all n we can compute z_{n+l} using z_{n+l−1}, …, z_{n+t}. In a similar manner we can compute z_{n+t} for almost all n using z_{n+l}, …, z_{n+t+1} if the trailing coefficient p_{t}(n) vanishes only for a finite set of values n. In contrast, in the case of a recurrent system S of the form

P_{l}z_{n+l} + P_{l−1}z_{n+l−1} + ⋯ + P_{t}z_{n+t} = b,

where z = (z^{1}, …, z^{m})^{T} is the column of unknown sequences, and P_{l}, …, P_{t} are m × m-matrices over K[n] (K is a ground field), the determinants of P_{l} and/or P_{t} (the leading and trailing matrices of the system) can vanish for all n, even though P_{l} and P_{t} themselves are nonzero.

Let a recurrence appear in the process of constructing the solutions of a given scalar functional equation. Then the leading and trailing coefficients of the recurrence are nonzero polynomials in K[n]. These polynomials have a finite number of integer roots, which give the singularities of the recurrence and the potential degrees of polynomial solutions of the initial scalar equation. In the case of a system, P_{l} and P_{t} can be singular, which prevents us from bounding the degree of the solutions. One faces the necessity to transform such a recurrent system S into an equivalent system S′ with nonsingular leading (or analogously, trailing) matrix.

We propose transformations based on a special process of EG-eliminations in the explicit matrix P = (P_{l}|P_{l−1}| … |P_{t}) of the system S. EG-eliminations bear similarities both to Euclidean algorithm and to Gaussian elimination (E = Euclidean, G = Gaussian), they either allow one to recognize the dependency of the recurrences, or lead to a system S′ of the desired form. It is possible that the transformation of S to S′ additionally results in a finite set of relations for the initial values of each of the solutions of the new system S′.

∗Supported in part by grant from French-Russian Lyapunov institute, Project No 98-03.

For example, using EG-elimination one can prove that the system

( \begin{matrix}n − 1 & 0 \\ −2 & 0 \end{matrix} ) z_{n} + ( \begin{matrix}0 & 0 \\ 0 & n − 2 \end{matrix} ) z_{n−1} + ( \begin{matrix}0 & −1 \\ 0 & 0 \end{matrix} ) z_{n−3} = 0

for z = (z^{1}, z^{2})^{T}, whose leading and trailing matrices are singular, is equivalent to the system

( \begin{matrix}−2 & 0 \\ 0 & n(n − 1) \end{matrix} ) z_{n} + ( \begin{matrix}0 & n − 2 \\ 0 & 0 \end{matrix} ) z_{n−1} + ( \begin{matrix}0 & 0 \\ 0 & −2 \end{matrix} ) z_{n−2} = 0  (1)

with the additional relation 2z_{1}^{1}+z_{0}^{2} = 0. EG-eliminations also show that the original system is equivalent to

( \begin{matrix}0 & 0 \\ (n − 5)(n − 2) & 0 \end{matrix} ) z_{n−1} + ( \begin{matrix}−2 & 0 \\ 0 & 0 \end{matrix} ) z_{n−2} + ( \begin{matrix}0 & n − 4 \\ −2 & 0 \end{matrix} ) z_{n−3} = 0  (2)

with the additional relation 3z_{4}^{1} − z_{1}^{2} = 0. The determinants of the leading matrix of (1) and of the trailing matrix of (2) are nonzero polynomials.

The problem of computing polynomial and series solutions of linear functional system with polynomial coefficients can be solved by the EG-generalization of [4, 8, 7]. Computing the rational function solutions can be done as follows: in the differential case, we can bound the order of the poles of the solutions using the nonsingular leading matrix of the corresponding system of recurrences, which has to be done at all the singularities of the given differential system. In the difference case we can use the methods of [3, 20] instead. The problem is then reduced to computing polynomial solutions of a functional system. In the q-difference case we can combine the differential and difference approaches.

EG-eliminations as a tool of solving functional systems give an alternative for the method of reducing such systems to super-irreducible form [19, 12]. Yet, that method is not optimal, because the superirreducible form contains more information than necessary for the purpose of computing the indicial equations of the systems at their singularities. Furthermore, the reduction algorithms are different for different cases of equations under consideration (differential, difference, q-difference). Moreover, the EG-method is applicable to systems of equations of arbitrary order (rather than first-order for the super-irreducible form) and does not require change of variables in the unknowns.

The traditional computer algebra approach to solving those systems is via the cyclic-vector, or some other similar elimination method [14, 9], that converts the systems to scalar equations (such a procedure is called uncoupling). Gröbner bases technique also can be used to reduce a recurrent system to the uncoupled form [16]. The major, and well-known, problem of this approach is the increase in size of the coefficients of equations which makes those approaches applicable only to systems of very small dimension.

Direct methods (i.e., methods which work without uncoupling of systems) for solving systems of linear functional equations are needed to design effective algorithms to compute invariants of Galois groups of linear ordinary differential equations [22], which in turn would allow the efficient computation of closed-form (Liouvillian) solutions of those equations [21]. Other applications would include efficient factoring algorithms for completely reducible linear ordinary differential operators [23]. In addition, such methods are helpful in designing effective Gröbner bases algorithms in multivariate Ore rings [15, 17]. Those bases are important in generalizing Zeilberger’s definite summation and integration algorithm.

Recurrent systems appear also in combinatorics and in numerical analysis, but their leading and trailing matrices can be singular, which makes them useless in computing the elements of the solution sequences. Applying the EG-eliminations allows us to desingularize the leading and trailing matrices, thereby computing the sequences.

2 Linear recurrent systems

2.1 Preliminaries

Let K be a field of characteristic zero. For a sequence

c = {c_{n}}_{n≥k},  (3)

k ∈ Z, with elements in K we write ν(c) = k (note that it is possible that c_{k} = 0). If c has the form (3) and f is a function Z → K then the product fc is equal to {f(n)c_{n}}_{n≥k}. The sum of two sequences c and c′ is defined elementwise and ν(c + c′) = max{ν(c), ν(c′)}. The sequences d = Ec and d′ = E^{−1}c are defined so that ν(d) = ν(c) − 1, ν(d′) = ν(c) + 1 with d_{n} = c_{n+1} for n = ν(d), ν(d) + 1, … and d′_{n} = c_{n−1} for n = ν(d′), ν(d′) + 1, … We say that c of the form (3) satisfies the equation (the recurrence)

R(z) = 0,  (4)

where R is a recurrent operator of the form

p_{l}(n)E^{l} + p_{l−1}(n)E^{l−1} + ⋯ + p_{t}(n)E^{t}  (5)

(l, t ∈ Z; p_{t}(n), p_{l}(n) ≠ 0) if applying R to c gives the sequence {d_{n}}_{n≥k−l} with zero elements. For R of the form (5) we write

ord^{∗}(R) = l.  (6)

From here on we will assume the coefficients of the recurrent operators to be belonging to a ring I of functions Z → K such that

I1. The shift operator E is an automorphism of I.

I2. If f(n) ∈ I and f(n) is not equal to zero identically then the equation f(n) = 0 has only finite set (possibly empty) of integer solutions.

Example 1

a) The ring K[n] of polynomial functions satisfies I1, I2.

b) Let K = K_{0}(q), where q is transcendental over K_{0}, then the ring K[q^{n}], considered as a ring of functions Z → K, satisfies I1, I2.

In the applications we will mainly consider the case I = K[n]. In the examples of the nearest sections we assume I to be K[n] and K to be the rational number field Q. An example with I = K[q^{n}] will be considered in Section 3.6.

The ring of all recurrent operators of the form (5) over I will be denoted by E_{I}.

Let b = {b_{n}}_{n≥s} and let R be an operator of the form (5). A solution of recurrence

R(z) = b  (7)

is any sequence of the form (3) such that (R(c))_{n} = b_{n} for n ≥ max{k − t, s}.

We will consider also sequences of the form {c_{n}}_{n∈Z}, ν(c) = −∞. We denote the set of such sequences by A. For any R of the form (5) and any sequence c ∈ A we have R(c) ∈ A. We will be interested in the solutions of (7) in the case b ∈ A and especially when b_{n} = 0 for all n < 0 (the set of all such sequences is denoted by A_{0}). For example, we will consider the problem of looking for solutions in A_{0} (or, equivalently, A_{0}-solutions) of a recurrence (7) with b ∈ A_{0}. Observe that the last problem can be reduced to the search for solutions c̃ such that ν(c̃) = t − l and c̃ = 0 for n = −1, −2, …, t − l. When such a sequence c̃ is found one can extend it for n = t − l − 1, t − l − 2, … by zero elements.

Now we concentrate on recurrent systems of the form

P_{l}z_{n+l} + P_{l−1}z_{n+l−1} + ⋯ + P_{t}z_{n+t} = b  (8)

where z = (z^{1}, …, z^{m})^{T} is the column of unknown sequences;

b = (b^{1}, …, b^{v})^{T}  (9)

(the right-hand side) is the column of known sequences; v is a nonnegative integer; P_{l}, …, P_{t} are v × mmatrices over I, with non-zero P_{l} and P_{t} (the leading and trailing matrices of the system). For a column b of the form (9) we set ν(b) = max{ν(b^{1}), …, ν(b^{v})}. If ν(b) = −∞ then one can consider b_{n} for any n, but if ν(b) > −∞ then only for n ≥ ν(b). The substitution of a column c into the left-hand side of (8) gives a column d such that ν(d) = ν(c) − t. If d_{n} = b_{n} for all n ≥ max{ν(d), ν(b)} (in the case max{ν(d), ν(b)} > −∞) or for all n (in the case ν(d), ν(b) = −∞), then column (8) is a solution of recurrent system (8). In the case when b ∈ A_{0}^{v} and one is looking for a solution in A_{0}^{m} (as in the scalar case we will call such a solution A_{0}-solution) the problem can be reduced to the search for solutions c̃ such that

ν(c̃) = t − l  (10)

and

t − l ≤ n ≤ −1 ⇒ c̃_{n} = 0.  (11)

We call the matrix

P = (P_{l}|P_{l−1}| … |P_{t})  (12)

the explicit matrix of the system S. The matrix P_{l} is the leading part of the matrix P, and the matrix P_{t} is its trailing part. One can consider the leading and trailing parts of any row of P.

We will say that b ∈ A^{v} is of finite degree if there exists an integer w ≥ −1 such that b_{n} is the zero column for all n > w. The minimal such w will be denoted by deg b.

Lemma 1 Let v = m in system (8) and the right-hand side of (8) be of finite degree. Let p(n) = det P_{t}(n) be a nonzero element of I. Let n_{0} be the largest integer root of p(n) = 0 if such roots exist and n_{0} = −1 otherwise. Let c ∈ A^{m} be a solution of (8) of finite degree. Then

deg c ≤ max{n_{0} + t, deg b + t, −1}.  (13)

Proof: Let s = max{n_{0} +t, −1} and s′ ≥ s. The columns c_{s′}, c_{s′+1}, … allow one to compute the columns c_{s+1}, c_{s+2}, …, c_{s′−1} (the matrix P_{t}(n) is an invertible matrix with elements belonging to K for all n > n_{0}, it lets one use (8) for computing columns c_{s′−1}, c_{s′−2}, …, c_{s+1}). If all the columns c_{n} are zero for all large enough n then c_{n} will be zero for all n > max{s, deg b + t}. Therefore deg c ≤ max{s, deg b + t} which gives (13). □

Using this bound one can construct all A_{0}-solution of finite degree of system (8). It can be done by undetermined coefficient method or by the method from [4].

If the leading part P_{l} of P is nonsingular (i.e., its determinant is a nonzero element of I) and b ∈ A_{0}^{m} then we can describe all A_{0}-solutions of system (8). In the general case the following lemma is useful.

Lemma 2 Let v = m in system (8) and p(n) = det P_{l}(n) be a nonzero element of I. Let n_{1} be the largest integer root of p(n) = 0 if such roots exist and n_{1} = −1 otherwise. Let c ∈ A^{m} be a solution of (8) and let c_{n} be known for n ≤ n_{1} + l. Then c_{n} are uniquely defined for all

n > n_{1} + l.  (14)

Proof: The matrix P_{l}(n) is an invertible matrix with elements belonging to K for all n > n_{1}, it lets one use (8) for computing columns c_{n_{1}+l+1}, c_{n_{1}+l+2}, . . .. □

Therefore it is sufficient to know the set c_{[0:s]} = (c_{0}, …, c_{s}), s = n_{1} + l, to find all other columns of an A_{0}-solution c of system (8) using recurrent system (8). Such c_{[0:s]} can be considered as an m × (s + 1)- matrix over K. Thanks to (10) and (14) we can find all c_{[0:s]} as follows. Taking into account (10) substitute step by step n = −l, −l + 1, …, n_{1} + l into (8). It gives us a system of linear algebraic equations for the matrix c_{[0:s]}. The method from [4] can be used also to solve this problem. Both methods allow one to obtain a basis of the affine space of all the matrices c_{[0:s]}.

But unfortunately in the case v = m the matrices P_{t} and P_{l} are singular very often. The situation is more difficult in the case v ≠ m. In the next section we will describe some transformations which allow one to get a system with non-singular square P_{l} or, resp., with non-singular square P_{t}.

2.2 Constraints, equivalent transformations and dependence of equations

We start with definitions.

We call a constraint for the sequence (3) any relation of the form

α_{M}c_{M} + α_{M−1}c_{M−1} + ⋯ + α_{N}c_{N} = σ,

where α_{M}, …, α_{N}, σ ∈ K, M ≥ N ≥ k. We stress that M and N are not variables but are concrete integer numbers. We will also consider constraints for elements of different sequences c^{1}, …, c^{m}:

αc_{M}^{1} + ⋯ + βc_{N}^{1} + ⋯ + γc_{M}^{m} + ⋯ + δc_{N}^{m} = σ.

Any recurrent system is equal to an infinite set of constraints.

System (8) can be rewritten in the form

( \begin{matrix}R_{11} & … & R_{1m} \\ ⋮ & & ⋮ \\ R_{v1} & … & R_{vm} \end{matrix} ) ( z^{1} ; ⋮ ; z^{m} ) = ( b^{1} ; ⋮ ; b^{v} ) ,  (15)

R_{ij} ∈ E_{I}, i = 1, …, v, j = 1, …, m. We will call the equations of this system dependent if there exist S_{1}, …, S_{v} ∈ E_{I} such that

S_{1} ∘ R_{1j} + ⋯ + S_{v} ∘ R_{vj} = 0  (16)

for all j = 1, …, m. If the last equalities take place but the sequence S_{1}b^{1} + ⋯ + S_{v}b^{v} is not zero then the original recurrent system is not compatible.

Consider the following particular case of the recurrent system: a first order system in the canonical form:

p_{1}(n)z_{n+1}^{1} = a_{11}(n)z_{n}^{1} + ⋯ + a_{1m}(n)z_{n}^{m} + b_{n}^{1}
· · · · · · · · ·  (17)
p_{v}(n)z_{n+1}^{v} = a_{v1}(n)z_{n}^{1} + ⋯ + a_{vm}(n)z_{n}^{m} + b_{n}^{v},

v ≤ m; p_{1}(n), …, p_{v}(n) ≠ 0.

Lemma 3 Equations (17) are independent.

Proof: System (17) can be rewritten in the form (15) with the matrix

( \begin{matrix}p_{1}E − a_{11} & −a_{12} & … & −a_{1v} & … & −a_{1m} \\ −a_{21} & p_{2}E − a_{22} & … & −a_{2v} & … & −a_{2m} \\ ⋮ & ⋮ & & ⋮ & & ⋮ \\ −a_{v1} & −a_{v2} & … & p_{v}E − a_{vv} & … & −a_{vm} \end{matrix} ) .

If its rows are dependent over E_{I} then there exist S_{1}, …, S_{v} ∈ E_{I} such that

− ∑_{i=1}^{v} S_{i} ∘ a_{ij} + S_{j} ∘ p_{j} ∘ E = 0,  (18)

j = 1, …, v. Choose j such that

ord^{∗}(S_{j}) = max{ord^{∗}(S_{1}), …, ord^{∗}(S_{v})}  (19)

(see (6)). But equality (18) does not hold for such j because ord^{∗}(S_{j} ∘ p_{j} ∘ E) > ord^{∗}(S_{i} ∘ a_{ij}) for i = 1, …, v. □

From here on we consider only the case when ν(b) = −∞ and the solutions c which we are looking for are also such that ν(c) = −∞. By equivalent transformations of systems we mean transformations which preserve these solutions.

Let’s concentrate on two types of equivalent transformations.

1. Applying E^{s} to the i-th equation of the system, s ∈ Z, 1 ≤ i ≤ v, i.e., going from (15) to

( \begin{matrix}R_{11} & … & R_{1m} \\ ⋮ & & ⋮ \\ R_{i−1,1} & … & R_{i−1,m} \\ E^{s} ∘ R_{i1} & … & E^{s} ∘ R_{im} \\ R_{i+1,1} & … & R_{i+1,m} \\ ⋮ & & ⋮ \\ R_{v1} & … & R_{vm} \end{matrix} ) ( z^{1} ; ⋮ ; ⋮ ; ⋮ ; ⋮ ; z^{m} ) = ( b^{1} ; ⋮ ; b^{i−1} ; E^{s}b^{i} ; b^{i+1} ; ⋮ ; b^{m} ) .

This transformation will be called the transformation of type T1 or, simpler, the transformation T1. The equivalence of this transformation can be proven by considering the sets of constraints which are equivalent for both systems.

2. Now let g(n) be a nonzero element of I. The multiplication of the i-th equation of the system by g gives the system

( \begin{matrix}R_{11} & … & R_{1m} \\ ⋮ & & ⋮ \\ R_{i−1,1} & … & R_{i−1,m} \\ gR_{i1} & … & gR_{im} \\ R_{i+1,1} & … & R_{i+1,m} \\ ⋮ & & ⋮ \\ R_{v1} & … & R_{vm} \end{matrix} ) ( z^{1} ; ⋮ ; ⋮ ; ⋮ ; ⋮ ; z^{m} ) = ( b^{1} ; ⋮ ; b^{i−1} ; gb^{i} ; b^{i+1} ; ⋮ ; b^{m} ) .

The equation

gR_{i1}z^{1} + ⋯ + gR_{im}z^{m} = gb^{i}

is equivalent to the original equation

R_{i1}z^{1} + ⋯ + R_{im}z^{m} = b^{i}  (20)

if the equation g(n) = 0 has no integer roots. Otherwise to get equivalence we have to add to the new recurrent system a finite set of constraints which are obtained by substitution of all the roots of g(n) = 0 for n in (20). If we are interested only in A_{0}-solutions then we can consider only such roots g(n) = 0 that are ≥ −l (see (8)).

In a similar way we can replace the i-th equation of the system by the sum of the product of g(n) by the i-th equation and the product of h(n) by the j-th equation, 1 ≤ j ≤ v, g(n), h(n) ∈ I, adding again to the recurrent system the finite set of constraints which are obtained by substitution of all those integer roots of g(n) which are ≥ −l, for n in (20) (transformation T2).

Example 2 Consider the homogeneous recurrence R(z) = 0, where R = (n + 1)E. Construct R_{1} = E ∘ R = (n + 2)E^{2} (transformation T1). As mentioned above the sets of A_{0}-solutions of the recurrences R(z) = 0 and R_{1}(z) = 0 are equal. Now construct R_{2} = (n + 1)R_{1} (transformation T2). It is easy to see that l = 2 for R_{1}. We have −1 ≥ −2, and we write down the constraint which is the result of the substitution of n = −1 into (n + 2)z_{n+2} = 0. It results in z_{1} = 0. Therefore the set of A_{0}-solutions of the recurrence

(n + 1)z_{n+1} = 0  (21)

coincides with the set of A_{0}-solutions of the recurrence

(n + 1)(n + 2)z_{n+2} = 0,  (22)

supplemented by the constraint z_{1} = 0. This is in accordance with the fact that A_{0}-solutions of recurrence (21) are described by the formula

z_{n} = \{ \begin{matrix}C, & if n = 0, \\ 0, & otherwise, \end{matrix}

where C is an arbitrary constant. In turn A_{0}-solutions of recurrence (22) are described by the formula

z_{n} = \{ \begin{matrix}C_{0}, & if n = 0, \\ C_{1}, & if n = 1, \\ 0, & otherwise, \end{matrix}

where C_{0}, C_{1} are arbitrary constants.

Now we consider the corresponding transformations of the explicit matrix P of the original recurrent system. Let the i-th row of the matrix P_{t} be zero. Then applying E^{−1} to the i-th equation is equivalent to the transformation of the i-th row of the matrix P

(p_{i1}(n), …, p_{i,(l−t)m}(n), 0, …, 0)  (23)

to

(0, …, 0, p_{i1}(n − 1), …, p_{i,(l−t)m}(n − 1)).  (24)

We will call this transformation (which is a particular case of T1) transformation t1^{−}. In turn transformation t1^{+} can be applied in the case where the i-th row of the matrix P_{l} is zero. Then applying E to the i-th equation is equivalent to the transformation of the i-th row of P

(0, …, 0, p_{i,m+1}(n), …, p_{i,(l−t+1)m}(n)),

to

(p_{i,m+1}(n + 1), …, p_{i,(l−t+1)m}(n + 1), 0, …, 0).

Therefore transformations t1^{−}, t1^{+} correspond to applying E^{−1} or E in the case where this neither decreases t nor increases l.

Additionally we will consider transformation t2. This replaces the i-th row of P by the sum of the product of g(n) by the i-th row and the product of h(n) by the j-th row, 1 ≤ j ≤ v, g(n), h(n) ∈ I, g(n) ≠ 0, and writing down the constraint

p_{i1}(n_{0})z_{n_{0}+l}^{1} + ⋯ + p_{im}(n_{0})z_{n_{0}+l}^{m}+
· · · · · · · · ·  (25)
+p_{i,(l−t)m+1}(n_{0})z_{n_{0}+t}^{1} + ⋯ + p_{i,(l−t+1)m}(n_{0})z_{n_{0}+t}^{m} = 0

for any integer roots n_{0} ≥ −l of g(n) = 0. We call t1^{−}, t1^{+}, t2 the elementary transformations. With any elementary transformation of the explicit matrix P a transformation of the operator matrix R (which is an m × v-matrix over E_{I}) can be associated. Such a new transformation is the left multiplication of some v × v-matrix over E_{I} by R. The multiplications by, resp., matrices

( \begin{matrix}1 \\ & . & . & . \\ & & & & 1 \\ & & & & & E \\ & & & & & & 1 \\ & & & & & & & . & . & . \\ & & & & & & & & & & 1 \end{matrix} ) , ( \begin{matrix}1 \\ & . & . & . \\ & & & & 1 \\ & & & & & E^{−1} \\ & & & & & & 1 \\ & & & & & & & . & . & . \\ & & & & & & & & & & 1 \end{matrix} ) , ( \begin{matrix}1 \\ & . & . & . \\ & & & & 1 \\ & & & & & g(n) & & h(n) \\ & & & & & & 1 \\ . & . & . \\ & & & & & & & & 1 \end{matrix} ) .

are associated with transformations t1^{+}, t1^{−}, t2. The entries E and E^{−1} of the first and the second matrices have indices (i, i), the entries g(n) and h(n) of the third matrix have indices (i, i) and (i, j) respectively. We will call the matrices of such kind also elementary.

We consider additionally the elementary transformation t3, which interchanges the i-th and the j-th rows of the explicit matrix P (though this transformation can be reduced to t2, it is more convenient to use itself). The corresponding elementary matrix is of the form

( \begin{matrix}1 \\ & . & . & . \\ & & & & 1 \\ & & & & & 0 & & & & & & 1 \\ & & & & & & 1 \\ & & & & & & & . & . & . \\ & & & & & & & & & & 1 \\ & & & & & 1 & & & & & & 0 \\ & & & & & & & & & & & & 1 \\ & & & & & & & & & & & & & . & . & . \\ & & & & & & & & & & & & & & & & 1 \end{matrix} ) .

The non-diagonal unity entries have indices (i, j) and (j, i).

The right-hand side of a recurrent system under an elementary transformation is multiplied by the corresponding elementary matrix.

Now we consider the problem of reducing the explicit matrix P of the given recurrent system to a matrix P′ with the trailing matrix in the trapezoidal form during which if the trailing part of a row of P′ is zero then the row will be completely zero (this is the t-problem of reducing the explicit matrix). It will be done by elementary transformations t1^{−}, t2, t3, i.e., by t1^{−} and by some “delicate” Gaussian eliminations. We also consider the similar problem of reducing the explicit matrix P to a matrix P′′ with the leading matrix in the trapezoidal form during which if the leading part of a row of P′′ is zero then the row will be completely zero (the l-problem of reducing the explicit matrix). It will be done by elementary transformations t1^{+}, t2, t3, i.e., by t1^{+} and again by some “delicate” Gaussian eliminations. It is obvious that if zero rows appeared in the process of solving the t- or l-problem then the equations of the original recurrent system are dependent; the corresponding row of the product of the elementary matrices involved gives S_{1}, …, S_{v} from (16).

Observe that adding t1^{−} to the “usual” Gaussian eliminations often allows one to solve the t-problem of reducing a given matrix.

Example 3 Consider the following recurrent system

( \begin{matrix}(n + 1)E − 1 & −1 \\ −1 & (n + 1)E − 1 \end{matrix} ) ( \begin{matrix}z^{1} \\ z^{2} \end{matrix} ) = ( \begin{matrix}0 \\ 0 \end{matrix} )  (26)

or

\begin{matrix}(n + 1)z_{n+1}^{1} − z_{n}^{1} − z_{n}^{2} = 0 \\ −z_{n}^{1} + (n + 1)z_{n+1}^{2} − z_{n}^{2} = 0. \end{matrix}  (27)

The explicit matrix of this system is of the form

( \begin{matrix}n + 1 & 0 & −1 & −1 \\ 0 & n + 1 & −1 & −1 \end{matrix} ) ,

where the first two columns form the leading part and the last two ones form the trailing part of the explicit matrix. Transform the explicit matrix:

( \begin{matrix}n + 1 & 0 & −1 & −1 \\ 0 & n + 1 & −1 & −1 \end{matrix} ) ⟶^{1} ( \begin{matrix}n + 1 & 0 & −1 & −1 \\ −n − 1 & n + 1 & 0 & 0 \end{matrix} ) ⟶^{2} ( \begin{matrix}n + 1 & 0 & −1 & −1 \\ 0 & 0 & −n & n \end{matrix} ) ⟶^{3}

( \begin{matrix}n + 1 & 0 & −1 & −1 \\ −n(n + 1) & 0 & 0 & 2n \end{matrix} ) .

Here ⟶^{1} and ⟶^{3} are transformations t2; ⟶^{2} is transformation t1^{−}. As in ⟶^{1} , ⟶^{3} the changed row is multiplied by 1, we write down no additional constraints. The determinant of the trailing part of the last matrix is equal to −2n, the equation −2n = 0 has the only nonnegative integer root n_{0} = 0. By (13) this means that if there exists a finite degree A_{0}-solution of the original system then its degree is 0. It is easy to see that such a solution really exists:

z_{n}^{1} = \{ \begin{matrix}C, & if n = 0, \\ 0, & otherwise, \end{matrix}

z^{2} = −z^{1},

where C is an arbitrary constant. We can construct the system whose explicit matrix is equal to the matrix which we got as the result of transformation ⟶^{3} :

\begin{matrix}(n + 1)z_{n+1}^{1} − z_{n}^{1} − z_{n}^{2} = 0 \\ −n(n + 1)z_{n+1}^{1} − 2nz_{n}^{2} = 0, \end{matrix}

or

( \begin{matrix}(n + 1)E − 1 & −1 \\ −n(n + 1)E & −2n \end{matrix} ) ( \begin{matrix}z^{1} \\ z^{2} \end{matrix} ) = ( \begin{matrix}0 \\ 0 \end{matrix} ) .  (28)

But sometimes this approach does not lead to the desired result.

Example 4 Consider the following system of two recurrences with one unknown sequence z:

( \begin{matrix}E + 1 \\ E + 2 \end{matrix} ) ( z) = ( \begin{matrix}0 \\ 0 \end{matrix} ) .

The explicit matrix of this system is

( \begin{matrix}1 & 1 \\ 1 & 2 \end{matrix} ) ,  (29)

its leading and, resp., trailing parts are

( \begin{matrix}1 \\ 1 \end{matrix} ) and ( \begin{matrix}1 \\ 2 \end{matrix} ) .

Gaussian eliminations in the second row of (29) and shifts of this row give

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

( \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ) ⟶ ( \begin{matrix}1 & 1 \\ −1 & 0 \end{matrix} ) ⟶ ⋯

There are nonzero elements in the second row of all these matrices although the two equations are obviously dependent. One has to use some “delicate” Gaussian eliminations.

2.3 EG-eliminations

The Euclidean algorithm suggests a way out. In this algorithm the degree or the order of the divisor (i.e., the polynomial or the operator by which one eliminates something) does not exceed the degree or the order of the dividend (i.e., the polynomial or the operator in which one eliminates something). We will use a matrix analog of the rule, introducing for this purpose the notion of the width of a row and the notion of EG-elimination (E=Euclidean, G=Gaussian). The width of a row is an analog of the degree of a polynomial and of the order of an operator. Strictly speaking we will use two types of the width, namely l-width and t-width.

Let as usual the original system have the form (8). If the i-th row of the explicit matrix P is zero then its l-width and t-width are both equal to 0. Let in (8) the i-th rows of all the matrices P_{l}, P_{l−1}, …, P_{s+1} be zero, while the i-th row of P_{s} is nonzero. Then s − t + 1 is the t-width of the i-th row of P. Similarly, let the i-th rows of P_{t}, P_{t+1}, …, P_{s−1} be zero, while the i-th row of P_{s} is nonzero. Then l − s + 1 is the l-width of the i-th row of P.

We will consider two versions of EG-eliminations:

a) EG_{t}-eliminations which are based on the notion of t-width and can be used for reducing P to a matrix with the trailing matrix in the trapezoidal form,

b) EG_{l}-eliminations which are based on the notion of l-width and can be used for reducing P to a matrix with the leading matrix in the trapezoidal form. We will detail EG_{t}-eliminations. As for EG_{l}-eliminations, they can be described in a similar way.

EG_{t}-elimination of p_{ik} in P by the j-th row is possible in the case where j < i and p_{jk} ≠ 0. If p_{ik} = 0 then P will not be changed. If p_{ik} ≠ 0 then:

1. Compare the t-width of the j-th row with the t-width of the i-th row. If the first value is bigger than the second one then interchange the i-th and the j-th rows.

2. Now the t-width of the j-th row does not exceed the t-width of the i-th row and p_{jk} ≠ 0. Eliminate p_{ik} by the j-th row: find g, h ∈ I such that gp_{ik} + hq_{jk} = 0 and replace the i-th row by the sum of the product of g by the i-th row and the product of h by the j-th row. Write down the constraint (25) for every integer root n_{0} of g(n) = 0 such that n_{0} ≥ −l. It is easy to see that any EG_{t}-elimination is a chain of elementary transformations.

Reducing P to a matrix with the trailing part in the trapezoidal form can be done in v steps (as above v is the number of the system equations or, equivalently, the number of the rows of P). The w-th step results in the integer numbers w_{1}, w_{2} such that w_{1} + w_{2} = w (initially w_{1} = w_{2} = 0) during which the matrix P will be transformed by elementary transformations in such a way that

• the first w_{1} rows of the trailing part form a trapezoidal matrix;

• the last w_{2} rows of matrix P are zero.

It is clear that if all v steps are performed then the matrix P will be reduced to the desirable form.

Assume that the steps with numbers ≤ w have been performed (at the beginning w = w_{1} = w_{2} = 0). Then the algorithm for step no. w + 1 is as follows:

1) Use EG_{t}-eliminations by the rows with numbers 1, …, w_{1} to eliminate the first w_{1} entries of the trailing part of the (w_{1} + 1)-st row (interchanging the rows when the t-width of the row by which elimination is done exceeds that of row w_{1} + 1).

2) If there is a nonzero entry in the trailing part of the (w_{1} + 1)-th row of P then change enumeration of the unknowns z^{w_{1}+1}, …, z^{m} in such a way that this entry (or one of such entries) will have position no. w_{1} + 1 in the trailing part. All entries in positions 1, …, w_{1} are zero. The step is finished by increasing w_{1} by 1.

3) If the (w_{1} + 1)-st row of P is completely zero then interchange this row and the row with number v − w_{2}. The step is finished by increasing w_{2} by 1.

4) If the (w_{1} + 1)-st row of P is non-zero but at the same time the trailing part of this row has only zero entries then apply t1^{−} to the (w_{1} + 1)-st row of P and go to 1).

The step no. w + 1 is now completely described.

Theorem 1 The above procedure implementing step no. w + 1 always terminates.

Proof: The algorithm terminates when it reaches cases 2) or 3). In case 4) the (w_{1} + 1)-st row of P is non-zero but with zero trailing part. Consider the integer vector (s_{1}, …, s_{w_{1}+1}), where s_{j} is equal to the t-width of the j-th row before EG_{t}-eliminations in the (w_{1} + 1)-st row. Let this vector be transformed into (s′_{1}, …, s′_{w_{1}+1}) as the result of EG_{t}-eliminations in the (w_{1} + 1)st row and of applying t1^{−} to this row. Then it is easy to see that s′_{1} ≤ s_{1}, …, s′_{w_{1}} ≤ s_{w_{1}} and if s′_{1} = s_{1}, …, s′_{w_{1}} = s_{w_{1}} then s′_{w_{1}+1} < s_{w_{1}+1}. Therefore (s′_{1}, …, s′_{w_{1}+1}) <_{lex} (s_{1}, …, s_{w_{1}+1}), where <_{lex} is the natural lexicographic order. The components of the vectors under consideration are nonnegative integers and the theorem follows. □

EG_{t}-eliminations allow one to solve the problem formulated in Example 4:

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

( \begin{matrix}0 & −1 \\ 1 & 0 \end{matrix} ) ⟶ ( \begin{matrix}0 & −1 \\ 0 & 1 \end{matrix} ) ⟶ ( \begin{matrix}0 & −1 \\ 0 & 0 \end{matrix} ) .

We get the matrix with a zero row and a zero column. By deleting them we obtain the explicit matrix of the system which consists of one recurrence −z_{n} = 0. The only solution of this recurrence is the zero sequence. The zero row appearance means that the rows of matrix (29) are dependent over E_{I}. In our case this fact is trivial. But the given approach allows one to recognize dependency (or independency) of the rows of the explicit matrix of any system. The additional computation of the product of the corresponding elementary matrices allows (as was mentioned) to find a relationship (16) for the original operator matrix. It allows also to construct the right-hand sides of the new system in the non-homogeneous case. In practice, of course, we can transform the right-hand side together with the row of the explicit matrix of the system.

Example 5 The recurrent system

n(n + 1)z_{n+1}^{1} + z_{n}^{1} + z_{n}^{2} = 0
n(n + 1)z_{n+1}^{1} − (n + 1)z_{n}^{2} = −δ_{n2}
z_{n+1}^{1} + (n + 3)z_{n+1}^{2} = δ_{n1},

has the explicit matrix

( \begin{matrix}n(n + 1) & 0 & 1 & 1 \\ n(n + 1) & 0 & 0 & −n − 1 \\ 1 & n + 3 & 0 & 0 \end{matrix} ) .

We extend this matrix by the right-hand side column and perform EG_{l}-eliminations:

( \begin{matrix}n(n + 1) & 0 & 1 & 1 & 0 \\ n(n + 1) & 0 & 0 & −n − 1 & −δ_{n2} \\ 1 & n + 3 & 0 & 0 & δ_{n1} \end{matrix} ) ⟶ ( \begin{matrix}1 & n + 3 & 0 & 0 & δ_{n1} \\ n(n + 1) & 0 & 1 & 1 & 0 \\ n(n + 1) & 0 & 0 & −n − 1 & −δ_{n2} \end{matrix} ) ⟶

( \begin{matrix}1 & n + 3 & 0 & 0 & δ_{n1} \\ 0 & −n(n + 1)(n + 3) & 1 & 1 & −2δ_{n1} \\ 0 & −n(n + 1)(n + 3) & 0 & −n − 1 & −δ_{n2} − 2δ_{n1} \end{matrix} ) ⟶

( \begin{matrix}1 & n + 3 & 0 & 0 & δ_{n1} \\ 0 & −n(n + 1)(n + 3) & 1 & 1 & −2δ_{n1} \\ 0 & 0 & −1 & −n − 2 & −δ_{n2} \end{matrix} ) ⟶

( \begin{matrix}1 & n + 3 & 0 & 0 & δ_{n1} \\ 0 & −n(n + 1)(n + 3) & 1 & 1 & −2δ_{n1} \\ −1 & −n − 3 & 0 & 0 & −δ_{n1} \end{matrix} ) ⟶ ( \begin{matrix}1 & n + 3 & 0 & 0 & δ_{n1} \\ 0 & −n(n + 1)(n + 3) & 1 & 1 & −2δ_{n1} \\ 0 & 0 & 0 & 0 & 0 \end{matrix} ) .

We have reduced the original system to an equivalent one:

\begin{matrix}z_{n+1}^{1} + (n + 3)z_{n+1}^{2} = δ_{n1} \\ −n(n + 1)(n + 3)z_{n+1}^{2} + z_{n}^{1} + z_{n}^{2} = −2δ_{n1}, \end{matrix}

which can be rewritten in the form

( \begin{matrix}1 & n + 3 \\ 0 & −n(n + 1)(n + 3) \end{matrix} ) ( \begin{matrix}z_{n+1}^{1} \\ z_{n+1}^{2} \end{matrix} ) + ( \begin{matrix}0 & 0 \\ 1 & 1 \end{matrix} ) ( \begin{matrix}z_{n}^{1} \\ z_{n}^{2} \end{matrix} ) = ( \begin{matrix}δ_{n1} \\ −2δ_{n1} \end{matrix} ) .  (30)

The largest nonnegative integer root of the determinant of the leading matrix is 0; in correspondence with (14) it is sufficient to find

z_{0}^{1}, z_{0}^{2}, z_{1}^{1}, z_{1}^{2},  (31)

and all other elements of A_{0}-solutions will be found by recurrence (30). The values (31) together with z_{−1}^{1} = z_{−1}^{2} = 0 have to satisfy (30) when n = −1, 0; this gives the linear algebraic equations

z_{0}^{1} + 2z_{0}^{2} = 0, z_{1}^{1} + 3z_{1}^{2} = 0, z_{0}^{1} + z_{0}^{2} = 0.  (32)

We get z_{0}^{1} = z_{0}^{2} = 0, z_{1}^{1} = −3C, z_{1}^{2} = C where C is an arbitrary constant. Taking into account that for n > 1 the right-hand side of (30) is zero, we consider the case n = 1 separately: z_{2}^{1} = C, z_{2}^{2} = \frac{1−C}{4} . Due to

( \begin{matrix}1 & n + 3 \\ 0 & −n(n + 1)(n + 3) \end{matrix} )^{−1} = ( \begin{matrix}1 & \frac{1}{n(n+1)} \\ 0 & \frac{−1}{n(n+1)(n+3)} \end{matrix} ) ,

we have for n = 2, 3, …

( \begin{matrix}z_{n+1}^{1} \\ z_{n+1}^{2} \end{matrix} ) = ( \begin{matrix}\frac{−1}{n(n+1)} & \frac{−1}{n(n+1)} \\ \frac{1}{n(n+1)(n+3)} & \frac{1}{n(n+1)(n+3)} \end{matrix} ) ( \begin{matrix}z_{n}^{1} \\ z_{n}^{2} \end{matrix} ) .

It is possible that the leading part after EG_{l}-eliminations will have the non-triangular trapezoidal form. Let the part have the size k × v, k < v. The unknowns z^{k+1}, …, z^{v}, in contrast to the case of linear algebraic systems, cannot be taken arbitrary in the general situation. The obtained linear algebraic equalities for the initial values (like equalities (32)) can keep the initial values of the sequences z^{k+1}, …, z^{v} from being arbitrary.

Now the last remark of the section. Let b be in A_{0}^{v} in the original system. After the performance of EG_{t}-eliminations we can get the new right-hand side b̂ = (b̂^{1}, …, b̂^{v})^{T}, b̂ ∉ A_{0}^{v}. But there exists a negative integer k such that we have

n < k ⇒ b̂_{n}^{i} = 0

for all i = 1, …, v. We can apply E^{k} to any equation of the new system. It gives a system whose right-hand side is E^{k}b̂ = (E^{k}b̂^{1}, …, E^{k}b̂^{v})^{T} ∈ A_{0}^{v}.

2.4 Three computing remarks

2.4.1 Ccancelling common factors from rows

When we perform EG-eliminations it is desirable to cancel common factors from the rows of the explicit matrices P. But linear recurrences can lose solutions due to such cancelling. The simplest example is the recurrence nz_{n} = 0 of order zero. Its A_{0}-solutions have the form

z_{n} = \{ \begin{matrix}0, & if n ≠ 0 \\ C, & if n = 0, \end{matrix}

where C is an arbitrary constant. After cancelling, the recurrence z_{n} = 0 has only the zero solution. From here on we will follow the principle that if the factor p(n) ∈ I with an integer root n_{0}, n_{0} ≥ −l, was cancelled from a recurrence and if we use the obtained recurrence with n = n_{0} then we have to assume all the recurrence coefficients to be zero.

Working with the explicit matrix P we will store a finite set (possible empty) A of pairwise different integers ≥ −l together with each of its rows. Such a row will be called a row with a removed factor. Let us have two rows with removed factors:

A, (a_{1}(n), a_{2}(n), …),  (33)

where a_{1}(n) ≠ 0, and

Ã, (ã_{1}(n), ã_{2}(n), …).  (34)

Let in the row

(ã_{1}(n), ã_{2}(n), …)  (35)

an entry (e.g., ã_{1}(n)) be eliminated by the help of the row

(a_{1}(n), a_{2}(n), …).  (36)

To do this the rows were multiplied by factors f(n), f̃(n) ∈ I such that f(n)a_{1}(n) = −f̃(n)ã_{1}(n), and added together. Let the result be (0, ˜ ã_{2}(n), ˜ ã_{3}(n), …). This means that if elimination is performed in the row (34) with the help of (33) then we will get

A ∪ Ã, (0, ˜ ã_{2}(n), ˜ ã_{3}(n), …),

and new constraints can appear here. These constraints are generated by the recurrence corresponding to the row (35) and by values n belonging to the set (A \ Ã) ∪ F, where F is the set of integer roots ≥ −l of f(n) = 0.

Applying t1^{−}, t1^{+} to (33) changes the row (a_{1}(n), a_{2}(n), …) in the usual way, while each of the elements of the set A either increases (in the case t1^{−}) or decreases (in the case t1^{+}) by 1.

2.4.2 G-steps

Here we suppose the coefficient ring I to be equal to K[n] or K[q^{n}] (see Example 1) and will call the elements of I polynomials. One can eliminate elements of the explicit matrix by, e.g., the Bareiss’ method, which gives an opportunity to factor out some polynomials from the matrix rows [11, 18]. But in the version of the EG-method described above the choice of the eliminating row is based on its width, not on the degree of its element in the corresponding column.

In fact the process of reducing the explicit matrix to the desired form (e.g., with a trapezoidal leading matrix which has no zero row) can be organized by alternating EG_{l}-steps (where one needs to consider the width of the rows) with G-steps (i.e., usual Gaussian elimination steps where one does not need to do that).

The process starts with a G-step. By Gaussian eliminations one transforms the explicit matrix P (having v rows) to P′ such that its leading part P′_{l} has the trapezoidal form, but its last rows can be zero. Suppose that in P′_{l} the zero rows have numbers k + 1, …, v, k ≤ v. We assume these last rows are not completely zero in the matrix P (otherwise one can drop these rows, decreasing v).

Then perform an EG_{l}-step. Apply to each of the rows with numbers k + 1, …, v transformation t1^{+} until a nonzero element appears in the leading part of the row. Then by EG_{l}-eliminations (i.e., being careful with the l-widths of the rows) get the matrix whose leading matrix is such that its rows with the numbers 1, …, k make up again a trapezoidal form matrix while succeeding rows do not have nonzero elements in the columns with the numbers 1, …, k + 1. If a row of this matrix is completely zero then drop the row decreasing v. If a row of the matrix has zero leading part, then apply t1^{+} to the row and EG_{l}-eliminate all the elements of the row that have numbers 1, …, k + 1 and so on. By Theorem 1 this leads to a matrix P′′ whose leading part P′′_{l} is such that

1) the matrix P′′_{l} has no zero row;

2) the first k + 1 rows of P′′_{l} make up a matrix of the trapezoidal form;

3) if i ≥ k + 1 and j ≤ k + 1 then p′′_{ij}(n) = 0.

Then again perform a G-step. One has to apply it to the rows with the numbers that are > k + 1. By this strategy we increase the number of the rows of the trapezoidal part of the leading matrix. Then again perform an EG_{l}-step and so on.

We can extend some of the G-steps. Let’s mark the rows which have only zero elements in the leading part after a G-step. If such a row was not marked before then after applying t1^{+} we use Gaussian eliminations in its leading part (it is an extension of the G-step). But if the row has been marked before then we use EG-eliminations. (This extension is optional.)

2.4.3 Sharpening the degree bound for A_{0}-solutions

Now we return to finding an upper bound for the degree of A_{0}-solutions of a recurrent system. If the trailing part is triangular with nonzero diagonal entries from the outset then we consider integer roots and use formula (13). But if the triangular form was obtained by a sequence of transformations t1^{−}, t2, t3 then it makes sense to invoke some additional reasoning. First some of the explicit matrix rows can have removed factors, i.e., be of the form (33), (34), and so on. Then the elements of the set A ∪ Ã ∪ … should be considered together with the integer roots of the determinant of the trailing part of the explicit matrix. It gives us a set U = {u_{1}, …, u_{τ} }, u_{1} > u_{2} > ⋯ > u_{τ} ≥ 0. Set n_{0} = −1 if U is empty and n_{0} = u_{1} otherwise. Now formula (13) can be used. But the bound given by (13) can be excessive. The reason is the following. Let U be nonempty. Then the substitution n = n_{0} = u_{1} into the recurrent system transforms it into a system S_{n_{0}} of linear algebraic equations for z_{n_{0}+t}^{1}, …, z_{n_{0}+t}^{m}. If the rank of the system is < m then we cannot uniquely determine z_{n_{0}+t}^{1}, …, z_{n_{0}+t}^{m} via z^{1}, …, z^{m} with larger indices using S_{n_{0}}. Add to S_{n_{0}} the constraints that were obtained in the process of the transformation of the trailing part into the triangular form. This gives a linear algebraic system S′_{n_{0}}. We substitute into S′_{n_{0}} zeros for all z^{1}, …, z^{m} with the indices > n_{0} + t. If the rank of the trailing part of the new system S′′_{n_{0}} is equal to m then u_{1} can be deleted from U. Then we can set n_{0} = u_{2} and so on until the moment when we either exhaust all elements of U which are > max{deg b + t, −1} or obtain a system whose rank is < m. In the first two cases the original recurrent system has no A_{0}-solution. In the third case we can use (13) with n_{0} = u_{λ} to find an upper bound.

Example 6

\begin{matrix}(n − 1)z_{n}^{1} − z_{n−3}^{2} = 0 \\ −2z_{n}^{1} + (n − 2)z_{n−1}^{2} = 0. \end{matrix}  (37)

The trailing part of the explicit matrix of the system

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

is equal to

( \begin{matrix}0 & −1 \\ 0 & 0 \end{matrix} ) .

Put the trailing part into triangular form by EG_{t}-eliminations:

( \begin{matrix}n − 1 & 0 & 0 & 0 & 0 & 0 & 0 & −1 \\ −2 & 0 & 0 & n − 2 & 0 & 0 & 0 & 0 \end{matrix} ) ⟶^{1} ( \begin{matrix}n − 1 & 0 & 0 & 0 & 0 & 0 & 0 & −1 \\ 0 & 0 & 0 & 0 & −2 & 0 & 0 & n − 4 \end{matrix} ) ⟶^{2}

( \begin{matrix}0 & 0 & 0 & 0 & −2 & 0 & 0 & n − 4 \\ n − 1 & 0 & 0 & 0 & 0 & 0 & 0 & −1 \end{matrix} ) ⟶^{3} ( \begin{matrix}0 & 0 & 0 & 0 & −2 & 0 & 0 & n − 4 \\ (n − 4)(n − 1) & 0 & 0 & 0 & 0 & 0 & 0 & −n + 4 \end{matrix} ) ⟶^{4}

( \begin{matrix}0 & 0 & 0 & 0 & −2 & 0 & 0 & n − 4 \\ 0 & 0 & (n − 5)(n − 2) & 0 & 0 & 0 & −2 & 0 \end{matrix} ) .

The determinant of the trailing part vanishes when n = 4 and therefore U = {4}. But the step ⟶^{3} used the multiplication of the second row by n − 4. This gives the constraint 3z_{4}^{1} − z_{1}^{2} = 0. We can add the row (0 − 1) to the matrix

( \begin{matrix}0 & n − 4 \\ −2 & 0 \end{matrix} )_{n=4} .

It gives the matrix

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

whose rank is equal to 2. We exclude the integer 4 from U. Now U is empty and formula (13) says that (37) has no A_{0}-solution of finite degree.

2.5 Right-hand sides in the form of sequences with nonnegative indices

To this point the assumption has been made that the right-hand sides b^{i} are double-sided sequences: b^{i} = {b_{n}^{i}}_{n∈Z} ∈ A_{0}, ν(b^{i}) = −∞. Recurrent relations are supposed to be satisfied for all n ∈ Z (or at least for almost all n as in the case (33)). In what follows we will consider the case ν(b^{i}) = 0 as well. In such a case transformation T1 (respectively transformations t1^{+}, t1^{−}) is not, in general, an equivalent transformation: let a recurrence

R_{i1}z^{1} + ⋯ + R_{im}z^{m} = b^{i}  (38)

belonging to system (15) be satisfied for n ≥ 0 and let E be applied to (38). For the new system to be equivalent to (15), the constraint which one gets after the substitution n = −1 into (38) must be supplied (if the i-th row of the explicit matrix of the new system has the form (33) and −1 ∈ A then replace A by A \ {−1}). Similarly, if E^{−1} was applied to (37) then there is a need to exclude the value n = 0 for the new recurrence (if the i-th row of the explicit matrix of the new system has the form (33) then replace A by A ∪ {0}) otherwise set A = {0} for the new row).

An approach to find an upper bound for the solution degree in the case where EG_{t}-eliminations gave the explicit matrix with the rows of the form (33) has been demonstrated in Section 2.4.3.

3 Solving systems of linear functional equations

3.1 Compatible bases

The methods described in the previous sections can be used for the investigation of linear functional equation systems and for finding their solutions in the form of polynomials and series. We mean in particular systems of linear differential, difference and q-difference equations with polynomial coefficients and solutions in the form of series in the bases such as

P = {x^{n}}_{n≥0},  (39)
C = \{\binom{x}{n}\}_{n≥0}.  (40)

We will use the general notion of the compatibility of an operator and a basis which was introduced in [4, 8]. In this section we briefly outline the needed facts concerning this theory (see [4, 8] for details).

Denote by L_{K[x]} the ring (the K-algebra) of linear operators L : K[x] → K[x]. Let B = {J_{n}(x)}_{n≥0} be a sequence of polynomials from K[x] such that

P1. deg J_{n} = n,

P2. J_{n} | J_{m} for 0 ≤ n < m.

From P1 it follows that {J_{0}, J_{1}, …} is a basis of K[x]. A basis B of K[x] satisfying P1, P2, and an operator L ∈ L_{K[x]} are compatible if there are nonnegative integers A, B, and elements α_{i,n} ∈ K for n ≥ 0 and −A ≤ i ≤ B, such that

LJ_{n} = ∑_{i=−A}^{B} α_{i,n} J_{n+i},  (41)

with J_{k} = 0 when k < 0.

We fix a basis B = {J_{n}(x)}_{n≥0} of K[x] having properties P1, P2. Any such a basis is compatible with multiplication by the independent variable.

If an operator L is compatible with B then L can be extended to K[[B]] i.e., to the space of formal series of the form

c_{0}J_{0}(x) + c_{1}J_{1}(x) + · · · .  (42)

Denote by L_{B} the set of operators L compatible with B. This set is a ring. Any L ∈ L_{B} induces a recurrent operator R_{B}L ∈ E, where E is the ring of the linear operators of the form (5), whose coefficients are arbitrary functions p_{i}(n) : Z → K. The operator R_{B}L is such that if y, f ∈ K[[B]] and y has the form (42) while f is a series

b_{0}J_{0}(x) + b_{1}J_{1}(x) + ⋯,  (43)

then Ly = f iff

R_{B}L({c_{n}}_{n∈Z}) = {b_{n}}_{n≥0},

where {c_{n}}_{n∈Z} is the sequence {c_{n}}_{n≥0} extended by taking c_{n} = 0 for all n < 0 (therefore {c_{n}}_{n∈Z} ∈ A_{0}). In this case the map

R_{B} : L_{B} → E

is a ring isomorphism.

Of the operators D, E_{x}, Q:

Dp(x) = \frac{d}{dx}p(x), E_{x}p(x) = p(x + 1), Qp(x) = p(qx),

the first and the third are compatible with P while the second one is compatible with C (see (39), (40)). Considering operator Q we assume K to be having a subfield K_{0} and an element q transcendental over K_{0} such that K = K_{0}(q). Let’s consider the rings K[x, D], K[x, E_{x}] and K[x, Q] i.e, the rings of differential, difference and q-difference operators with polynomial coefficients. To describe transformation R_{P} on K[x, Q], it suffices to give it on the two generators Q and x. It can be shown that

R_{P}Q = q^{n}, R_{P}x = E^{−1},  (44)

and, resp.,

R_{P}^{−1}q^{n} = Q, R_{P}^{−1}E^{−1} = x.

Therefore R_{P} is an isomorphism K[x, Q] onto K[q^{n}, E^{−1}]. It is convenient pass from K[x, D] to K[x, x^{−1}, D] because R_{P} is an isomorphism K[x, x^{−1}, D] onto K[n, E, E^{−1}]. The complete list of formulas is the following:

R_{P}x = E^{−1}, R_{P}x^{−1} = E, R_{P}D = (n + 1)E,
R_{P}^{−1}n = xD, R_{P}^{−1}E = x^{−1}, R_{P}^{−1}E^{−1} = x.

It is easy to see that Q and D are not compatible with C. In turn E_{x} is not compatible with P but it is compatible with C. We have

R_{C}E_{x} = E + 1, R_{C}x = n(1 + E^{−1});

E_{x}^{−1} and C are not compatible, but we can easily find R_{C}^{−1}M ∈ K[x, E_{x}, E_{x}^{−1}] if M ∈ K[n, E], using

R_{C}^{−1}E = E_{x} − 1, R_{C}^{−1}n = x(1 − E_{x}^{−1}).

3.2 Functional equations systems

The papers [4, 8] were devoted to the case of a scalar equation, i.e., to the case of a single equation Ly = b with one unknown y. But similarly to the scalar case the following theorem can be easily proven

Theorem 2 Let L_{1}, …, L_{m} ∈ L_{B} and b, y_{1}, …, y_{n} ∈ K[[B]]. Let b be of the form (43) and

y_{i} = z_{0}^{i}J_{0}(x) + z_{1}^{i}J_{1}(x) + ⋯,  (45)

i = 1, …, m. Then

L_{1}y_{1} + ⋯ + L_{m}y_{m} = b  (46)

iff

R_{B}L_{1}({z_{n}^{1}}_{n∈Z}) + ⋯ + R_{B}L_{m}({z_{n}^{m}}_{n∈Z}) = {b_{n}}_{n≥0},  (47)

where {z_{n}^{i}}_{n∈Z} are the sequences {z_{n}^{i}}_{n≥0} extended by taking z_{n}^{i} = 0 for all n < 0 (therefore {z_{n}^{i}}_{n∈Z} ∈ A_{0}, for i = 1, …, m).

The proof, similarly to the scalar case, is based on the observation that if L ∈ L_{B} and y_{i} = ∑_{n=0}^{∞} z_{n}^{i}j_{n} ∈ K[[B]] then the equality

L ∑_{n=0}^{∞} z_{n}^{i}J_{n} = ∑_{n=0}^{∞}((R_{B}L)z^{i})_{n}J_{n}

takes place.

The equivalence of (46) and (47) implies that if L_{ij} ∈ L_{B} ,

f_{i} = ∑_{n=0}^{∞} b_{n}^{i}J_{n}, y_{j} = ∑_{n=0}^{∞} z_{n}^{j}J_{n},

i = 1, …, v, j = 1, …, m, then

( \begin{matrix}L_{11} & … & L_{1m} \\ ⋮ & & ⋮ \\ L_{v1} & … & L_{vm} \end{matrix} ) ( y^{1} ; ⋮ ; y^{m} ) = ( f^{1} ; ⋮ ; f^{v} ) ,  (48)

iff

( \begin{matrix}R_{11} & … & R_{1m} \\ ⋮ & & ⋮ \\ R_{v1} & … & R_{vm} \end{matrix} ) ( z^{1} ; ⋮ ; z^{m} ) = ( b^{1} ; ⋮ ; b^{v} ) ,  (49)

where R_{ij} = R_{B}L_{ij} and z^{j} = {z_{n}^{j}}_{n∈Z}, b^{i} = {b_{n}^{i}}_{n≥0}, z^{j} ∈ A_{0}.

We reduce in this way the problem of solving a system of linear functional equations in the class of series to the problem of solving a recurrent system. If we consider differential and difference operators with coefficients from K[x] then we get recurrences with coefficients from I = K[n]. In turn q-difference operators give recurrences with coefficients from I = K[q^{n}].

We point out an essential distinction between the differential and q-difference cases on the one hand and the difference case on the other. In the differential and q-difference cases we use the basis P to which the negative powers of x can be added: J_{−1} = x^{−1}, J_{−2} = x^{−2}, … Correspondingly, the right-hand sides of the original systems and solutions of the systems can be assumed to belong to the class of double-sided power series of the form

⋯ + z_{−2}x^{−2} + z_{−1}x^{−1} + z_{0} + z_{1}x + z_{2}x^{2} + ⋯  (50)

This leads to the same recurrent systems of the form (49) but with ν(b^{j}) = −∞, j = 1, …, v. The recurrences of such systems are satisfied for all n ∈ Z (or for almost all if the rows of the explicit matrix have the form (33)). Let b^{j}, z^{j} ∈ A_{0}, j = 1, …, v. Apply E^{−1} to a recurrence of the system (or t1^{−} to a row of the explicit matrix). Then it is possible not to exclude the case n = 0 for the new recurrence (i.e., not to add n to A in (33)). The reason is that if we do not exclude it then we in fact add a constraint to our system, but this constraint is necessarily true.

The simple example

y(x + 1) − y(x) = 0  (51)

shows that the difference case is not so nice. The recurrence for coefficients of a series

z_{0} + z_{1}x + z_{2} \frac{x(x − 1)}{2} + ⋯

satisfying (51) is z_{n+1} = 0. This recurrence has to be satisfied for n = 0, 1, … The general solution is

z_{n} = \{ \begin{matrix}C, & if n = 0, \\ 0, & otherwise, \end{matrix}

where C is an arbitrary constant. But we cannot consider the recurrence for all n because in such a case we have only the identically zero solution.

3.3 Independence of equations

Lemma 4 If the recurrences corresponding to a linear differential, difference or q-difference system with polynomial coefficients are dependent over E_{I} (I = K[n] in the differential and difference cases and I = K[q^{n}] in the q-difference case) then the equations of the original system are dependent over, resp., K[x, x^{−1}, D],K[x, E_{x}, E_{x}^{−1}] or K[x, Q].

Proof : trivial in the differential and q-difference cases since if (16) together with R_{ij} = R_{P}L_{ij} hold then applying R_{P}^{−1} to all S_{i} we will obtain S̃_{i} such that

S̃_{1} ∘ L_{1j} + ⋯ + S̃_{v} ∘ L_{vj} = 0,  (52)

j = 1, …, m. In the difference case we deal with R_{C}; we know that R_{C}^{−1}S_{i} is not necessary a difference operator because we can face with E with a negative exponent in S_{i}. Let w be the lowest exponent of E in S_{1}, …, S_{v}. We obtain

E^{−w} ∘ S_{1} ∘ R_{1j} + ⋯ + E^{−w} ∘ S_{v} ∘ R_{1v} = 0,

j = 1, …, m. If we now set S̃_{i} = R_{C}^{−1}(E^{−w} ∘ S_{i}) (a difference operator) then we get (52). □

We proved earlier Lemma 3 which claims that the equations of recurrent system (17) are independent over E_{I}. An analogous proposition can be proven for systems of the form

p_{1}(x)y′_{1}(x) = a_{11}(x)y_{1}(x) + ⋯ + a_{1m}(x)y_{m}(x) + b_{1}(x)
· · · · · · · · · · · ·  (53)
p_{v}(x)y′_{v}(x) = a_{v1}(x)y_{1}(x) + ⋯ + a_{vm}(x)y_{m}(x) + b_{v}(x),

p_{1}(x)y_{1}(x + 1) = a_{11}(x)y_{1}(x) + ⋯ + a_{1m}(x)y_{m}(x) + b_{1}(x)
· · · · · · · · · · · ·  (54)
p_{v}(x)y_{v}(x + 1) = a_{v1}(x)y_{1}(x) + ⋯ + a_{vm}(x)y_{m}(x) + b_{v}(x),

p_{1}(x)y_{1}(qx) = a_{11}(x)y_{1}(x) + ⋯ + a_{1m}(x)y_{m}(x) + b_{1}(x)
· · · · · · · · · · · ·  (55)
p_{v}(x)y_{v}(qx) = a_{v1}(x)y_{1}(x) + ⋯ + a_{vm}(x)y_{m}(x) + b_{v}(x),

where v ≤ m and p_{1}, …, p_{v} ≠ 0, using the same approach as in Lemma 3 (where in differential and q-difference cases ord must be used instead of ord^{∗}). We have therefore

Theorem 3 Let a linear recurrent, differential, difference or q-difference system of the first order with polynomial coefficients be given. Let the system be in the canonical form, resp., (17), (53), (54) or (55). Let the number of the equations in this system do not exceed the number of the unknowns. Then the equations of the system are independent over, resp., E_{I}, K[x, x^{−1}, D], K[x, E_{x}, E_{x}^{−1}] or K[x, Q].

Theorem 3 and Lemma 4 show that if we apply EG-eliminations to solve a system of the form (53),

(54) or (55) then it is not possible to get zero rows in the explicit matrix; if additionally the number of the equations in this system is equal to the number of the unknowns then we can find an upper bound for the degrees of polynomial solutions of each of the systems. We can also construct all the polynomial solutions and describe the set of all solutions in the form of series in the basis P (the differential and q-difference case) or in the basis C (the difference case).

3.4 Laurent series; rational solutions of differential systems

Considering recurrent systems we discussed the existence of their solutions in A_{0}. We can consider together with A_{0} the classes A_{−1}, A_{−2}, …:

c ∈ A_{k} ⇔ (c_{i} ≠ 0 ⇒ i ≥ k),

k = 0, −1, −2, … (we have noted that in the differential and q-difference cases we can consider the recurrences of system (49) for all n). It is obvious that A_{0} ⊂ A_{−1} ⊂ A_{−2} ⊂ · · ·. Set A_{≤0} = A_{0} ∪ A_{−1} ∪ A_{−2} ∪ … All constructive processes described above for A_{0} can be performed for A_{−1}, A_{−2}, … (but if we deal with A_{−k} and consider the integer roots of the determinant of the trailing part of the explicit matrix then we have to take into account all the roots ≥ −k − l instead of the roots ≥ −l as in the case of A_{0}). The only question which must be previously answered is: let system (8) have the right-hand side in A_{k}^{m}, where k is a fixed non-positive integer; how to find a non-positive k_{0} such that if the system has a solution in A_{≤0}^{m}, then this solution is in A_{k_{0}}^{m}? If the leading part of the explicit matrix of (8) is square and non-singular then the question can be answered similarly to the question on the upper bound of the degree of solution (considered in Section 2.1). Similarly to Lemma 1 the following lemma can be proven:

Lemma 5 Let v = m in system (8) and the right-hand side of (8) is in A_{s}, s ≤ 0. Let p(n) = det P_{l}(n) be a nonzero element of I. Let n_{1} be the lowest integer root of p(n) = 0 if such roots exist and n_{1} = 1 otherwise. Let c ∈ A_{≤0}^{m} be a solution of (8). Then c ∈ A_{k_{0}}^{m}, where k_{0} ≥ min{n_{1} + l, s + l, 0}.

If the leading part of the explicit matrix is singular then we use EG_{l}-eliminations to make it nonsingular. We can use the approach described in 2.4.3 to compute a bound more accurately. Of course, we have to deal with the leading part instead of the trailing one.

The basis P can be extended by negative powers of x, and the approach just mentioned can be applied to search for Laurent series solutions of the differential and q-difference systems (near zero). The problem of finding Laurent or Taylor series solutions near a fixed point a can be reduced in the differential case to that at the point 0 by changing any operators

L_{ij} = p_{ijd}(x)D^{d} + ⋯ + p_{ij0}(x)

of (16) by

L_{ij}^{a} = p_{ijd}(x + a)D^{d} + ⋯ + p_{ij0}(x + a)  (56)

and, resp., taking f_{i}(x + a) instead of f_{i}(x) in the right-hand side, i = 1, …, v; j = 1, …, m. It gives in particular rational function solutions of a linear differential system with polynomial coefficients and polynomial right-hand sides if we know all the singularities of the system. The set of the singularities of a first order system written in the canonical form (53) is a subset of the set of the roots of the polynomials p_{1}(x), …, p_{v}(x). For every of these roots we can find (as was explained above) an upper bound for the order of the pole at this point. It allows one to obtain a common denominator U(x) of all rational function solutions of the system. Then one can substitute

y_{i}(x) = \frac{ỹ_{i}(x)}{U(x)},  (57)

i = 1, …, m, into (16) and find the polynomial solutions of the obtained system.

3.5 Rational solutions of difference and q-difference systems

The search for rational solutions of difference system has been considered in [3, 20]. As in the differential case one can construct a universal denominator U(x) and then reduce the search to the problem of finding polynomial solutions. The last problem can be solved using EG_{t}-eliminations (in [3, 20] the usage of the super-irreducible form of the given system was considered).

The methods [3, 20] are applicable to the q-difference case, excepting the situation where the denominators of some of y_{i}(x) are divisible by x, and, consequently, U(x) also has to be divisible by x.

The general case was considered in [2], where the difference approach was combined with the search for an upper bound for the pole order at x = 0. EG_{l}-eliminations allow one to find such a bound.

3.6 Examples of correspondence of functional equation systems and recurrent systems

Those examples of recurrent systems that we considered above can be interpreted as corresponding to linear systems of differential and difference equations with polynomial coefficients (in the differential case the polynomial can be Laurentian).

Example 7 The recurrent system (26) corresponds to the differential system

\begin{matrix}y′_{1} = y_{1} + y_{2} \\ y′_{2} = y_{1} + y_{2} \end{matrix}

and to the difference system

\begin{matrix}xy_{1}(x + 1) − (2x + 1)y_{1}(x) − y_{2}(x) + xy_{1}(x − 1) = 0 \\ xy_{2}(x + 1) − y_{1}(x) − (2x + 1)y_{2}(x) + xy_{2}(x − 1) = 0. \end{matrix}

The reasoning given in Example 3 shows that the polynomial solutions of these systems are y_{1}(x) = −y_{2}(x) = C, where C is an arbitrary constant. Note that in the difference case we have to pay attention to transformation ⟶^{2} in Example 3. It is transformation t1^{−} and the value n = 0 should be considered as a candidate for the solution degree. But the determinant of the trailing part of the matrix after transformation ⟶^{3} has the root 0 and therefore this value is already under consideration as a possible solution degree.

In a similar manner, the recurrent systems in Examples 4 and 5 can be considered as corresponding to some differential and difference systems. The correspondence to difference systems is possible since these systems can be written without the operator E with a negative exponent.

Example 8 The recurrent system considered in Example 6 uses E^{−1} and E^{−3}. We can consider it as corresponding to the differential system

\begin{matrix}xy′_{1} − y_{1} − x^{3}y_{2} = 0 \\ x^{2}y′_{2} − 2y_{1} − xy_{2} = 0. \end{matrix}  (58)

This system was considered by M.A.Barkatou [13] and was reduced to the super-irreducible form (this form is convenient, e.g., for recognizing existence of polynomial solutions) by a special change of independent variables y_{1}, y_{2}. In Example 6 we established the non-existence of polynomial solutions by EG_{t}-eliminations. When looking for power series solutions of (58) one can use the EG_{l}-eliminations and transform the leading part of the explicit matrix to the triangular form:

( \begin{matrix}n − 1 & 0 & 0 & 0 & 0 & 0 & 0 & −1 \\ −2 & 0 & 0 & n − 2 & 0 & 0 & 0 & 0 \end{matrix} ) ⟶ ( \begin{matrix}−2 & 0 & 0 & n − 2 & 0 & 0 & 0 & 0 \\ n − 1 & 0 & 0 & 0 & 0 & 0 & 0 & −1 \end{matrix} ) ⟶

( \begin{matrix}−2 & 0 & 0 & n − 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & (n − 1)(n − 2) & 0 & 0 & 0 & −2 \end{matrix} ) ⟶ ( \begin{matrix}−2 & 0 & 0 & n − 2 & 0 & 0 & 0 & 0 \\ 0 & n(n − 1) & 0 & 0 & 0 & −2 & 0 & 0 \end{matrix} ) .

Using the obtained matrix we can write down the recurrent system in the form

( \begin{matrix}−2 & 0 \\ 0 & n(n − 1) \end{matrix} ) ( \begin{matrix}z_{n}^{1} \\ z_{n}^{2} \end{matrix} ) + ( \begin{matrix}0 & n − 2 \\ 0 & 0 \end{matrix} ) ( \begin{matrix}z_{n−1}^{1} \\ z_{n−1}^{2} \end{matrix} ) + ( \begin{matrix}0 & 0 \\ 0 & −2 \end{matrix} ) ( \begin{matrix}z_{n−2}^{1} \\ z_{n−2}^{2} \end{matrix} ) = ( \begin{matrix}0 \\ 0 \end{matrix} ) .  (59)

The largest integer root of the determinant of the leading part of the explicit matrix is equal to 1, it is sufficient to find z_{0}^{1}, z_{0}^{2}, z_{1}^{1}, z_{1}^{2}. The substitution z_{−2}^{1} = z_{−2}^{2} = z_{−1}^{1} = z_{−1}^{2} = 0 into (59) with n = 0, 1 gives −2z_{0}^{1} = 0, −2z_{1}^{1} = z_{0}^{2}. Therefore z_{1}^{1}, z_{1}^{2} can be taken arbitrary, while z_{0}^{1} = 0, z_{0}^{2} = 2z_{1}^{1}. The successive coefficients of the series

y_{1} = z_{0}^{1} + z_{1}^{1}x + z_{2}^{1}x^{2} + ⋯, y_{2} = z_{0}^{2} + z_{1}^{2}x + z_{2}^{2}x^{2} + ⋯

can be computed by (59). Due to

( \begin{matrix}−2 & 0 \\ 0 & n(n − 1) \end{matrix} )^{−1} = ( \begin{matrix}−\frac{1}{2} & 0 \\ 0 & \frac{1}{n(n−1)} \end{matrix} )

we have the formula

( \begin{matrix}z_{n}^{1} \\ z_{n}^{2} \end{matrix} ) = ( \begin{matrix}0 & \frac{n−2}{2} \\ 0 & 0 \end{matrix} ) ( \begin{matrix}z_{n−1}^{1} \\ z_{n−1}^{2} \end{matrix} ) + ( \begin{matrix}0 & 0 \\ 0 & \frac{2}{n(n−1)} \end{matrix} ) ( \begin{matrix}z_{n−2}^{1} \\ z_{n−2}^{2} \end{matrix} )

or

z_{n}^{1} = \frac{n − 2}{2}z_{n−1}^{2}, z_{n}^{2} = \frac{2}{n(n − 1)}z_{n−2}^{2}

for n = 2, 3, … Example 9 Consider now the q-difference system

\begin{matrix}qxG_{0}(q^{3}x) = G_{0}(x) − G_{1}(x) \\ qx^{2}G_{1}(q^{3}x) = G_{1}(x) − G_{2}(x) \\ qx^{3}G_{3}(q^{3}x) = G_{2}(x) − G_{3}(x) \\ G_{0}(q^{3}x) = G_{3}(x). \end{matrix}  (60)

which was originally considered in [10] in connection with a partitions theory problem. The rewriting of the system in the operator form (48) leads to

( \begin{matrix}qxQ^{3} − 1 & 1 & 0 & 0 \\ 0 & qx^{2}Q^{3} − 1 & 1 & 0 \\ 0 & 0 & −1 & qx^{3}Q^{3} + 1 \\ Q^{3} & 0 & 0 & −1 \end{matrix} ) ( \begin{matrix}G_{0} \\ G_{1} \\ G_{2} \\ G_{3} \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ 0 \\ 0 \end{matrix} ) .

Applying R_{P} to the operator matrix gives

( \begin{matrix}q^{3n−2}E^{−1} − 1 & 1 & 0 & 0 \\ 0 & q^{3n−5}E^{−2} − 1 & 1 & 0 \\ 0 & 0 & −1 & q^{3n−8}E^{−3} + 1 \\ q^{3n} & 0 & 0 & −1 \end{matrix} ) .

The explicit matrix is

( \begin{matrix}−1 & 1 & 0 & 0 & q^{3n−2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & −1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & q^{3n−5} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & −1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & q^{3n−8} \\ q^{3n} & 0 & 0 & −1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix} ) .  (61)

Let’s investigate the existence of polynomial solutions of this system. The determinant of the trailing part of (61) is zero. EG_{t}-eliminations give us

( \begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & q^{3n−9} & 0 & −1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & −q^{3n−9} & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & −1 & 1 & 0 & 0 & 0 & 0 & q & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & −1 & 1 & 0 & 0 & 0 & 0 & 0 & q^{3n−8} \end{matrix} )

The determinant of the trailing part of the last matrix is −q^{9n−25}. The equation q^{9n−25} = 0 has no integer root. Therefore the system (60) has no polynomial solution. Observe that the determinant of the leading part of (61) is 1 − q^{3n} with integer root n = 0. This and the equalities l = 0, t = −2 show that (60) has a power series solution. This solution is unique up to a constant factor since the rank of the leading part of (61) is equal to 3 when n = 0. This solution was found in a different way (by using the “right” substitution) in [10].

Example 10 Now we construct all rational solutions of the differential system:

\begin{matrix}(x^{2} − 100)y′_{1} + 4xy_{1} + 2y_{2} = 0 \\ y′_{2} − y_{1} = 0 \end{matrix}  (62)

with singularities −10, 10. This system transformed by means of (56) with a = −10 looks like the following:

\begin{matrix}(x^{2} − 20x)y′_{1} + (4x − 40)y_{1} + 2y_{2} = 0 \\ y′_{2} − y_{1} = 0. \end{matrix}

The corresponding recurrent system has the explicit matrix

( \begin{matrix}0 & 0 & −20 n − 40 & 2 & n + 3 & 0 \\ 0 & n + 1 & −1 & 0 & 0 & 0 \end{matrix} )

with l = 1, t = −1. EG_{l}-eliminations give

( \begin{matrix}0 & n + 1 & −1 & 0 & 0 & 0 \\ −20 n − 60 & 2 & n + 4 & 0 & 0 & 0 \end{matrix} ) .

The equation (n + 1)(−20n − 60) = 0 has roots −3, −1, and since −3 + l = −3 + 1 = −2 the order of the solution pole at a = −10 is bounded by 2.

Analogously the same bound can be obtained in the case a = 10. We have the common denominator U(x) = (x−10)^{2}(x+10)^{2} for all rational solutions of (62). Substitution (57) gives the differential system

\begin{matrix}(x^{2} − 100)y′_{1} + 2ỹ_{2} = 0 \\ (x^{2} − 100)y′_{2} + (x^{2} − 100)ỹ_{1} − 4xỹ_{2} = 0. \end{matrix}  (63)

The corresponding recurrent system has the explicit matrix

( \begin{matrix}−100 n − 100 & 0 & 0 & 2 & n − 1 & 0 & 0 & 0 \\ 0 & −100 n − 100 & 100 & 0 & 0 & n − 5 & −1 & 0 \end{matrix} )

with l = 1, t = −2. EG_{t}-eliminations give

( \begin{matrix}0 & 0 & −100 n & 0 & 0 & 2 & n − 2 & 0 \\ 0 & 0 & 0 & 100 (n − 3) n & 200 & 0 & 0 & −n^{2} + 9n − 20 \end{matrix} )

with the additional constraint in the case n = 2:

−300z_{3}^{2} + 100z_{2}^{1} − 3z_{1}^{2} − z_{0}^{1} = 0.

The equation (−n^{2} + 9n − 20)(n − 2) = 0 has the roots 2, 4, 5. So the degree of any polynomial solution of the system is bounded by 5 + t = 5 + (−2) = 3.

The original explicit matrix of the system has the leading part already in the triangular form and could be used for computing coefficients of the polynomial solutions. It leads to the following general polynomial solution of (63):

[−\frac{C_{1}}{100}x^{3} − \frac{C_{2}}{2} x^{2} + C_{1}x + 50C_{2}, \frac{C_{1}}{100}x^{2} + C_{2}x + C_{1}] ,

and correspondingly to the following general rational solution of (62):

[−\frac{1}{100} \frac{C_{1}x + 50C_{2}}{(x + 10) (x − 10)}, \frac{1}{100} \frac{C_{1}x^{2} + 100C_{2}x + 100C_{1}}{(x − 10)^{2} (x + 10)^{2}} ]

where C_{1} and C_{2} are arbitrary constants.

4 Shifts of unknown sequences in recurrent systems

Let the i-th column of the leading matrix of a linear recurrent system (8) be zero, 1 ≤ i ≤ m. Then one can substitute z^{i} = Ez^{i} without increasing l or decreasing t in the system. Similarly if the i-th column of the trailing matrix of the system is zero then one can substitute z^{i} = E^{−1}z^{i}. Going back to Example 6 we see that it is possible to set instantaneously z^{1} = E^{−3}z^{1} (thus z_{n}^{1} = z_{n−3}^{1}). The system (37) will be rewritten in the form

\begin{matrix}(n − 1)z_{n−3}^{1} − z_{n−3}^{2} = 0 \\ −2z_{n−3}^{1} + (n − 2)z_{n−1}^{2} = 0, \end{matrix}

and the explicit matrix will be equal to

( \begin{matrix}0 & 0 & 0 & 0 & n − 1 & −1 \\ 0 & n − 2 & 0 & 0 & −2 & 0 \end{matrix} ) ,

l = −1, t = −3. The determinant of the trailing part of the explicit matrix is equal to −2. In particular this means that the original system has no polynomial solutions.

Now consider the leading part of the explicit matrix of system (37). After substitution z^{2} = Ez^{2}, i.e., z_{n}^{2} = z_{n+1}^{2}, the system (37) will be rewritten in the form

\begin{matrix}(n − 1)z_{n}^{1} − z_{n−2}^{2} = 0 \\ −2z_{n}^{1} + (n − 2)z_{n}^{2} = 0, \end{matrix}

and the explicit matrix will be equal to

( \begin{matrix}n − 1 & 0 & 0 & 0 & 0 & −1 \\ −2 & n − 2 & 0 & 0 & 0 & 0 \end{matrix} ) ,

l = 0, t = −2. The leading part of the explicit matrix is a non-singular matrix.

Certainly such transformations of the explicit matrix can be done without rewriting the original system.

In the examples given in this section we did not use EG-eliminations at all. In general, the combination of shifts and EG-eliminations seems to be useful.

Acknowledgement

I would like to thank M. Bronstein, D. E. Khmelnov and H. Q. Le for their discussions about this problem.

References

[1] S.A. Abramov (1995): Rational solutions of linear difference and q-difference equations with polynomial coefficients, Programming and Comput. Software 21, No 6, 273–278. Transl. from Programmirovanie, No 6, 3–11.

[2] S.A. Abramov (1999): Rational solutions of first order q-difference systems, submitted to FPSAC’99.

[3] S.A. Abramov, M.A. Barkatou (1998): Rational solutions of first order difference systems, Proc. IS-SAC’98, 124–131.

[4] S. Abramov, M. Bronstein, M. Petkovšek (1995): On polynomial solutions of linear operator equations, Proc. ISSAC’95, 290–296.

[5] S. A. Abramov, P. Paule, M. Petkovšek (1998): q-Hypergeometric solutions of q-difference equations, Discrete Math. 180 (1998) 3–22.

[6] S.A. Abramov, M. Petkovšek (1994): D’Alembertian solutions of linear differential and difference equations, Proc. ISSAC’94, 169–174.

[7] S. Abramov, M. Petkovšek (1996): Special power series solutions of linear differential equations, Proc. FPSAC’96, 1–7.

[8] S. Abramov, M. Petkovšek, A. Ryabenko (1996): Special formal series solutions of linear operator equations, Discrete Math. (to appear).

[9] S.A. Abramov, E.V. Zima (1997): A universal program to uncouple linear systems, Proc. of CMCP’96 (International Conference on Computational Modeling and Computing in Physics, Dubna, Russia, Sept. 16-21, 1996)7, 16–26.

[10] G.E. Andrews (1976): The Theory of Partitions, Encyclopedia of Mathematics and its Applications, Vol. 2, Addison-Wesley, Reading, Mass.

[11] E.H. Bareiss (1968): Sylvester’s identity and multistep integer-preserving Gaussian elimination, Math. Comp. 22, 565–578.

[12] M.A. Barkatou (1989): Contribution à l’étude des équations différentielles et aux différences dans le champ complexe, Thèse soutenue le 6 juin 1989 à l’Institut national polytechnique de Grenoble.

[13] M.A. Barkatou (1997): On rational solutions of systems of linear differential solutions, Journal of Symbolic Computation (to appear).

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

[15] F. Chyzak (1997): An extension of Zeilberger’s fast algorithm to general holonomic functions, Proc. FPSAC’97 , 172–183.

[16] F. Chyzak (1998): Gröbner bases, symbolic summation and symbolic integration, in: Gröbner Bases and Applications, B. Buchberger and F. Winkler, Eds., London Mathematical Society Lecture Note Series 251, Cambridge University Press, Cambridge UK, 1998, 32–60.

[17] F. Chyzak (1998): Fonctions holonomes en calcul formel, Thèse d’informatique, Ecole Polytechnique, Palaiseau.

[18] J. H. Davenport, Y. Siret, E. Tournier (1988): Computer Algebra, Academic Press.

[19] A. Hilali, A. Wazner (1987): Formes super-irréductibles des systèmes différentiels linéaires, Num. Math. 50, 429–449.

[20] M. van Hoeij (1998): Rational solutions of linear difference equations, Proc. ISSAC’98, 120–123.

[21] M. van Hoeij, J.-F. Ragot, F. Ulmer, J.-A. Weil (1998): Liouvillian solutions of linear differential equations of order three and higher, Journal of Symbolic Computation (to appear).

[22] M. van Hoeij, J.-A. Weil (1997): An algorithm for computing invariants of differential Galois groups, J. Pure and Applied Alg. 117 & 118, 353–379.

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