Fast matrix computation of subresultant polynomial remainder sequences ^{∗}

Alkiviadis G. Akritas and Gennadi I. Malaschonok ^{†} University of Thessaly (Greece), Tambov University (Russia)

Abstract

We present an improved (faster) variant of the matrix-triangularization subresultant prs method for the computation of a greatest common divisor of two polynomials A and B (of degrees d_{A} and d_{B}, respectively) along with their polynomial remainder sequence [1]. The computing time of our fast method is O(n^{2+β} log ∥C∥^{2}), for standard arithmetic and O(((n^{1+β} +n^{3} log ∥C∥)(log n+ log ∥C∥)^{2}) for the Chinese remainder method, where n = d_{A} + d_{B}, ∥C∥ is the maximal coefficient of the two polynomials and the best known β < 2.356. By comparison, the computing time of the old version is O(n^{5} log ∥C∥^{2}).

Our improvement is based on the work of Malaschonok [11] who proposed a new, recursive method for the solution of systems of linear equations in integral domains with complexity O(n^{β}) over the integers (same as the complexity of matrix multiplication).

In this paper we present an overview of the two methods mentioned above and show how they are combined into a fast matrix method for computing polynomial remainder sequences. An example is also presented to clarify the concepts.

1 Introduction

Let Z be an integral domain, and let

A_{i} = ∑_{j=1}^{m} c_{ij}x^{m−j},

where c_{ij} ∈ Z, i = 1, 2, …, n; then

mat(A_{1}, A_{2}, …, A_{n})

∗This paper was published in the book Computer Algebra in Scientific Computing - CASC 2000, (Ed. by V.G. Ganzha, E.W. Mayr, E.V. Vorozhtsov), Springer, 2000, 1–11. No part of this materials may be reproduced, stored in retrieval system, or transmitted, in any form without prior permission of the copyright owner.

†Akritas is on leave from the University of Kansas, USA

denotes the matrix (c_{ij}) of order n × m. Moreover, let A, B ∈ Z[x], deg A = d_{A}, deg B = d_{B} and let

M_{k} = mat(x^{d_{B}−k−1}A, x^{d_{B}−k−2}A, …, A, x^{d_{A}−k−1}B, x^{d_{A}−k−2}B, …, B),
0 ≤ k < min(d_{A}, d_{B})

be the matrix of order (d_{A} + d_{B} − 2k) × (d_{A} + d_{B} − k), where M_{0} is the well-known Sylvester’s matrix. Then, kth subresultant polynomial of A and B is called the polynomial

S_{k} = ∑_{i=0}^{k} M_{k}^{i}x^{i},

of degree ≤ k, where M_{k}^{i} is a minor of the matrix M_{k} of order d_{A} + d_{B} − 2k, formed by the elements of columns 1, 2, …, d_{A} + d_{B} − 2k − 1 and column d_{A} + d_{B} − k − i. Habicht’s known theorem [8] establishes a relation between the subresultant polynomials S_{0}, S_{1}, …, S_{min(d_{A},d_{B})−1} and the polynomial remainder sequence(prs) of A and B, and also demonstrates the so-called gap structure.

According to the matrix-triangularization subresultant prs method (see for example Akritas’ book [2] or papers [3], [4]) all the subresultant polynomials of A and B can be computed within sign by transforming the matrix (suggested by Sylvester [13])

mat(x^{max(d_{A},d_{B})−1}A, x^{max(d_{A},d_{B})−1}B, x^{max(d_{A},d_{B})−2}A, x^{max(d_{A},d_{B})−2}B, …, A, B),

of order 2 · max(d_{A}, d_{B}), into its upper triangular form with the help of Dodgson’s integer preserving transformations [7]; they are then located using an extension of a theorem by Van Vleck [3], [14]. (We depart from established practice and we give credit to Dodgson, and not to Bareiss [5], for the integer preserving transformations; see also the work of Waugh and Dwyer [15] where they use the same method as Bareiss, but 23 years earlier, and they name Dodgson as their source–differing from him only in the choice of the pivot element ([15], page 266). Charles Lutwidge Dodgson (1832–1898) is the same person widely known for his writing Alice in Wonderland under the pseudonym Lewis Carroll.)

In Section 2 we present an overview of the matrix-triangularization subresultant prs method allowing us to exactly compute and locate the members of the prs (without using Van Vleck’s theorem [14]) by applying Dodgson’s integer preserving transformations to a matrix of order d_{A} + d_{B}.

In Section 3 we present an overview of the recursive method for the solution of systems of linear equations in integral domains. Its complexity is the same as that of matrix multiplication.

In Section 4 we prove that the recursive method allowing us to compute and locate the members of the prs and gcd.

In Section 5 we present an example to clarify the concepts.

2 The ”slow” matrix computation of subresultant polynomial remainder sequences

We assume that deg A = d_{A} ≥ deg B = d_{B} and we denote by M the following matrix

M = mat(x^{d_{A}−1}B, x^{d_{A}−2}B, …, x^{d_{B}−1}B,^{d_{B}−1} A, x^{d_{B}−2}B, x^{d_{B}−2}A, …, B, A)

of order d_{A} + d_{B} with elements a_{ij}(j, i = 1, 2, …, d_{A} + d_{B}). (This matrix can be obtained from Sylvester’s matrix M_{0} after a rearrangement of its rows.)

Dodgson’s integer preserving transformations

a_{ij}^{k+1} = \frac{(a_{ij}^{k}a_{kk}^{k}−a_{ik}^{k}a_{kj}^{k})}{a_{k−1,k−1}^{k−1}}  (D)

(see [5], [7], [9] or [15]) where we set a_{00}^{0} = 1, a_{ij}^{1} = a_{ij} and it is assumed that a_{kk}^{k} ≠ 0, k = 1, 2, …, d_{A} + d_{B}, are applied to the matrix M = (a_{ij}) and transform it to the upper-triangular matrix M_{D} = (b_{ij}), (i, j = 1, 2, …, d_{A} + d_{B}), where

b_{ij} = \{ \begin{matrix}0 & for i > j \\ a_{ij}^{i} & for i ≤ j \end{matrix}

and, in general,

a_{ij}^{k} = | \begin{matrix}a_{11} & … & a_{1,k−1} & a_{1j} \\ ⋮ & ... & ⋮ & ⋮ \\ a_{k−1,1} & … & a_{k−1,k−1} & a_{k−1,j} \\ a_{i1} & … & a_{i,k−1} & a_{ij} \end{matrix} |

with 1 ≤ k ≤ d_{A} + d_{B}, and k ≤ i, j ≤ d_{A} + d_{B}.

The upper-triangular matrix M_{D} = (b_{ij}) will be called Dodgson matrix.

The following two theorems can be used to locate the members of the prs in the rows of M_{D}. The correct sign is computed.

Case 1: If none of the diagonal minors of the matrix M is equal to zero, then we have:

Theorem 1. Dodgson’s integer preserving transformation will transform matrix M to the upper triangular matrix M_{D}, which contains all n subresultants (located in rows d_{A} + d_{B} − 2k, k = 0, 1, 2, …, d_{B} − 1) S_{k} = ∑_{i=0}^{k} M_{k}^{i}x^{i},

M_{k}^{i} = (−1)^{σ(k)}a_{d_{A}+d_{B}−2k,d_{A}+d_{B}−k−i}^{d_{A}+d_{B}−2k}

and σ(k) = (d_{A} − d_{B} + 1) + ⋯ + (d_{A} − k) = (1/2)(d_{B} − k)(2d_{A} − d_{B} − k + 1), k = 0, 1, …, d_{B} − 1.

Case 2: If not all diagonal minors of the matrix M are nonzero, then we have the following theorem (the term bubble pivot, used below, means that, after pivoting, row i_{p} is immediately below row j_{p}):

Theorem 2. Dodgson’s integer preserving transformations with bubble pivot and choice of the pivot element by column, will transform matrix M to the upper triangular matrix M_{D}, and at the same time will compute all subresultants S_{k}; if, in the process, s row replacements take place, namely row j_{1} replaces row i_{1}, j_{2} replaces i_{2}, …, j_{s} replaces i_{s}, (and after each replacement row i_{p} is immediately below row j_{p}, p = 1, 2, …, s), then (a) S_{k} = 0, for all k such that \frac{(d_{A}+d_{B}−i_{p})}{2} > k > \frac{(d_{A}+d_{B}−j_{p})}{2} and for all p = 1, 2, …, s. (b) for all p = 1, 2, …, s, if k = \frac{(d_{A}+d_{B}−i_{p})}{2} is an integer number not in (a), S_{k} is located in row i_{p} before it is replaced by row j_{p}. (c) for the remaining k, (k = 0, 1, …, d_{B} − 1 and those not in (a) or (b)) S_{k} is located in row j = d_{A} + d_{B} − 2k.

Moreover, in (b) and (c) the subresultant S_{k} = ∑_{i=0}^{k} M_{k}^{i}x^{i}, is located in row j in such a way that

M_{k}^{i} = (−1)^{σ(k)+σ(j)}a_{j,j+k−i}^{j}

σ(k) = (1/2)(d_{B} − k)(2d_{A} − d_{B} − k + 1), σ(j) = ∑_{p=1}^{s} j_{p} − ∑_{p=1}^{s} i_{p}, j_{p} ≤ j, i_{p} ≤ j.

Note that in cases (b) and (c) Theorem 2 reduces to Theorem 1 in the case of a complete prs, and due to the fact that rows above row j change places, the sign changes by a factor (−1)^{σ(j)}.

For the given polynomials over the integers in this discussion, the complexity of this method is

O(n^{5} log ∥C∥^{2}))

where n = d_{A} + d_{B} and ∥C∥ is the maximal coefficient of the two polynomials.

3 The recursive method

Let ∑_{j=1}^{m−1} a_{ij}x_{j} = a_{im}, i = 1, 2, …, n be the system of linear equations with extended coefficients matrix

A = (a_{ij}), i = 1, …, n, j = 1, …, m.

whose coefficients are in integral domain Z: A ∈ Z^{n×m}.

The solution of such a system may be written according to Cramer’s rule

x_{j} = \frac{δ_{jm}^{n} − ∑_{p=n+1}^{m−1} x_{p}δ_{jp}^{n}}{δ^{n}} , j = 1, …, n,

where x_{p}, p = n + 1, …, m, are free variables and δ^{n} ≠ 0. δ^{n} = |a_{ij}|, i = 1, …, n, j = 1, …, n, - denote the corner minors of the matrix A of order n, δ_{ij}^{n} - denote the minors obtained by a substitution of the column j of the matrix A instead of the column i in the minors δ^{n}, i = 1, …, n, j = n + 1, …, m. So we need to construct an algorithm for computing the minor δ^{n} and the matrix G = ( δ_{ij}^{n}), i = 1, …, n, j = n + 1, n + 2, …, m.

This means that we must make the reduction of the matrix A to the diagonal form

A → (δ^{n}I_{n}, G).

I_{n} denotes the unit matrix of order n.

For the extended coefficients matrix A we use the following notation:

A_{ij}^{k} = ( \begin{matrix}a_{11} & a_{12} & ⋯ & a_{1,k−1} & a_{1j} \\ a_{21} & a_{22} & ⋯ & a_{2,k−1} & a_{2j} \\ ⋮ & ⋮ & ... & ⋮ & ⋮ \\ a_{k−1,1} & a_{k−1,2} & ⋯ & a_{k−1,k−1} & a_{k−1,j} \\ a_{i1} & a_{i2} & ⋯ & a_{i,k−1} & a_{ij} \end{matrix} )

is the matrix, formed by surrounding the sub-matrix of order k − 1 in the upper left corner by row i and column j,

a_{ij}^{k} = det A_{ij}^{k},

a_{ij}^{1} = a_{ij}, δ^{0} = 1, δ^{k} = a_{kk}^{k}, δ_{ij}^{k} is the determinant of the matrix, that is obtained from the matrix A_{kk}^{k} after the substitution of column i by column j.

We shall use the minors δ_{ij}^{k} and a_{ij}^{k} for the construction of the matrices

A_{k,c}^{r,l,(p)} = ( \begin{matrix}a_{r+1,k+1}^{p} & a_{r+1,k+2}^{p} & ⋯ & a_{r+1,c}^{p} \\ a_{r+2,k+1}^{p} & a_{r+2,k+2}^{p} & ⋯ & a_{r+2,c}^{p} \\ ⋮ & ⋮ & ... & ⋮ \\ a_{l,k+1}^{p} & a_{l,k+2}^{p} & ⋯ & a_{l,c}^{p} \end{matrix} )

and

G_{k,c}^{r,l,(p)} = ( \begin{matrix}δ_{r+1,k+1}^{p} & δ_{r+1,k+2}^{p} & ⋯ & δ_{r+1,c}^{p} \\ δ_{r+2,k+1}^{p} & δ_{r+2,k+2}^{p} & ⋯ & δ_{r+2,c}^{p} \\ ⋮ & ⋮ & ... & ⋮ \\ δ_{l,k+1}^{p} & δ_{l,k+2}^{p} & ⋯ & δ_{l,c}^{p} \end{matrix} )

G_{k,c}^{r,l,(p)}, A_{k,c}^{r,l,(p)} ∈ Z^{(l−r)×(c−k)}, 0 ≤ k < n, k < c ≤ n, 0 ≤ r < m, r < l ≤ m, 1 ≤ p ≤ n.

We describe one recursive step, that makes the following reduction of the matrix Ã to the diagonal form

Ã → (δ^{l}I_{l−k}, G_{l,c}^{k,l,(l)}), 0 ≤ k < c ≤ m, k < l ≤ n, l < c.

Note that if k = 0, l = n and c = m, then we obtain the solution of the system.

We choose the arbitrary integer number s: k < s < l and write the matrix Ã = A_{k,c}^{s,l,(k+1)} as follows:

Ã = ( \begin{matrix}A_{k,c}^{k,s(k+1)} \\ A_{k,c}^{s,l,(k+1)} \end{matrix} )

where A_{k,c}^{k,s,(k+1)} is the upper part of the matrix Ã consisting of s − k rows and A_{k,c}^{s,l,(k+1)} is the lower part of the matrix Ã.

The first step

We make the following reduction of the matrix A_{k,c}^{k,s(k+1)} ∈ Z^{(s−k)×(c−k)} to the diagonal form

A_{k,c}^{k,s(k+1)} → (δ^{s}I_{s−k}, G_{s,c}^{k,s,(s)}).

The second step We write the matrix A_{k,c}^{s,l,(k+1)} in the following way:

A_{k,c}^{s,l,(k+1)} = (A_{k,s}^{s,l,(k+1)}, A_{s,c}^{s,l,(s+1)})

where A_{k,s}^{s,l,(k+1)} consists of the first s − k columns and A_{s,c}^{s,l,(k+1)} consists of the last c − s columns of the matrix A_{k,c}^{s,l,(k+1)}.

The matrix A_{s,c}^{s,l,(s+1)} is obtained from the matrix identity:

δ^{k} · A_{s,c}^{s,l,(s+1)} = δ^{s} · A_{s,c}^{s,l,(k+1)} − A_{k,s}^{s,l,(k+1)} · G_{s,c}^{k,s,(s)}.

The minors δ^{k} must not equal zero.

The third step As the next recursive step we make the following reduction of the matrix A_{s,c}^{s,l,(s+1)} ∈

Z^{(l−s)×(c−s)} to the diagonal form

A_{s,c}^{s,l,(s+1)} → (δ^{l}I_{l−s}, G_{l,c}^{s,l,(l)}).

The fourth step We write the matrix G_{s,c}^{k,s,(s)} in the following way:

G_{s,c}^{k,s,(s)} = (G_{s,l}^{k,s,(s)}, G_{l,c}^{k,s,(s)})

where G_{s,l}^{k,s,(s)} consists of the first l − s columns and G_{l,c}^{k,s,(s)} consists of the last c − l columns of the matrix G_{s,c}^{k,s,(s)}.

The matrix G_{l,c}^{k,s,(l)} is obtained from the matrix identity:

δ^{s} · G_{l,c}^{k,s,(l)} = δ^{l} · G_{l,c}^{k,s,(s)} − G_{s,l}^{k,s,(s)} · G_{l,c}^{s,l,(l)}.

The minors δ^{s} must not equal zero.

So we get

G_{l,c}^{k,l,(l)} = ( \begin{matrix}G_{l,c}^{k,s,(l)} \\ G_{l,c}^{s,l,(l)} \end{matrix} )

and δ^{l}.

3.1 Representation of the one recursive step

We can represent one recursive step as the following reduction of the matrix A_{kc}^{kl(k+1)}:

A_{kc}^{kl(k+1)} = ( \begin{matrix}A_{k,c}^{k,s(k+1)} \\ A_{k,c}^{s,l,(k+1)} \end{matrix} ) →_{1} ( \begin{matrix}δ^{s}I_{s−k} & G_{sc}^{ks(s)} \\ A_{k,s}^{s,l,(k+1)} & A_{s,c}^{s,l,(k+1)} \end{matrix} ) →_{2} ( \begin{matrix}δ^{s}I_{s−k} & G_{sc}^{ks(s)} \\ 0 & A_{sc}^{sl(s+1)} \end{matrix} ) →_{3}

→_{3} ( \begin{matrix}δ^{s}I_{s−k} & G_{sl}^{ks(s)} & G_{lc}^{ks(s)} \\ 0 & δ^{l}I_{l−s} & G_{lc}^{sl(l)} \end{matrix} ) →_{4} ( \begin{matrix}δ^{s}I_{s−k} & 0 & G_{lc}^{ks(l)} \\ 0 & δ^{l}I_{l−s} & G_{lc}^{sl(l)} \end{matrix} ) → ( δ^{l}I_{s−k} G_{lc}^{kl(l)} )

4 Dichotomous recursion

If the process of partition is dichotomous, and the number of rows in the upper and lower submatrixes is the same and equal to a power of 2 in every step, then we will call such process - the dichotomous recursion process.

Let us denote

A_{k}^{2^{p}} = A_{kc}^{k,k+2^{p}(k+1)}, G_{k}^{2^{p}} = G_{kc}^{k−2^{p},k(k)}

matrices with 2^{p} rows and c−k columns. For any natural number s we denote by b[s] the maximal power of 2, that devide the nubmer s, i.e. for s = ∑_{i=0}^{j} c_{i}2^{i}, c_{i} ∈ {0, 1} b[s] = 2^{w}, where w = min{i|c_{i} ≠ 0}.

Then one dichotomous recursive step is as follows:

A_{t2^{p}}^{2^{p}} →_{1} ( \begin{matrix}δ^{(2t+1)2^{p−1}}I_{2^{p−1}} & G_{(2t+1)2^{p−1}}^{2^{p−1}} \\ ∗ & ∗ \end{matrix} ) →_{2} ( \begin{matrix}δ^{(2t+1)2^{p−1}}I_{2^{p−1}} & G_{(2t+1)2^{p−1}}^{2^{p−1}} \\ 0 & A_{(2t+1)2^{p−1}}^{2^{p−1}} \end{matrix} )

→_{3} ( \begin{matrix}∗ & ∗ & ∗ \\ 0 & δ^{(t+1)2^{p}}I_{2^{p−1}} & G_{(t+1)2^{p}}^{2^{p−1}} \end{matrix} ) →_{4} ( δ^{(t+1)2^{p}}I_{2^{p}} G_{(t+1)2^{p}}^{2^{p}})

p = 1, 2, …, N, t = 0, 1, …, 2^{N−p} − 1, n = 2^{N}.

We obtain the following sequense of matrices:

A_{0}^{2}; G_{2}^{2}, A_{2}^{2}; G_{4}^{2}, G_{4}^{4}, A_{4}^{4}; G_{6}^{2}, A_{6}^{2}, G_{8}^{2}, G_{8}^{4}, G_{8}^{8}, A_{8}^{8}; G_{10}^{2}, A_{10}^{2}, G_{12}^{2}, G_{12}^{4}, A_{12}^{4}, G_{14}^{2}, A_{14}^{2}, G_{16}^{2}, G_{16}^{4}, G_{16}^{8}, G_{16}^{16}, … A_{4s}^{b[4s]}, G_{4s+2}^{2}, A_{4s+2}^{2}, G_{4s+4}^{2}, G_{4s+4}^{4}, …, G_{4s+4}^{b[4s+4]}, … which we call the dichotomous sequence of matrices.

Theorem 3. The dichotomous sequence of matrices contains all the rows of the Dodgson matrix M_{D}. Proof: The dichotomous sequence of matrices includes the matrices A_{t}^{b[t]}, t = 2, 4, …, n − 2, G_{t}^{2}, t = 2, 4, …, n.

The matrix A_{t}^{b[t]}, t = 2, 4, …, n − 2, has the next first row

(a_{t+1,t+1}^{t+1}, a_{t+1,t+2}^{t+1}, …, a_{t+1,m}^{t+1}),

the matrix G_{t}^{2}, t = 2, 4, …, n − 2, has the next last row

(δ_{t,t}^{t}, δ_{t,t+1}^{t}, …, δ_{t,m}^{t}) = (a_{t,t}^{t+1}, a_{t,t+1}^{t+1}, …, a_{t,m}^{t+1})

.

There are rows of the Dodgson matrix M_{D} with the numbers t+1 (t = 2, 4, …, n− 2) and t (t = 2, 4, …, n) respectively. So the dichotomous sequence of matrices contains all the rows of the Dodgson matrix M_{D} . □

Therefore, we can use the dichotomous sequence of matrices for computing the Dodgson matrix with the prs’s and gcd.

Comparing with Dodgson’s method, care must be taken not to make bubblepivot in each row with zero diagonal element but just in the even rows, because just the diagonal elements in the even rows will be the pivots. Details can be seen in the example that follows – taken from our previous paper [1].

5 Complexity issue

Let ∥C∥ denote the maximal coefficient of the two polynomials, n^{β} denote the complexity of matrix multiplication of order n, n = d_{A} + d_{B}, where d_{A} and d_{B} are the degrees of the two polynomials. The best known β < 2.356. Since the fast method of matrix transformation required n^{β} operations with integer numbers so it is easy to obtain the complexity of our method.

Our method is required

O(n^{2+β} log^{2} ∥C∥)

bit operations for standard arithmetic and

O(((n^{1+β} + n^{3} log ∥C∥)(log n + log ∥C∥)^{2})

for the Chinese remainder method. The cost of reducing modulo a prime is included.

6 Example of Fast Matrix Transformation

We present an example to clarify the concepts discussed above. It is the same example taken from our previous paper [1].

For the polynomials

A = 2x^{4} + 5x^{3} + 5x^{2} − 2x + 1
B = 3x^{3} + 3x^{2} + 3x − 4

we obtain the following matrices:

( \begin{matrix}3 & 3 & 3 & −4 & 0 & 0 & 0 \\ 2 & 5 & 5 & −2 & 1 & 0 & 0 \\ 0 & 3 & 3 & 3 & −4 & 0 & 0 \\ 0 & 2 & 5 & 5 & −2 & 1 & 0 \\ 0 & 0 & 3 & 3 & 3 & −4 & 0 \\ 0 & 0 & 2 & 5 & 5 & −2 & 1 \\ 0 & 0 & 0 & 3 & 3 & 3 & −4 \end{matrix} ) →_{1} ( \begin{matrix}9 & 0 & 0 & −14 & −3 & 0 & 0 \\ 0 & 9 & 9 & 2 & 3 & 0 & 0 \\ 0 & 3 & 3 & 3 & −4 & 0 & 0 \\ 0 & 2 & 5 & 5 & −2 & 1 & 0 \\ 0 & 0 & 3 & 3 & 3 & −4 & 0 \\ 0 & 0 & 2 & 5 & 5 & −2 & 1 \\ 0 & 0 & 0 & 3 & 3 & 3 & −4 \end{matrix} ) →_{2}

( \begin{matrix}9 & 0 & 0 & −14 & −3 & 0 & 0 \\ 0 & 9 & 9 & 2 & 3 & 0 & 0 \\ 0 & 0 & 0 & 21 & 45 & 0 & 0 \\ 0 & 0 & 27 & 41 & −24 & 9 & 0 \\ 0 & 0 & 3 & 3 & 3 & −4 & 0 \\ 0 & 0 & 2 & 5 & 5 & −2 & 1 \\ 0 & 0 & 0 & 3 & 3 & 3 & −4 \end{matrix} ) →_{3} ( \begin{matrix}9 & 0 & 0 & −14 & −3 & 0 & 0 \\ 0 & 9 & 9 & 2 & 3 & 0 & 0 \\ 0 & 0 & −63 & 0 & −149 & −21 & 0 \\ 0 & 0 & 0 & −63 & 135 & 0 & 0 \\ 0 & 0 & 3 & 3 & 3 & −4 & 0 \\ 0 & 0 & 2 & 5 & 5 & −2 & 1 \\ 0 & 0 & 0 & 3 & 3 & 3 & −4 \end{matrix} )

→_{4} ( \begin{matrix}−63 & 0 & 0 & 0 & 231 & 0 & 0 \\ 0 & −63 & 0 & 0 & 98 & 21 & 0 \\ 0 & 0 & −63 & 0 & −149 & −21 & 0 \\ 0 & 0 & 0 & −63 & 135 & 0 & 0 \\ 0 & 0 & 3 & 3 & 3 & −4 & 0 \\ 0 & 0 & 2 & 5 & 5 & −2 & 1 \\ 0 & 0 & 0 & 3 & 3 & 3 & −4 \end{matrix} ) →_{5}

( \begin{matrix}−63 & 0 & 0 & 0 & 231 & 0 & 0 \\ 0 & −63 & 0 & 0 & 98 & 21 & 0 \\ 0 & 0 & −63 & 0 & −149 & −21 & 0 \\ 0 & 0 & 0 & −63 & 135 & 0 & 0 \\ 0 & 0 & 0 & 0 & −147 & 315 & 0 \\ 0 & 0 & 0 & 0 & −692 & 168 & −63 \\ 0 & 0 & 0 & 0 & −594 & −189 & 252 \end{matrix} ) →_{6}

( \begin{matrix}−63 & 0 & 0 & 0 & 231 & 0 & 0 \\ 0 & −63 & 0 & 0 & 98 & 21 & 0 \\ 0 & 0 & −63 & 0 & −149 & −21 & 0 \\ 0 & 0 & 0 & −63 & 135 & 0 & 0 \\ 0 & 0 & 0 & 0 & −3068 & 0 & −315 \\ 0 & 0 & 0 & 0 & 0 & −3068 & −147 \\ 0 & 0 & 0 & 0 & −594 & −189 & 252 \end{matrix} ) →_{7}

( \begin{matrix}−63 & 0 & 0 & 0 & 231 & 0 & 0 \\ 0 & −63 & 0 & 0 & 98 & 21 & 0 \\ 0 & 0 & −63 & 0 & −149 & −21 & 0 \\ 0 & 0 & 0 & −63 & 135 & 0 & 0 \\ 0 & 0 & 0 & 0 & −3068 & 0 & −315 \\ 0 & 0 & 0 & 0 & 0 & −3068 & −147 \\ 0 & 0 & 0 & 0 & 0 & 0 & 15683 \end{matrix} )

The obtained polynomial remainder sequence is incomplete and we only have the remainder −63x + 135 of degree 1 after the thired step and remainder 15683 of degree 0 after the seventh step.

Matrix operations in the each step:

1. δ_{2} = | \begin{matrix}3 & 3 \\ 2 & 5 \end{matrix} | = 9; | \begin{matrix}3 & 3 \\ 2 & 5 \end{matrix} | = 9; | \begin{matrix}3 & −4 \\ 2 & −2 \end{matrix} | = 2; | \begin{matrix}3 & 0 \\ 2 & 1 \end{matrix} | = 3; | \begin{matrix}3 & 0 \\ 2 & 0 \end{matrix} | = 0;

| \begin{matrix}3 & 0 \\ 2 & 0 \end{matrix} | = 0; | \begin{matrix}3 & 3 \\ 5 & 5 \end{matrix} | = 0; | \begin{matrix}−4 & 3 \\ −2 & 5 \end{matrix} | = −14; | \begin{matrix}0 & 3 \\ 1 & 5 \end{matrix} | = −3; | \begin{matrix}0 & 3 \\ 0 & 5 \end{matrix} | = 0; | \begin{matrix}0 & 3 \\ 0 & 5 \end{matrix} | = 0.

2. 9 ( \begin{matrix}3 & 3 & −4 & 0 & 0 \\ 5 & 5 & −2 & 1 & 0 \end{matrix} ) − ( \begin{matrix}0 & 3 \\ 0 & 2 \end{matrix} ) ( \begin{matrix}0 & −14 & −3 & 0 & 0 \\ 9 & 2 & 3 & 0 & 0 \end{matrix} ) =

( \begin{matrix}0 & 21 & 45 & 0 & 0 \\ 27 & 41 & −24 & 9 & 0 \end{matrix} ) .

3. δ_{4} = | \begin{matrix}0 & 21 \\ 27 & 41 \end{matrix} | \frac{1}{9} = −63; | \begin{matrix}0 & −45 \\ 27 & −24 \end{matrix} | \frac{1}{9} = 135; | \begin{matrix}0 & 0 \\ 27 & 9 \end{matrix} | \frac{1}{9} = 0; | \begin{matrix}0 & 0 \\ 27 & 0 \end{matrix} | \frac{1}{9} = 0;

| \begin{matrix}−45 & 21 \\ −24 & 41 \end{matrix} | \frac{1}{9} = −149; | \begin{matrix}0 & 21 \\ 9 & 41 \end{matrix} | \frac{1}{9} = −21; | \begin{matrix}0 & 21 \\ 0 & 41 \end{matrix} | \frac{1}{9} = 0.

4. ((−63) ( \begin{matrix}−3 & 0 & 0 \\ 3 & 0 & 0 \end{matrix} ) − ( \begin{matrix}0 & −14 \\ 9 & 2 \end{matrix} ) ( \begin{matrix}−149 & −21 & 0 \\ 135 & 0 & 0 \end{matrix} )) \frac{1}{9} = ( \begin{matrix}231 & 0 & 0 \\ 98 & 21 & 0 \end{matrix} ) .

5. (−63) ( \begin{matrix}3 & −4 & 0 \\ 5 & −2 & 1 \\ 3 & 3 & −4 \end{matrix} ) − ( \begin{matrix}0 & 0 & 3 & 3 \\ 0 & 0 & 2 & 5 \\ 0 & 0 & 0 & 3 \end{matrix} ) ( \begin{matrix}231 & 0 & 0 \\ 98 & 21 & 0 \\ −149 & −21 & 0 \\ 135 & 0 & 0 \end{matrix} ) =

( \begin{matrix}−147 & 315 & 0 \\ −692 & 168 & −63 \\ −594 & −189 & 252 \end{matrix} ) .

6. δ_{6} = \frac{| \begin{matrix}−147 & 315 \\ −692 & 168 \end{matrix} |}{(−63)} = −3068; \frac{| \begin{matrix}−147 & 0 \\ −692 & −63 \end{matrix} |}{(−63)} = −147; \frac{| \begin{matrix}0 & 315 \\ −63 & 168 \end{matrix} |}{(−63)} = −315;

7. δ_{7} = (252(−3068) − (−315)(−594) − (147)(−189))/(−63) = 15683.

References

[1] Akritas, A.G., Akritas, E.K. and Gennadi I. Malaschonok. Matrix computation of subresultant polynomial remainder sequences in integral domains. Reliable Computing 1 (1995), pp. 375–381.

[2] Akritas, A. G. Elements of Computer Algebra with Applications. J. Wiley Interscience, New York, 1989.

[3] Akritas, A. G. A new method for computing polynomial greatest common divisors and polynomial remainder sequences. Numerische Mathematik 52, 119– 127, 1988.

[4] Akritas, A. G. Exact algorithms for the matrix-triangularization subresultant prs method. Proceedings of the Conference on Computers and Mathematics, Boston, Massachusetts, 145–155, June 1989.

[5] Bareiss, E. H. Sylvester’s identity and multistep integer-preserving Gaussian elimination. Mathematics of Computation 22, 565–578, 1968.

[6] Coppersmith, D. and S. Winograd. in Proc. 19th Annu ACM Symp. on Theory of Comput., 1987, 1–6.

[7] Dodgson, C. L. Condensation of determinants. Proceedings of the Royal Society of London 15. 150–155, 1866.

[8] Habicht, W. Eine Verallgemeinerung des Sturmschen Wurzelzaelverfahrens. Commentarii Mathematici Helvetici 21, 99–116, 1948.

[9] Malaschonok, G. I. Solution of a system of linear equations in an integral domain. USSR Journal of Computational Mathematics and Mathematical Physics 23, 1497–1500, 1983 (in Russian).

[10] G.I. Malaschonok. Algorithms for the solution of systems of linear equations in commutative rings. In Effective Methods in Algebraic Geometry, Edited by T. Mora and C. Traverso, Progress in Mathematics 94, Birkhauser, Boston-Basel-Berlin, 1991, 289–298.

[11] G. I. Malaschonok. Recursive Method for the Solution of systems of Linear Equations,Computational Mathematics (A. Sydow Ed, Proceedings of the 15th IMACS World Congress, Vol. I, Berlin, August 1997 ), Wissenschaft & Technik Verlag, Berlin 1997, 475–480.

[12] Strassen, V. Gaussian Elimination is not optimal. Numerische Mathematik 13, 354–356, 1969.

[13] Sylvester, J. J. On a theory of the syzygetic relations of two rational integral functions, comprising an application to the theory of Sturm’s functions, and that of the greatest common measure. Philoshophical Transactions 143, 407– 548, 1853.

[14] Van Vleck, E. B. On the determination of a series of Sturm’s functions by the calculation of a single determinant. Annals of Mathematics 1, Second Series, 1–13, 1899–1900.

[15] Waugh, F. V. and P. S. Dwyer Compact computation of the inverse of a matrix. Annals of Mathematical Statistics 16, 259–271, 1945.