Improved Polynomial Remainder Sequences for Ore Polynomials

Polynomial remainder sequences contain the intermediate results of the Euclidean algorithm when applied to (non-)commutative polynomials. The running time of the algorithm is dependent on the size of the coefficients of the remainders. Different ways have been studied to make these as small as possible. The subresultant sequence of two polynomials is a polynomial remainder sequence in which the size of the coefficients is optimal in the generic case, but when taking the input from applications, the coefficients are often larger than necessary. We generalize two improvements of the subresultant sequence to Ore polynomials and derive a new bound for the minimal coefficient size. Our approach also yields a new proof for the results in the commutative case, providing a new point of view on the origin of the extraneous factors of the coefficients.


Introduction
When given a system of differential equations, one might be interested in finding the common solutions of these equations. In order to do so, one can compute another differential equation whose solution space is the intersection of the solution spaces of the equations in the original system. One way to do this is to translate the equations into operators and use the Euclidean algorithm to compute their greatest common right divisor. The solution space of the greatest common right divisor then consists of the desired elements.
Similarly, given a sequence of numbers (t n ) n∈{0,1,... } that satisfies two different recurrence equations, the Euclidean algorithm is used in applications to find a reasonable candidate for the least order equation of which (t n ) n∈{0,1,... } is a solution.
Carrying out Euclid's algorithm applied to two polynomials over a domain usually requires a prediction of the denominators that might appear in the coefficients of the remainders in order to bypass costly computations in the quotient field of . While such a prediction can be done easily, the growth of the coefficients of the remainders can be tremendous, which might result in an unnecessary high running time. This can be avoided by dividing out possible content of the remainders to make their coefficients as small as possible. For commutative polynomials as well as for non-commutative operators, different ways have been extensively studied to find factors of the content in the sequence of remainders without computing the GCD of the coefficients of each element of the sequence. Most notably in this respect are subresultant sequences, where the growth of the coefficients can be reduced from exponential to linear in the number of reduction steps in the Euclidean algorithm. When taking generic, randomly generated input, the coefficient size in the subresultant sequence is usually optimal, but when taking the input from applications in e.g. combinatorics or physics, the remainders still have non-trivial content in many cases.
For commutative polynomials, some ways are known to improve on subresultants. In this article we generalize two of these results to Ore polynomials and we also give a new proof for the commutative case that is based on the structure of subresultants as matrix determinants. Furthermore, we use these results to derive a new bound for the coefficient size of the content-free remainders.
In Section 2 the basic notions of Ore polynomial rings are stated. A precise definition and examples of polynomial remainder sequences are given in Section 3 and further details on the subresultant sequence are then presented in Section 4. The main results of this article can be found in Sections 5 and 6, where we first describe how additional content in the subresultant sequence can emerge and then use these results to improve on the Euclidean algorithm and to get a new bound for the size of the coefficients.

Preliminaries
The algebraic framework for different kinds of operators that we consider here are Ore polynomial rings, which were introduced by Øystein Ore in the 1930's. We provide an overview of some basic facts that suffice our needs and that can be found in Ore (1933) and Bronstein and Petkovšek (1996).
2. Suppose that δ is a pseudo-derivation w.r.t. σ. We define the Ore polynomial ring ( [x], +, ·) with componentwise addition and the unique distributive and associative extension of the multiplication rule xa = σ(a)x + δ(a) for any a ∈ , to arbitrary polynomials in [x]. To clearly distinguish this ring from the standard polynomial ring over , we denote it by [x; σ, δ].
Elements of an Ore polynomial ring are called operators and are denoted by capital letters. We refer to the leading coefficient of an operator A as lc(A), to the coefficient of x 0 in A as tc(A) and to the polynomial degree of A in x as the order d A of A.
In this article, we consider the following situation: Let be a Euclidean domain with degree function deg and let [x; σ, δ] be an Ore polynomial ring where σ is an automorphism. For any operator A ∈ [x; σ, δ], we define A to be the maximal coefficient degree of A. The content cont(A) of A is the greatest common divisor of all the coefficients of A and it is defined to be lc(A) if is a field. It is possible to extend [x; σ, δ] to an Ore polynomial ring over the quotient field Ã of by setting σ(a −1 ) = σ(a) −1 and δ(a/b) = (bδ(a) − aδ(b))/(bσ(b)) for a, b ∈ , b = 0 (see Li (1996) Definition 2. For a ∈ and n ∈ AE, σ n (a) is obtained by applying n times σ to a and σ −n (a) := (σ −1 ) n (a), where σ −1 is the inverse map of σ. The nth σ-factorial of a ∈ is defined as the product

Polynomial Remainder Sequences for Ore Polynomials
The greatest common right divisor of A and B can be computed by using the Euclidean algorithm. If we multiply any intermediate result that appears during the execution of the algorithm by an element of Ã \ {0}, the final output will be a Ã-left multiple of G. This amount of freedom allows us to optimize the running time by choosing these factors appropriately. In order to be able to formulate improvements of this kind, the notion of polynomial remainder sequences has been introduced. Each element of such a sequence corresponds to a remainder computed in one iteration of the Euclidean algorithm.
A PRS of A and B is uniquely determined by specifying the α i and β i . Whenever we talk about a PRS (R i ) i∈{0,...,ℓ+1} , we allow ourselves to refer to the related sequences (Q i ) i∈{1,...,ℓ} , (d i ) i∈{0,...,ℓ} etc. as in the above definition without explicitly introducing them.
In order to efficiently compute G, one wants to make sure that all the remainders are elements of [x; σ, δ] rather than Ã[x;σ,δ]. This can be achieved by choosing the α i in a way such that the quotient of any two consecutive remainders has coefficients in . To this extent, for 1 ≤ i ≤ ℓ set and division with remainder yields Q i and R i+1 in [x; σ, δ] with: We call pquo(R i−1 , R i ) := Q i the pseudo-quotient of R i−1 and R i and prem(R i−1 , R i ) := R i+1 the pseudo-remainder of R i−1 and R i . The α i are used to make sure that computations can be done in [x; σ, δ] and the β i control the coefficient growth in a PRS. We want β i to contain as many factors of the content of R i+1 as possible without much computational overhead needed to obtain these factors.
1. β i = 1. This is called the pseudo PRS of A and B. Here, no content will be divided out. 2. β i = cont(R i+1 ). This is called the primitive PRS of A and B. The coefficients of the remainders will be as small as possible, but it is necessary to compute the GCD of the coefficients of each remainder in order to get the β i .

The subresultant PRS of A and B (see Section 4) is given by
In this PRS, the content that is generated systematically by pseudo-remaindering will be cleared from the remainders.
While in all of the above PRSs the remainders are elements of [x; σ, δ], the degrees of the coefficients differ drastically, as illustrated in the following example. It can be shown that the degrees of the coefficients in the pseudo PRS grow exponentially with i, which renders this PRS practically useless. The growth in the subresultant and primitive PRS is linear in i.
Example 3. Assume we are given a finite sequence of rational numbers that comes from a sequence (t n ) n∈{0,1,... } which admits a linear recurrence equation with polynomial coefficients. If the amount of data is sufficiently large, we are able to guess recurrence operators of some fixed order and maximal coefficient degree that annihilate (t n ) n∈{0,1,... } , i.e. the operators applied to the sequence give zero. (For details on guessing and a Mathematica implementation of the method, see Kauers (2009).) For example, consider Given the first 300 terms of this sequence, we can find two operators A and B in É[n][S;s n , 0] with d A = 14, d B = 13 and maximal coefficient degree A = 5, B = 6 resp. Both operators annihilate the given sequence, but none of them is of minimal order. To get an annihilating minimal order operator, we compute the GCRD of A and B in É(n)[S;s n , 0]. Table 1 shows the maximal coefficient degrees of the remainders for different PRSs of A and B.
PRS  The example confirms that the degrees in the pseudo PRS grow exponentially, whereas the subresultant PRS and the primitive PRS show linear growth. At the same time, the degrees in the subresultant PRS are not as small as possible. This behavior is typical not only for this pair A and B, but in general for operators coming from applications. For randomly generated operators, the subresultant PRS and the primitive PRS usually coincide. Our goal is to understand the difference between randomly generated input and the operators A and B as above and to identify the source of some (and most often all) of the additional content in the subresultant PRS. To make use of this knowledge, we will then adjust the formulas for α i and β i from Example 2.3 so that we get a PRS with smaller degrees without having to compute the content of every remainder.

Subresultant Theory for Ore Polynomials
For commutative polynomials, the theory of subresultants was intensively studied by Brown (1978), Brown and Traub (1971), Collins (1967) and Loos (1982). The main idea is to translate relations between the elements of a PRS like the Bézout relation or the (pseudo-)remainder formula into linear algebra. A central tool in this context is the Sylvester matrix, which, roughly speaking, contains the coefficients of all the monomial multiples of the input polynomials that are necessary to compute remainders of any possible degree. The remainders in the subresultant sequence turn out to be polynomials whose coefficients are determinants of certain submatrices of this matrix. Li (1998) generalized these results to Ore polynomials.
lc(x d B −1 A) · · · · · · · · · · · · · · · · · · tc(x d B −1 A) the entry in the ith row and jth column is removing the rows 1 to i, the rows d B + 1 to d B + i, the columns 1 to i and the last i + 1 columns except for the column Theorem 1 (Li (1998)). The polynomial remainder sequence given by α i and β i as in Example 2.3, the subresultant PRS, is equal to the subresultant sequence of A and B of the first kind.

Identifying Content of Polynomial Subresultants
The representation of subresultants in terms of determinants of the matrices Syl i,j (A, B) makes it possible to identify content by exploiting the special form of these matrices as well as the correspondence between rows of the Sylvester matrix and monomial multiples of A and B. For the case of commutative polynomials, some results are known for detecting such additional content. We generalize two results to the Ore setting. The first (Theorem 2) is a generalization of an observation mentioned in Brown (1978), which carries over quite easily to the Ore case. The second (Theorem 4) usually performs better in terms of coefficient size of the remainders, but a heuristic argument is necessary to use it algorithmically (see Section 6).
Proof. Let i be fixed. The coefficients of sres i (A, B) are the determinants of the matrices Syl i,j (A, B) for 0 ≤ j ≤ i. The first column of all of these matrices is Laplace expansion along this column proves the claim.
Not all of the subresultants of A and B are in the subresultant PRS of A and B. To make use of Theorem 2 for a new PRS, we need a minor specialisation of the statement: Corollary 1. Let (R i ) i∈{0,...,ℓ+1} be the subresultant PRS of A and B (not necessarily normal). If we choose Proof. Suppose R i is the jth subresultant of A and B. Then, by the definition of the subresultant sequence of the first kind and Theorem 1, the (j + 1)st subresultant of A and B is regular. Because of this and the subresultant block structure (see Li (1998)), R i−1 is of order j + 1 and so j is equal In the commutative case, a second source of additional content was determined, although this result is not widely known. The following theorem can be found in Knuth (1981): Theorem 3. Let A, B ∈ [x] be such that the subresultant PRS of A and B is normal, i.e. d i−1 = d i + 1 for 1 ≤ i ≤ ℓ, and let G be the GCD of A and B. Then lc(G) 2(i−1) | cont(R i ) for 2 ≤ i ≤ ℓ.
A generalization of Theorem 3 to Ore polynomials is not straightforward, as Example 4 shows.
Example 4 (Example 3 cont.). If we take A and B as in Example 3, then the leading coefficient of the GCRD of A and B is (n + 9)p(n), where p(n) is a polynomial of degree 17. The subresultant PRS of A and B turns out to be normal and R 2 is of order d 2 = 12. By Theorem 3, if the polynomials were elements of [x], cont(R 2 ) would be divisible by lc(G) 2 and a naive translation of the theorem to the non-commutative case suggests divisibility by a polynomial of degree at least 36. The (monic) content of R 2 , however, is only (n + 16)(n + 17), which is contained in, but not equal to, σ 7 (lc(G)) [2] .
Again in the commutative case, let Q A , Q B ∈ [x] be such that A = Q A G and B = Q B G. Knuth (1981) proves Theorem 3 by showing that if (R i ) i∈{0,...,ℓ+1} is the subresultant PRS of A and B and (R i ) i∈{0,...,ℓ+1} is the subresultant PRS of Q A , Q B , then Q i = lc(G) 2(i−1)R i . This approach is problematic for Ore polynomials, because there the Q i 's and theR i 's have coefficients in Ã and not necessarily in . This means that even after showing that a quotient Q i is a -left multiple of some subresultantR i of Q A and Q B , the left factor and the denominators in the coefficients ofR i might not be coprime and thus lead to cancellation. Therefore we will not only describe why in the non-commutative case only some factors of lc(G) appear as content, but we also present a new proof of Theorem 3 that makes it more explicit where the additional content comes from. Moreover, we won't require the remainder sequence to be normal.
In [x], if A is a multiple of the primitive polynomial G, then their quotient will always have coefficients in , and therefore, the leading coefficient of A contains all the factors of the leading coefficient of G. For Ore polynomials, this is not necessarily true, since the quotient of A and G might be an element of Ã[x;σ,δ] \ [x; σ, δ]. Still, different left multiples of G in [x; σ, δ] may share some common factors in their leading coefficients, as described in Lemma 1.
Lemma 1. Let d T ∈ AE be fixed, let I ⊳ [x; σ, δ] be a left ideal and let T be any element of I of order d T such that, among all the operators of order d T in I, its leading coefficient t is minimal with respect to the degree. Then t is independent of the choice of T (up to multiplication by units in ) and for any L Proof. Assume there are T, L ∈ I for which the claim σ d L −d T (t) | lc(L) does not hold. We let L ′ = x d T −d L L and get lc(L ′ ) = σ d T −d L (lc(L)), thus t ∤ lc(L ′ ) by assumption. Division with remainder yields nonzero q, r ∈ such that lc(L ′ ) = qt + r, deg(r) < deg(t).
Hence the operator L ′ − qT is an element of I whose leading coefficient has degree less than deg(t). This contradicts the choice of T . For the uniqueness, let T ′ ∈ I be any other operator of order d T with minimal leading coefficient degree. By what was just shown above, we get lc(T ′ ) | t and t | lc(T ′ ), so t and lc(T ′ ) are associates.
Definition 5. Consider I, T and t from Lemma 1. The shift σ −d T (t) of the leading coefficient of T is called the essential part of I at order d T . If there is no operator in I for some order n, the essential part of I at order n is defined to be 1. is the left ideal generated by L. We give an informal explanation of essential parts of I in terms of solutions of L, i.e. functions that are annihilated by L. Any non-removable singularity of a solution of L corresponds to a root of the leading coefficient of L, but not for any root of lc(L) there has to be a solution with a non-removable singularity at that point. Any solution of L is also a solution of every operator in I and it can happen that there are nonzero Ã-left multiples of L in I that have strictly smaller leading coefficient degree than L. If such a desingularized operator exists, it means that some of the roots of lc(L) can be removed by multiplying L with another operator from the left. These removable roots are called the apparent singularities of L. It is shown in Jaroschek (2013) that there exists a unique minimal (w.r.t. degree) essential part of I that appears in the essential parts of I at every order greater than d L . This minimal essential part of I is a polynomial whose roots are exactly the non-apparent singularities of L, and it turns out that for each root of the essential part of I, there is at least one solution of L that does not admit an analytic continuation at that point. A more detailed description of desingularization and apparent singularities of differential equations can be found in Ince (1926). Further references and recent results on desingularization of Ore operators can be found in Chen et al. (2013).
Note that for commutative polynomials, by Gauß' Lemma, the essential part of a nonzero ideal at any order is equal to the leading coefficient of the primitive greatest common divisor of the ideal elements.
For the remaining part of this article, let I ⊳ [x; σ, δ] be the left ideal generaed by A and B. We formulate our Ore generalization of Theorem 3, where now some of the essential parts of I play the role of the leading coefficient of the GCRD of A and B. (A, B)).
Proof. For any j ∈ {0, . . . , i}, Syl i,j (A, B) is of size ∆ × ∆ and if the last column is removed, the resulting matrix does not depend on j anymore. For n ∈ {1, . . . , ∆ − 1}, let M i,n be the set of all n × n matrices obtained by removing the last ∆ − n columns and any ∆ − n rows from Syl i,j (A, B). The jth coefficient of sres i (A, B) is the determinant of Syl i,j (A, B) and Laplace expansion along the last column shows that it is a -linear combination of the elements of M i,∆−1 . By induction on n we show that the determinant of any element of M i,n is divisible by t ∆+i−n t ∆+i−(n−1) . . . t ∆+i−1 . The theorem is then proven by setting n = ∆ − 1.
For n = 1, the only entry in a matrix in M i,1 is either zero or the leading coefficient of a monomial left multiple of A or B of order ∆ + i − 1, so the claim follows from Lemma 1. Now suppose the claim is true for 1 ≤ n < ∆ − 1 and let M be any element of M i,n+1 . If the determinant of M is zero, then there is nothing to show. Consider the case where det(M ) = 0.
Then there is a v ∈ Ã n+1 such that M T v = (0, . . . , 0, 1) T . By Cramer's rule, the jth component v j of v is of the form p j / det(M ) where p j ∈ is the determinant of some element of M i,n . By induction hypothesis it is divisible by t ∆+i−n t ∆+i−(n−1) . . . t ∆+i−1 . Every row in M corresponds to an operator of the form x k A or x k B for k ∈ AE, minus some of the lower order terms. For the jth row, 1 ≤ j ≤ n + 1, we denote the corresponding operator by L j . By the definition of v, the operator n+1 j=0 v j L j ∈ Ã[x;σ,δ] will have order ∆ + i − (n + 1) and leading coefficient 1. So if we and L = n+1 j=0 v ′ j L j , then L is an element in I of order ∆ + i − (n + 1) and its leading coefficient is det(M )/(t ∆+i−n t ∆+i−(n−1) . . . t ∆+i−1 ) ∈ . Lemma 1 yields that lc(L) is divisible by t ∆+i−(n+1) , so we get in total t ∆+i−(n+1) t ∆+i−n . . . t ∆+i−1 | det(M ).
In practice, the essential parts of I will most likely be the same at every order n with d G ≤ n ≤ d A + d B . In that case, Theorem 4 is equivalent to the following simplification, where only the essential part of I at order d A + d B needs to be known. (A, B)).
Proof. According to Lemma 1, σ j (t) divides the essential part of I at order j for any d G ≤ j ≤ d A + d B . If i < d G , then the ith subresultnat of A and B is zero. Otherwise, Theorem 4 yields that cont(sres i (A, B)) is divisible by Like for Theorem 2, an adjustment of Corollary 2 to the block structure of the subresultant sequence of the first kind is needed in order to construct a new PRS.
Corollary 3. Let (R i ) i∈{0,...,ℓ+1} be the subresultant PRS of A and B (not necessarily normal) and let t be the essential part of I at order Proof. Suppose R i is the jth subresultant of A and B. As in the proof of Corollary 1, we have that j is equal to d i−1 −1. So by Corollary 2, the content of R i is divisible by Simple hand calculation shows that this is equal to γ i .

Improved Polynomial Remainder Sequence
We now derive formulas for the α i and β i that take into account the potential additional content characterized by Theorems 2 and 4. For this we need the following lemma: Proof. By Lemma 2.3 in Li (1998), the pseudo-remainder of γ 1 A and γ 2 B is the (d B − 1)st subresultant of γ 1 A and γ 2 B (up to sign). Consequently, its coefficients are determinants of submatrices of Syl(γ 1 A, γ 2 B) that contain one row corresponding to the operator γ 1 A and d A − d B + 1 rows corresponding to operators of the form x i γ 2 B, 0 ≤ i ≤ d A − d B . Thus, by Lemma 2.2 in Li (1998), it follows that (up to sign) The pseudo-remainder formula (1) applied to γ 1 A and γ 2 B is Combining this with (2) and dividing the resulting equation by γ 1 γ from the left gives the desired result.
This now allows us to state α i and β i for improved polynomial remainder sequences: Theorem 5. Suppose (R i ) i∈{0,...,ℓ+1} is the subresultant PRS of A and B and (γ i ) i∈{0,...,ℓ+1} is any sequence in Ã \ {0} with γ 0 = γ 1 = 1. SetR i = 1 γ i R i . Then (R i ) i∈{0,...,ℓ+1} is a PRS of A and B with:α Proof. From the definition ofR i and the equations For the first summand on the right hand side, Lemma 2 yields For the second summand, observe that since γ i lc(R i ) equals lc(R i ), we have that ψ i equalsψ i for all 1 ≤ i ≤ ℓ. Thus The proof is concluded by combining (3), (4) and (5) and dividing the resulting equation by γ Two possible choices for (γ i ) i∈{i,...,ℓ+1} were presented in Corollary 1 and 3. The computation of γ i in Corollary 1 is straightforward, but in Corollary 3, the essential part of I (the ideal generated by A and B) at order d A + d B is usually not known. A simple heuristic can solve this problem in most cases: As was shown in Lemma 1, the essential part of I at order d A + d B appears in a shifted version in the leading coefficient of every nonzero ideal element with order less than or equal to d A + d B . In particular it is contained in lc(A) and lc(B). Thus, if t is the essential part of I at and in most cases, we not only have divisibility but equality. In fact, in all the examples we looked at that came from combinatorics or physics, this guess for the essential part turned out to be correct.
Example 5 (Example 4 cont.). We now use Theorem 5 and Corollaries 1 and 3 to compute new PRSs of A and B as in Example 3. The essential part of I at order d A + d B is (n + 3), so σ d A (n + 3) = (n + 17), which is also the guess given by the right hand side of (6). Applying Corollary 1 yields the factors γ 2 = n + 17, γ 3 = n + 18, . . . γ i = n + 16 + i − 1, . . .

A, B ∈ Q[y][x],
A = x 4 + yx 2 + yx + y, The second subresultant of A and B is sres 2 (A, B) = (y + y 2 )x 2 + yx + y, so cont(sres 2 (A, B)) = y, but in the improved PRS, no content will be found. As mentioned, it may also happen that the guess for the essential part of I at order d A + d B is too large, for example: Here, cont(R 3 ) in the subresultant PRS is (y + 1), but a factor (y + 1) 2 is predicted. The mistake in predicting the essential part can be noticed on the fly during the execution of the algorithm as soon as a remainder with coefficients in Q(y) appears. It is then possible to either switch to another PRS or to refine the guess of the essential part. One strategy to do so is to remove all the factors from the guess that could be responsible for the appearance of denominators. Let t be the guess for the essential part of I at d A + d B and let c be the non-trivial common denominator of the coefficients of a remainder R i in the improved PRS. Furthermore let M be the set of all integers m such that gcd(σ m (c), t) = 1. Update R i , γ i and t with and continue the computation with these new values. For differential operators in (y)[D; 1, d dy ], we have M = {0} and for recurrence operators in (n)[S n ; s n , 0], M contains all the integer roots of the polynomial res n (c(n + m), t).
The GCRD of A and B is of order 1 and the essential part of I at d A + d B is of degree 4. The essential part of I at order 11, however, is of degree 11, so here we are in the rare case where the essential part of I at order d A + d B is only contained but not equal to the essential part at lower orders. Formula (6)   It is possible to guess the essential part of I at lower orders and then use Theorem 4 to get the primitive remainders, but like in the direct computation of the primitive PRS, GCD computations in the base ring would be necessary after each division step.
As another consequence of Theorem 4, we can give a new bound for the coefficient degrees of the primitive PRS in terms of the essential parts of the left ideal generated by A and B.