Polynomial Matrices

Author

Adriano Pascoletti
Dipartimento di Matematica e Informatica
Università di Udine

v. 2.1. december 1999

Summary

This package provides some useful and efficient functions for treating polynomial matrices (i.e. matrices whose entries are univariate polynomials with rational or symbolic coefficients). The supplied functions compute the classical Smith, Hermite and McMillan forms. Each function has two versions: one returning only the form and one returning 
the form and the unimodular transformations leading to it. 
Several other functions are defined: computation of left and right GCD’s,  lcm’s, quotients, remainders, row/column proper forms, full rank and coprimality tests,  extended polynomial GCD, solution of diophantine equations. 
All functions work on matrices of polynomials in z or (z^(-1)).

Notebook version

2.1

Mathematica Version

3.x and 4.0

History

This package supersedes package 0207-751. 
Version 2.1 removes some incompatibilities between Mathematica 2.2  and Mathematica 3.x and 4.0.
Version 2.0 is the result of a work with M. Anziutti.

Keywords

Polynomials, polynomial matrices, non commutative algebra, Smith form, Hermite Form, McMillan Form, diophantine equations, coprimality, properness, greatest common divisors, least common multiples.

Interface

Package context

BeginPackage["PolynomialMatrices`"]

Usage messages

DiophantineSolve::usage = "DiophantineSolve[a,b,c,x] takes three \
matrices a,b,c of polynomials in the indeterminate x with the same number of \
rows and returns a minimum column-degree polynomial solution {X,Y} of the \
equation a.X+b.Y==c. When no solution exists the result is { {{}},{{}} }. In \
order to get all the minimum column degree solutions {X,Y}, the user must \
specify a name for the free variables, e.g. t in DiophantineSolve[a,b,c,x,t].\
";

EliminateDep::usage = "EliminateDep[m,x] eliminates from m (matrix \
of polynomials in x) the dependent rows and returns a matrix with full row \
rank.";

ExtendedHermiteForm::usage = "HermiteForm[a,x] yields the Hermite \
form h of a matrix a of polynomials in x. The result is a list {h,u} s.t. \
h=u.a with u unimodular.";

ExtendedMcMillanForm::usage = "ExtendedMcMillanForm[a,x] returns a \
list {m,u,v} s.t. m=u.a.v. m is the McMillan form of a matrix a of rational \
functions in x, and u,v are unimodular matrives.";

ExtendedSmithForm::usage = "ExtendedSmithForm[a,x] yields the Smith \
form s of the matrix a of polynomials in x, together with two unimodular \
matrices u, v such that s=u.a.v. The result is a list {s,u,v}.";

ExtPolQ::usage = "ExtPolQ[l,x] yields True if all the components of \
the list l are polynomials in x or all in x^-1.";

FullRowRankQ::usage = "FullRowRankQ[a,x] yields True if a \
(polynomial matrix in x) has full row rank.";

HermiteForm::usage = "HermiteForm[m,x] yields the Hermite form of a \
matrix m of polynomials in x.";

LCoprime::usage = "LCoprime[a,b,x], (a,b polynomial matrices in x \
with the same number of rows) returns True if the two matrices are left \
coprime.";

LDivision::usage = "LDivision[a,b,x] takes two matrices of \
polynomials in x with the same number of rows, b square and regular, and \
returns a list {q,r} s.t. a=b.q+r with deg[r]<deg[b].";

LGCD::usage = "LGCD[a,b,x] gives a left greatest common divisor of \
the matrices a and b of polynomials in x  (a and b must have the same number \
of rows).";

Llcm::usage = "Llcm[a,b,x] gives a left least common multiple of \
the matrices a and b whose elements are polynomials in x (a and b must have \
the same number of columns).";

LQuotient::usage = "LQuotient[a,b,x] gives the left quotient of the \
matrices a (l x m) and b (l x l proper) of polynomials in x. The result is a \
matrix q such that a=b.q+r and deg[r]<deg[b].";

LRemainder::usage = "LRemainder[a,b,x] returns the left remainder \
from dividing the matrices a (l x m) and b (l x l proper) of polynomials in \
x. The result is a matrix r such that a=b.q+r and deg[r]<deg[b].";

McMillanForm::usage = " McMillanForm[a,x] takes a matrix of \
rational forms in x and returns its  McMillan form.";

PolyExtendedGCD::usage = "PolyExtendedGCD[{a,b},x] takes two \
polynomials in x (a,b not both null) and returns a list {g,{r,s}}, where g is \
a monic GCD of a and b, and r,s s.t. g=ar+bs.";

PolynomialColumnReduce::usage = "PolynomialColumnReduce[m,x] takes \
a polynomial matrix in x with full column rank and returns a unimodular \
matrix u s.t. r=m.u is a column reduced polynomial matrix.";

PolynomialDegree::usage = "PolynomialDegree[m,x] returns the degree \
of a matrix m of polynomials in x.";

PolynomialRowReduce::usage = "PolynomialRowReduce[m,x] (m \
polynomial matrix in x with full row rank)rreturns a unimodular matrix u s.t. \
r=u.m is a row reduced polynomial matrix.";

Rank::usage = "Rank[a] returns the rank of a matrix a.";

RCoprime::usage = "RCoprime[a,b,x] (a,b polynomial matrices in x \
with the same number of rows) returns True if a and b are right \
coprime.";

RDivision::usage = "RDivision[a,b,x] returns a list {q,r} (quotient \
and remainder after right division of the matrices a (l x m) and b (m x m \
proper) whose elements are polynomials in x, s.t. a=qb+r and \
deg[r]<deg[b].";

RGCD::usage = "RGCD[a,b,x] gives a  right GCD of the matrices a and \
b of polynomials in x. (a and b must have the same number of \
columns).";

Rlcm::usage = "Rlcm[a,b,x] gives a right least common multiple of \
two matrices of polynomials in x. (a and b must have the same number of \
rows).";

    
RQuotient::usage = "RQuotient[a,b,x] gives the right quotient of a \
(l x m) and b (m x m proper) whose entries are polynomials in x. The result q \
satisfies a=q.b+r with deg[r]<deg[b].";

RRemainder::usage = "RRemainder[a,b,x] gives the right remainder \
from division of the matrices a (l x m) and b (m x m proper) of polynomials \
in x. The result r satisfies a=q.b+r and deg[r]<deg[b].";

SmithForm::usage = "SmithForm[a,x] yields the Smith form of the \
matrix a of polynomials in x.";

Error messages

\(\(\(General::mtrx = "The argument of `1` must be a matrix.";\)\(
\)
    \)\)
General::pol = "The elements of the argument(s) must be \
polynomials.";

    
PolyExtendedGCD::egcdz = "PolyExtendedGCD[{0,0}] has'nt a unique \
solution.";

    
General::badarg = "The argument has more than one variable. One \
more specifying an indeterminate is expected.";

    
General::gcd = "The dimensions of the two matrices a, b (`1`, `2` \
respectively), must satisfy the relation `3`.";

    
\(\(\(LDivision::badarg = "The arguments of `1` [a,b,x] must be two matrices a,b with dimensions l x m and `2`, respectivly, and their entries must be polynomials in x. Matrix b must be proper.";\)\(
\)
    \)\)
PolynomialColumnReduce::badarg = "PolynomialColumnReduce requires \
a full column rank matrix.";

    
PolynomialRowReduce::badarg = "PolynomialRowReduce requires a \
full row rank matrix.";

    
DiophantineSolve::badarg = "DiophantineSolve requires 3 \
polynomial matrices with the same number of rows."

Implementation

Begin["`Private`"]

code

The function is used when at has only one variable ore none.It returns the name of the only variable of at if any, or letter "d" if at is numeric.It is used to initialize the variable x  that represents the indeterminate in which the polynomials are defined

Variabile[at_] := 
        Module[{a = at, xt}, 
          xt = Variables[a]; 
If[xt === {}, d, 
            xt[[1]]]];

In all functions if the user does not specify the indeterminate defining the polynomials (xt==={}) and there is more than one variable in the argument,an error message is displayed,otherwise the functions call Variabile to initialize x

ExtPolQ[lt_List, xt_: {}] := 
        Module[{x = xt}, 
If[
            x === {} && Length[Variables[lt]] > 1, 
            Message[General::badarg]; 
Return[$Failed], 
            If[x === {}, 
              x = Variabile[
                  lt]]]; 
  
            Apply[And, PolynomialQ[#, x] & /@ 
                Flatten[lt]] ||   Apply[
              And, PolynomialQ[#, x] & /@ 
                Flatten[
                  Expand[lt] /. {(x^(n : _)) -> (x^-n), 
                      x -> x^-1}]]];

This function gives the extended polynomial GCD of a list of two given polynomials.It is the polynomial analogue of ExtendedGCD

PolyExtendedGCD[ft_, xt_: {}] := 
        Module[{f = ft, x = xt, v, l, n, min, max, dega, degb, gmax, gmin}, 
          If[x === {} && Length[Variables[f]] > 1, 
            Message[General::badarg]; 
Return[$Failed], 
            If[x === {}, x = Variabile[f]]]; 
If[
            f === {0, 0}, Message[PolyExtendedGCD::egcdz]; 
            Return[$Failed]]; 
v = 
            IdentityMatrix[2]; 
dega = 
            Exponent[f[[1]], x]; 
degb = 
            Exponent[f[[2]], x]; 
While[
            f[[1]] =!= 0 && f[[2]] =!= 0, 
            If[dega < degb, min = 1; max = 2; gmin = dega; gmax = degb, 
              min = 2; max = 1; gmin = degb; 
              gmax = dega]; 
l = 
              PolynomialQuotient[f[[max]], f[[min]], 
                x]; 
f = 
              ReplacePart[f, Expand[f[[max]] - f[[min]]*l], 
                max]; 
v = 
              Transpose[ReplacePart[#, 
                      Expand[(#[[max]] - #[[min]]*l)], max] &[
                  Transpose[v]]]; 
dega = 
              Exponent[f[[1]], x]; 
degb = 
              Exponent[f[[2]], x]]; 
Return[
            Expand[{f[[#]], Transpose[v][[#]]}/
                    Coefficient[f[[#]], x, Exponent[f[[#]], x]] &[
                If[f[[1]] === 0, 2, 1]], x]]];

The following function returns True if vector vt has exactly one nonzero component

Onenonzero[vt_, xt_] := 
        Module[{v = vt, x = xt}, 
          Count[v, n_ /; (n =!= 0)] == 1];

\(\(ExtendedSmithForm[at_, xt_: {}] := 
        Module[{a = at, x = xt, p, m, u, v, k, t, i, j, r, q, deg, esp, coef, 
              lpc, upc, g, g1, g2, ch = False}, 
            If[x === {} && Length[Variables[a]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[a]]]; 
If[
              Not[ExtPolQ[a, x]], Message[General::pol]; 
              Return[$Failed]]; 
If[
              Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]], 
              a = a /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = 
                True]; 
 (*Initialize\ the\ variables*) {p, 
                m} = Dimensions[a]; 
u = 
              IdentityMatrix[p]; 
v = 
              IdentityMatrix[m]; 
k = 1; 
            t = False; 
 (*Main\ loop\ on\ \(k : 
                  go\ on\ until\ the\ last\ p - k\ rows\ and\ m - 
                    k\ columns\ become\ zero\)*) 
              While[\(! MatrixQ[
                    a[\([Range[k, p], Range[k, m]]\)], # === 0 &]\) && 
                t == 
                  False,  (*Find\ the\ least\ degree\ element\ in\ the\ remaining\ submatrix*) {i, 
                  j} = \(\(Position[#, 
                        Min[Cases[#, n_ /; n =!= \(-Infinity\), {2}]]] &\)[
                    Exponent[a[\([Range[k, p], Range[k, m]]\)], 
                      x]]\)[\([1]\)]; 
{i, j} = {i, j} + 
                  k - 1; 
 (*bring\ it\ in\ position\ \((k, 
                    k)\), if\ not\ already\ there*) If[i ≠ k, 
                a = a /. {a[\([k]\)] -> a[\([i]\)], 
                      a[\([i]\)] -> a[\([k]\)]}; 
u = 
                  u /. {u[\([k]\)] -> u[\([i]\)], 
                      u[\([i]\)] -> u[\([k]\)]}]; 
If[
                j ≠ k, 
                a = Transpose[\(# /. {#[\([k]\)] -> #[\([j]\)], #[\([j]\)] -> #[\([k]\)]} &\)[Transpose[a]]]; 
v = 
                  Transpose[\(# /. {#[\([k]\)] -> #[\([j]\)], #[\([j]\)] -> #[\([k]\)]} &\)[
                      Transpose[
                        v]]]]; 
 (*The\ current\ element\ is\ a[\([k, k]\)] . Make\ it\ monic*) deg = 
                Exponent[a[\([k, k]\)], x]; 
coef = 
                Coefficient[a[\([k, k]\)], x, deg]; 
If[
                coef =!= 1, 
                a = ReplacePart[a, Expand[a[\([k]\)]/coef], 
                    k]; 
u = 
                  ReplacePart[u, Expand[u[\([k]\)]/coef], 
                    k]]; 
 (*Lower\ the\ degree\ of\ non - 
                  zero\ polynomials\ in\ column\ k\ below\ a[\([k, k]\)]*) If[
                Onenonzero[\(Transpose[a]\)[\([k]\)], x] == False, 
                For[r = k + 1, r ≤ p, \(r++\), 
                  If[a[\([r, k]\)] =!= 0, 
                    lpc = PolynomialQuotient[a[\([r, k]\)], a[\([k, k]\)], 
                        x]; 
a = 
                      ReplacePart[a, Expand[a[\([r]\)] - a[\([k]\)]*lpc], 
                        r]; 
u = 
                      ReplacePart[u, Expand[u[\([r]\)] - u[\([k]\)]*lpc], 
                        r]]],  (*When\ the\ elements\ in\ column\ k\ below\ a[\([k, k]\)]\ are\ all\ zero\ lower\ the\ degree\ of\ non - 
                    zero\ polynomials\ in\ position\ \((k, 
                        j)\), \((j ≥ k)\)*)  (*else*) If[
                  Onenonzero[a[\([k]\)], x] == False, 
                  For[r = k + 1, r ≤ m, \(r++\), 
                    If[a[\([k, r]\)] =!= 0, 
                      lpc = PolynomialQuotient[a[\([k, r]\)], a[\([k, k]\)], 
                          x]; 
a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([r]\)] - #[\([k]\)]*lpc], r] &\)[
                            Transpose[a]]]; 
v = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([r]\)] - #[\([k]\)]*lpc], r] &\)[
                            Transpose[v]]]]],  (*When\ the\ only\ non - 
                      zero\ element\ in\ row\ k\ and\ column\ k\ is\ a[\([k, 
                            k]\)]\ make\ a[\([k - 1, 
                            k - 1]\)]\ a\ divisor\ of\ a[\([k, 
                            k]\)]*)  (*else*) If[k ≠ 1, 
                    If[PolynomialRemainder[a[\([k, k]\)], 
                          a[\([k - 1, k - 1]\)], x] =!= 0, {g, {g1, g2}} = 
                        PolyExtendedGCD[{a[\([k - 1, k - 1]\)], 
                            a[\([k, k]\)]}, x]; 
a = 
                        ReplacePart[a, 
                          Expand[a[\([k]\)] + a[\([k - 1]\)]*g1], 
                          k]; 
u = 
                        ReplacePart[u, 
                          Expand[u[\([k]\)] + u[\([k - 1]\)]*g1], 
                          k]; 
a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k - 1]\)] + #[\([k]\)]*g2], 
                                k - 1] &\)[
                            Transpose[a]]]; 
v = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k - 1]\)] + #[\([k]\)]*g2], 
                                k - 1] &\)[
                            Transpose[v]]]; 
a = 
                        a /. {a[\([k]\)] -> a[\([k - 1]\)], 
                            a[\([k - 1]\)] -> 
                              a[\([k]\)]}; 
u = 
                        u /. {u[\([k]\)] -> u[\([k - 1]\)], 
                            u[\([k - 1]\)] -> 
                              u[\([k]\)]}; 
lpc = 
                        PolynomialQuotient[a[\([k, k - 1]\)], g, 
                          x]; 
a = 
                        ReplacePart[a, 
                          Expand[a[\([k]\)] - a[\([k - 1]\)]*lpc], 
                          k]; 
u = 
                        ReplacePart[u, 
                          Expand[u[\([k]\)] - u[\([k - 1]\)]*lpc], 
                          k]; 
lpc = 
                        PolynomialQuotient[a[\([k - 1, k]\)], g, 
                          x]; 
a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k]\)] - #[\([k - 1]\)]*lpc], 
                                k] &\)[Transpose[a]]]; 
v = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k]\)] - #[\([k - 1]\)]*lpc], 
                                k] &\)[
                            Transpose[v]]]; 
coef = 
                        Coefficient[a[\([k, k]\)], x, 
                          Exponent[a[\([k, k]\)], x]]; 
If[
                        coef =!= 1, 
                        a = ReplacePart[a, Expand[a[\([k]\)]/coef], 
                            k]; 
u = 
                          ReplacePart[u, Expand[u[\([k]\)]/coef], k]]; 
                      k = k - 2]]; 
If[
                    k == p || k == m, 
                    t = True, \(k++\)]]]
 (*endwhile*) ]; 
If[
              ch, {a, u, 
                  v} = {a, u, v} /. {x^\((n : _)\) -> x^\(-n\), 
                    x -> x^\(-1\)}]; 
Return[{a, u, 
                v}]] /; MatrixQ[at] || 
            Message[General::mtrx, "ExtendedSmithForm"];\)\)

Same as above with the only difference that in this case the unimodular matrices u and v are not computed

\(\(SmithForm[at_, xt_: {}] := 
        Module[{a = at, x = xt, p, m, k, t, i, j, r, q, deg, esp, coef, lpc, 
              upc, g, g1, g2, ch = False}, 
            If[x === {} && Length[Variables[a]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[a]]]; 
If[
              Not[ExtPolQ[a, x]], Message[General::pol]; 
              Return[$Failed]]; 
If[
              Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]], 
              a = a /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = True]; 
{p, m} = 
              Dimensions[a]; 
k = 1; 
            t = False; 
While[\(! MatrixQ[
                    a[\([Range[k, p], Range[k, m]]\)], # === 0 &]\) && 
                t == 
                  False,  (*\(\(\(=!=\)\(Table[
                        0, {p - k + 1}, {m - k + 1}]\ For[r = k; 
                        i = \(j = k\); deg = \(+Infinity\), 
                        r ≤ p, \(r++\), 
                        For[q = k, q ≤ m, \(q++\), 
                          esp = Exponent[a[\([r, q]\)], x]; 
                          If[esp ≠ \(-Infinity\) && 
                              esp < deg,  (*then*) deg = esp; i = r; 
                            j = q]]]\)\);\)*) {i, 
                  j} = \(\(Position[#, 
                        Min[Cases[#, n_ /; n =!= \(-Infinity\), {2}]]] &\)[
                    Exponent[a[\([Range[k, p], Range[k, m]]\)], 
                      x]]\)[\([1]\)]; 
{i, j} = {i, j} + 
                  k - 1; 
If[i ≠ k, 
                a = a /. {a[\([k]\)] -> a[\([i]\)], 
                      a[\([i]\)] -> a[\([k]\)]}]; 
If[
                j ≠ k, 
                a = Transpose[\(# /. {#[\([k]\)] -> #[\([j]\)], #[\([j]\)] -> #[\([k]\)]} &\)[Transpose[a]]]]; 
deg = 
                Exponent[a[\([k, k]\)], x]; 
coef = 
                Coefficient[a[\([k, k]\)], x, deg]; 
If[
                coef =!= 1, 
                a = ReplacePart[a, Together[a[\([k]\)]/coef], 
                    k]]; 
If[
                Onenonzero[\(Transpose[a]\)[\([k]\)], x] == False, 
                For[r = k + 1, r ≤ p, \(r++\), 
                  If[a[\([r, k]\)] =!= 0, 
                    lpc = PolynomialQuotient[a[\([r, k]\)], a[\([k, k]\)], 
                        x]; 
a = 
                      ReplacePart[a, Together[a[\([r]\)] - a[\([k]\)]*lpc], 
                        r]]],  (*else*) If[
                  Onenonzero[a[\([k]\)], x] == False, 
                  For[r = k + 1, r ≤ m, \(r++\), 
                    If[a[\([k, r]\)] =!= 0, 
                      lpc = PolynomialQuotient[a[\([k, r]\)], a[\([k, k]\)], 
                          x]; 
a = 
                        Transpose[\(ReplacePart[#, 
                                Together[#[\([r]\)] - #[\([k]\)]*lpc], r] &\)[
                            Transpose[a]]]]],  (*else*) If[k ≠ 1, 
                    If[PolynomialRemainder[a[\([k, k]\)], 
                          a[\([k - 1, k - 1]\)], x] =!= 0, {g, {g1, g2}} = 
                        PolyExtendedGCD[{a[\([k - 1, k - 1]\)], 
                            a[\([k, k]\)]}, x]; 
a = 
                        ReplacePart[a, 
                          Expand[a[\([k]\)] + a[\([k - 1]\)]*g1], 
                          k]; 
a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k - 1]\)] + #[\([k]\)]*g2], 
                                k - 1] &\)[
                            Transpose[a]]]; 
a = 
                        a /. {a[\([k]\)] -> a[\([k - 1]\)], 
                            a[\([k - 1]\)] -> 
                              a[\([k]\)]}; 
lpc = 
                        PolynomialQuotient[a[\([k, k - 1]\)], g, 
                          x]; 
a = 
                        ReplacePart[a, 
                          Expand[a[\([k]\)] - a[\([k - 1]\)]*lpc], 
                          k]; 
lpc = 
                        PolynomialQuotient[a[\([k - 1, k]\)], g, 
                          x]; 
a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k]\)] - #[\([k - 1]\)]*lpc], 
                                k] &\)[
                            Transpose[a]]]; 
coef = 
                        Coefficient[a[\([k, k]\)], x, 
                          Exponent[a[\([k, k]\)], x]]; 
If[
                        coef =!= 1, 
                        a = ReplacePart[a, Expand[a[\([k]\)]/coef], k]]; 
                      k = k - 2]]; 
If[
                    k == p || k == m, 
                    t = True, \(k++\)]]]
 (*endwhile*) ]; 
If[ch, 
              a = a /. {x^\((n : _)\) -> x^\(-n\), 
                    x -> x^\(-1\)}]; 
Return[
              Expand[a]]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "SmithForm"];\)\)

Compute the least degree element in vector at[[j]]/; (j≥k)

FindMin[at_, kt_, xt_, pt_] := 
        Module[{a = at, k = kt, x = xt, p = pt, i, rgmin, degmin}, 
          For[i = k; rgmin = k; degmin = +Infinity, 
            i ≤ p, i++, 
            If[a[[i]] =!= 0 && Exponent[a[[i]], x] < degmin, 
              rgmin = i; 
degmin = 
                Exponent[a[[i]], 
                  x]]]; 
rgmin];

\(\(ExtendedHermiteForm[at_, xt_: {}] := 
        Module[{a = at, p, m, x = xt, u, v, k, kcl, i, j, q, deg, esp, coef, 
              rg, rmin, upc, lpc, ch = False}, 
            If[x === {} && Length[Variables[a]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[a]]]; 
If[
              Not[ExtPolQ[a, x]], Message[General::pol]; 
              Return[$Failed]]; 
If[
              Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]], 
              a = a /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = 
                True]; 
 (*Initialize\ the\ variables*) {p, 
                m} = Dimensions[a]; 
u = 
              IdentityMatrix[
                p]; 
 (*Main\ loop\ on\ k\ \((and\ kcl)\)*) 
              For[k = 1; kcl = 1, 
              k ≤ p && 
                kcl ≤ m, \(k++\),  (*Find\ the\ first\ non - 
                  zero\ element\ in\ row\ k*) While[
                a[\([Range[k, p], kcl]\)] === Table[0, {p - k + 1}] && 
                  kcl < m, \(kcl++\)]; 
 (*With\ this\ loop\ we\ eliminate\ all\ the\ elements\ in\ column\ kcl\ below\ position\ k*) 
                While[If[k == p, False, 
                  a[\([Range[k + 1, p], kcl]\)] =!= 
                    Table[0, {p - 
                          k}]],  (*Put\ the\ least\ degree\ element\ within\ column\ kcl\ in\ position\ \((k, kcl)\)*) rmin = 
                  FindMin[\(Transpose[a]\)[\([kcl]\)], k, x, 
                    p]; 
If[rmin ≠ k, 
                  a = a /. {a[\([rmin]\)] -> a[\([k]\)], 
                        a[\([k]\)] -> 
                          a[\([rmin]\)]}; 
u = 
                    u /. {u[\([rmin]\)] -> u[\([k]\)], 
                        u[\([k]\)] -> 
                          u[\([rmin]\)]}]; 
 (*Lower\ the\ degree\ of\ non - zero\ polynomials\ in\ column\ kcl\ below\ a[\([k, kcl]\)]*) 
                  For[i = k + 1, i ≤ p, \(i++\), 
                  If[a[\([i, kcl]\)] =!= 0, 
                    lpc = PolynomialQuotient[a[\([i, kcl]\)], 
                        a[\([k, kcl]\)], x]; 
a = 
                      ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc], 
                        i]; 
u = 
                      ReplacePart[u, Expand[u[\([i]\)] - u[\([k]\)]*lpc], 
                        
                        i]]]];  (*Endwhile*)  (*In\ column\ kcl\ lower\ the\ degree\ of\ polynomials\ above\ a[\([k, 
                      kcl]\)]\ having\ the\ degree\ higher\ than\ that\ of\ a[\([k, kcl]\)]*) If[a[\([k, kcl]\)] =!= 0, 
                For[i = 1, i < k, \(i++\), 
                  If[a[\([i, kcl]\)] =!= 0, 
                    If[Exponent[a[\([i, kcl]\)], x] ≥ 
                        Exponent[a[\([k, kcl]\)], x], 
                      lpc = PolynomialQuotient[a[\([i, kcl]\)], 
                          a[\([k, kcl]\)], x]; 
a = 
                        ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc], 
                          i]; 
u = 
                        ReplacePart[u, Expand[u[\([i]\)] - u[\([k]\)]*lpc], 
                          i]]]]]; 
 (*Make\ a[\([k, 
                      kcl]\)]\ monic*) If[a[\([k, kcl]\)] =!= 0, 
                lpc = Coefficient[a[\([k, kcl]\)], x, 
                    Exponent[a[\([k, kcl]\)], x]]; 
If[
                  lpc =!= 1, 
                  a = ReplacePart[a, Expand[a[\([k]\)]/lpc], 
                      k]; 
u = 
                    ReplacePart[u, Expand[u[\([k]\)]/lpc], 
                      k]]]; 
\(kcl++\)];  (*Endfor*)  (*Put\ all\ the\ zero\ rows\ in\ the\ last\ positions*) For[k = 1, k < p, \(k++\), 
              If[a[\([k]\)] === Table[0, {m}], 
                i = p; 
While[
                  a[\([i]\)] === Table[0, {m}] && i > k, 
                  i = i - 1]; 
If[i > k, 
                  For[j = k, j ≤ i, \(j++\), 
                    a = ReplacePart[a, a[\([j + 1]\)], 
                        j]; 
u = 
                      ReplacePart[u, u[\([j + 1]\)], 
                        j]]];] (*endif*) ];  (*endfor*) If[
              ch, {a, u} = {a, u} /. {x^\((n : _)\) -> x^\(-n\), 
                    x -> x^\(-1\)}]; 
{a, u}] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "ExtendedHermiteForm"];\)\)

Same as above with the only difference that the unimodular matrix u is not computed

\(\(HermiteForm[at_, xt_: {}] := 
        Module[{a = at, p, m, x = xt, k, kcl, i, j, q, coef, rmin, lpc, 
              ch = False}, 
            If[x === {} && Length[Variables[a]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[a]]]; 
If[
              Not[ExtPolQ[a, x]], Message[General::pol]; 
              Return[$Failed]]; 
If[
              Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]], 
              a = a /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = True]; 
{p, m} = 
              Dimensions[a]; 
For[k = 1; kcl = 1, 
              k ≤ p && kcl ≤ m, \(k++\), 
              While[a[\([Range[k, p], kcl]\)] === Table[0, {p - k + 1}] && 
                  kcl < m, \(kcl++\)]; 
While[
                If[k == p, False, 
                  a[\([Range[k + 1, p], kcl]\)] =!= Table[0, {p - k}]], 
                rmin = FindMin[\(Transpose[a]\)[\([kcl]\)], k, x, 
                    p]; 
If[rmin ≠ k, 
                  a = a /. {a[\([rmin]\)] -> a[\([k]\)], 
                        a[\([k]\)] -> 
                          a[\([rmin]\)]}]; 
For[i = k + 1, 
                  i ≤ p, \(i++\), 
                  If[a[\([i, kcl]\)] =!= 0, 
                    lpc = PolynomialQuotient[a[\([i, kcl]\)], 
                        a[\([k, kcl]\)], x]; 
a = 
                      ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc], 
                        i]]]];  (*Endwhile*) 
If[
                a[\([k, kcl]\)] =!= 0, 
                For[i = 1, i < k, \(i++\), 
                  If[a[\([i, kcl]\)] =!= 0, 
                    While[Exponent[a[\([i, kcl]\)], x] ≥ 
                        Exponent[a[\([k, kcl]\)], x], 
                      lpc = PolynomialQuotient[a[\([i, kcl]\)], 
                          a[\([k, kcl]\)], x]; 
a = 
                        ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc], 
                          i];] (*endwhile*) ]]]; 
If[
                a[\([k, kcl]\)] =!= 0, 
                lpc = Coefficient[a[\([k, kcl]\)], x, 
                    Exponent[a[\([k, kcl]\)], x]]; 
If[
                  lpc =!= 1, 
                  a = ReplacePart[a, Expand[a[\([k]\)]/lpc], 
                      k]]]; 
\(kcl++\)];  (*Endfor*) For[
              k = 1, k < p, \(k++\), 
              If[a[\([k]\)] === Table[0, {m}], 
                i = p; 
While[
                  a[\([i]\)] === Table[0, {m}] && i > k, 
                  i = i - 1]; 
If[i > k, 
                  For[j = k, j ≤ i, \(j++\), 
                    a = ReplacePart[a, a[\([j + 1]\)], 
                        j]]];] (*endif*) ];  (*endfor*) If[ch, 
              a = a /. {x^\((n : _)\) -> x^\(-n\), 
                    x -> x^\(-1\)}]; 
Return[a]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "HermiteForm"];\)\)

ExtendedMcMillanForm[at_, xt_: {}] := 
        Module[{a = at, x = xt, s, u, v, lcm, coef}, 
            If[x === {} && Length[Variables[a]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, 
                x = Variabile[
                    a]]]; 
  lcm = 
              Apply[PolynomialLCM, 
                Flatten[Denominator[Together[a]]]]; 
coef = 
              Coefficient[lcm, x, Exponent[lcm, x]]; 
If[
              coef =!= 1, 
              lcm = lcm/
                  coef]; 
  {s, u, v} = 
              ExtendedSmithForm[Apart[a*lcm], 
                x]; 
  {Together[s/lcm], u, v}] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "ExtendedMcMillanForm"];

Same as above without unimodular matrices

McMillanForm[at_, xt_: {}] := 
        Module[{a = at, x = xt, s, lcm, coef}, 
            If[x === {} && Length[Variables[a]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[a]]]; 
lcm = 
              Apply[PolynomialLCM, 
                Flatten[Denominator[Together[a]]]]; 
coef = 
              Coefficient[lcm, x, Exponent[lcm, x]]; 
If[
              coef =!= 1, lcm = lcm/coef]; 
s = 
              SmithForm[Apart[a*lcm], x]; 
Together[
              s/lcm]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "McMillanForm"];

The following auxiliary function eliminates the dependent rows of the polynomial matrix at

\(\(EliminateDep[at_, xt_: {}] := 
        Module[{a = \(b = at\), x = xt, u, p, m, rmin, k, kcl, ch = False}, 
            If[x === {} && Length[Variables[a]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[a]]]; 
If[
              Not[ExtPolQ[a, x]], Message[General::pol]; 
              Return[$Failed]]; 
If[
              Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]], 
              a = a /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = True]; 
{p, m} = 
              Dimensions[a]; 
u = 
              IdentityMatrix[p]; 
For[k = \(kcl = 1\), 
              k ≤ p && kcl ≤ m, \(k++\), 
              While[a[\([Range[k, p], kcl]\)] === Table[0, {p - k + 1}] && 
                  kcl < m, \(kcl++\)]; 
While[
                If[k == p, False, 
                  a[\([Range[k + 1, p], kcl]\)] =!= Table[0, {p - k}]], 
                rmin = FindMin[\(Transpose[a]\)[\([kcl]\)], k, x, 
                    p]; 
If[rmin ≠ k, 
                  a = a /. {a[\([rmin]\)] -> a[\([k]\)], 
                        a[\([k]\)] -> 
                          a[\([rmin]\)]}; 
u = 
                    u /. {u[\([rmin]\)] -> u[\([k]\)], 
                        u[\([k]\)] -> 
                          u[\([rmin]\)]}]; 
For[i = k + 1, 
                  i ≤ p, \(i++\), 
                  If[a[\([i, kcl]\)] =!= 0, 
                    lpc = PolynomialQuotient[a[\([i, kcl]\)], 
                        a[\([k, kcl]\)], x]; 
a = 
                      ReplacePart[a, Expand[a[\([i]\)] - a[\([k]\)]*lpc, x], 
                        i];]]];  (*Endwhile*) \(kcl++\)];  (*endfor*) k = 
              Flatten[Position[Inverse[u] . a, 
                  n_ /; n === Table[0, {m}]]]; 
If[ch, 
              b = b /. {x^\((n : _)\) -> x^\(-n\), 
                    x -> 
                      x^\(-1\)}]; 
Return[\(IdentityMatrix[
                    p]\)[\([Complement[Range[1, p], k]]\)] . b]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "EliminateDep"];\)\)

The following auxiliary function computes a left GCD and a right lcm of matrices at and bt

AuxLGCDlcm[at_, bt_, xt_] := 
        Module[{a = at, b = bt, x = xt, l, m, q, v, f, i, j, done, krow, 
            lambda, esp, coef, deg, ch = False}, 
          If[Not[ExtPolQ[{a, b}, x]], Message[General::pol]; 
            Return[$Failed]]; 
If[
            Not[Apply[And, PolynomialQ[#, x] & /@ Flatten[{a, b}]]], {a, 
                b} = {a, 
                  b} /. {(x^(n : _)) -> (x^-n)}; 
ch = 
              True]; 
 (*Initialize the variables*) 
{l, m} = 
            Dimensions[
              a]; 
q = Dimensions[
                b][[2]]; 
v = 
            IdentityMatrix[m + q]; 
f = 
            Transpose[
              Flatten[{Transpose[a], Transpose[b]}, 
                1]]; 
 (*Main loop on k, 
            columns of f*) 
For[k = krow = 1, 
            krow ≠ l + 1 && k ≠ m + q, k++, 
            done = False; 
While[
              f[[krow, Range[k, m + q]]] === Table[0, {m + q - k + 1}] && 
                krow < l, krow++]; 
 (*Eliminate all the elements of row krow, starting from column k + 1*) 
              While[f[[krow, Range[k + 1, m + q]]] =!= 
                    Table[0, {m + q - k}] && Not[done], 
                j = FindMin[f[[krow]], k, x, 
                    m + q]; 
If[j ≠ k, 
                  f = Transpose[# /. {#[[k]] -> #[[j]], #[[j]] -> #[[k]]} &[Transpose[f]]]; 
v = 
                    Transpose[# /. {#[[k]] -> #[[j]], #[[j]] -> #[[k]]} &[
                        Transpose[
                          v]]]]; 
 (*If on row krow all the elements following a[[krow, k]] are zero we are done, 
                  otherwise lower the degree of these non - 
                    zero elements*) If[
                  Onenonzero[f[[krow, Range[k, m + q]]], x], done = True, 
                  For[i = k + 1, i ≤ m + q, i++, 
                    If[f[[krow, i]] =!= 0, 
                      lambda = 
                        PolynomialQuotient[f[[krow, i]], 
                          f[[krow, k]], x]; 
f = 
                        Transpose[ReplacePart[#, 
                                Expand[#[[i]] - #[[k]] (lambda), 
                                  x], i] &[
                            Transpose[f]]]; 
v = 
                        Transpose[ReplacePart[#, 
                                Expand[#[[i]] - #[[k]] (lambda), 
                                  x], i] &[
                            Transpose[
                              v]]]]] (*endfor*) ] (*endif*) ] (*endwhile*)  krow++];  (*endfor*)  (*Count the number of last zero columns in f*) j = q + m; 
While[
            j ≥ 1 && Transpose[f][[j]] === 
                Table[0, {l}], 
            j = j - 1]; 
 (*Return the left GCD and if Flatten[lambda] =!= {} return also a right lcm*) If[
            ch, {a, v, 
                f} = {a, v, f} /. {x^(n : _) -> x^-n, 
                  x -> x^-1}]; 
lambda = 
            v[[Range[1, m], Range[j + 1, q + m]]]; 
{If[
              l > (m + q), 
              Transpose[
                Flatten[{Transpose[f[[Range[1, l], Range[1, m + q]]]], 
                    Table[0, {l - m - q}, {l}]}, 1]], 
              f[[Range[1, l], Range[1, l]]]], 
            If[Flatten[lambda] =!= {}, 
              If[l < m + q, 
                Transpose[EliminateDep[Transpose[Expand[a . lambda, x]], x]], 
                Expand[a . lambda, x], {}], {}]}];

LGCD[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt}, 
            If[x === {} && Length[Variables[{a, b}]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, 
                x = Variabile[{a, 
                      b}]]]; 
If[# =!= $Failed, #[[1]], #] &[AuxLGCDlcm[a, b, x]]] /; 
          MatrixQ[at] && 
            MatrixQ[bt] && Dimensions[at][[1]] == Dimensions[
                  bt][[1]];

Rlcm[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt}, 
            If[x === {} && Length[Variables[{a, b}]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, 
                x = Variabile[{a, 
                      b}]]]; 
If[# =!= $Failed, #[[2]], #] &[AuxLGCDlcm[a, b, x]]] /; 
          MatrixQ[at] && 
            MatrixQ[bt] && Dimensions[at][[1]] == Dimensions[
                  bt][[1]];

The following auxiliary function computes a right GCD and a left lcm of matrices at and bt

AuxRGCDlcm[at_, bt_, xt_] := 
        Module[{a = at, b = bt, x = xt, l, m, p, v, f, i, j, done, kcl, 
            lambda, esp, coef, deg, ch = False}, 
          If[Not[ExtPolQ[{a, b}, x]], Message[General::pol]; 
            Return[$Failed]]; 
If[
            Not[Apply[And, PolynomialQ[#, x] & /@ Flatten[{a, b}]]], {a, 
                b} = {a, 
                  b} /. {(x^(n : _)) -> (x^-n)}; 
ch = 
              True]; 
 (*Initialize the variables*) {l, 
              m} = Dimensions[
              a]; 
p = Dimensions[
                b][[1]]; 
v = 
            IdentityMatrix[l + p]; 
f = 
            Flatten[{a, b}, 1]; 
 (*Main loop on k, 
            rows of f*) For[k = kcl = 1, 
            kcl ≠ m + 1 && k ≠ l + p, k++, 
            done = False; 
While[
              f[[Range[k, l + p], kcl]] === Table[0, {l + p - k + 1}] && 
                kcl < m, kcl++]; 
 (*Eliminate all the elements of column kcl, starting from row k + 1*) 
              While[f[[Range[k + 1, l + p], kcl]] =!= 
                    Table[0, {l + p - k}] && Not[done], 
                j = FindMin[Transpose[f][[kcl]], k, x, 
                    l + p]; 
If[j ≠ k, 
                  f = f /. {f[[k]] -> f[[j]], 
                        f[[j]] -> 
                          f[[k]]}; 
v = 
                    v /. {v[[k]] -> v[[j]], 
                        v[[j]] -> 
                          v[[k]]}]; 
 (*If on column kcl all the elements below a[[k, kcl]] are zero we are done, 
                  otherwise lower the degree of these non - 
                    zero elements*) If[
                  Onenonzero[f[[Range[k, l + p], kcl]], x], done = True, 
                  For[i = k + 1, i ≤ l + p, i++, 
                    If[f[[i, kcl]] =!= 0, 
                      lambda = 
                        PolynomialQuotient[f[[i, kcl]], f[[k, kcl]], 
                          x]; 
f = 
                        ReplacePart[f, 
                          Expand[f[[i]] - f[[k]] (lambda), x], 
                          i]; 
v = 
                        ReplacePart[v, 
                          Expand[v[[i]] - v[[k]] (lambda), x], 
                          i]]] (*endfor*) ] (*endif*) ] (*endwhile*)  kcl++];  (*endfor*)  (*Count the number of last zero rows in f*) j = 
            l + p; 
While[
            j ≥ 1 && f[[j]] === Table[0, {m}], 
            j = j - 1]; 
If[
            ch, {a, v, 
                f} = {a, v, f} /. {x^(n : _) -> x^-n, 
                  x -> x^-1}]; 
lambda = 
            v[[Range[j + 1, l + p], 
                Range[1, 
                  l]]]; 
 (*Return the right GCD and if[Flatten[lambda] =!= {} return also a left lcm*) {If[
              m > (p + l), 
              Transpose[
                Flatten[{f[[Range[1, p + l], Range[1, m]]], 
                    Table[0, {m - p - l}, {m}]}, 1]], 
              f[[Range[1, m], Range[1, m]]]], 
            If[Flatten[lambda] =!= {}, 
              If[m < p + l, EliminateDep[Expand[lambda . a, x], x], 
                Expand[lambda . a, x], {}], {}]}];

RGCD[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt}, 
            If[x === {} && Length[Variables[{a, b}]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, 
                x = Variabile[{a, 
                      b}]]]; 
If[# =!= $Failed, #[[1]], #] &[AuxRGCDlcm[a, b, x]]] /; 
          MatrixQ[at] && 
            MatrixQ[bt] && Dimensions[at][[2]] == Dimensions[
                  bt][[2]];

Llcm[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt}, 
            If[x === {} && Length[Variables[{a, b}]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, 
                x = Variabile[{a, 
                      b}]]]; 
If[# =!= $Failed, #[[2]], #] &[AuxRGCDlcm[a, b, x]]] /; 
          MatrixQ[at] && 
            MatrixQ[bt] && Dimensions[at][[2]] == Dimensions[
                  bt][[2]];

RCoprime[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt, d}, 
            If[x === {} && Length[Variables[{a, b}]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[{a, b}]]]; 
If[
              Not[ExtPolQ[{a, b}, x]], Message[General::pol]; 
              Return[$Failed]]; 
If[
              Not[Apply[
                  And, PolynomialQ[#, x] & /@ Flatten[{a, b}]]], {a, 
                    b} = {a, 
                      b} /. {(x^(n : _)) -> (x^-n)};]; 
d = 
              Det[RGCD[a, b]]; 
(FreeQ[d, x] && 
                d ≠ 0)] /; 
          MatrixQ[at] && 
            MatrixQ[bt] && Dimensions[at][[2]] == Dimensions[
                  bt][[2]];

LCoprime[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt, d}, 
            If[x === {} && Length[Variables[{a, b}]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              If[x === {}, x = Variabile[{a, b}]]]; 
If[
              Not[ExtPolQ[{a, b}, x]], Message[General::pol]; 
              Return[$Failed]]; 
If[
              Not[Apply[
                  And, PolynomialQ[#, x] & /@ Flatten[{a, b}]]], {a, 
                    b} = {a, 
                      b} /. {(x^(n : _)) -> (x^-n)};]; 
d = 
              Det[LGCD[a, b]]; 
(FreeQ[d, x] && 
                d ≠ 0)] /; 
          MatrixQ[at] && 
            MatrixQ[bt] && Dimensions[at][[1]] == Dimensions[
                  bt][[1]];

\(\(LDivision[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt, l, m, u, esp, lambda, delta, 
            ch = False}, 
          If[x === {} && Length[Variables[{a, b}]] > 1, 
            Message[General::badarg]; 
Return[$Failed], 
            If[x === {}, x = Variabile[{a, b}]]]; 
If[
            Not[MatrixQ[at] && 
                MatrixQ[
                  bt] && \(Dimensions[at]\)[\([1]\)] === \(Dimensions[
                      bt]\)[\([1]\)] === \(Dimensions[bt]\)[\([2]\)] && 
                Det[Coefficient[b, x, Max[Exponent[b, x]]]] =!= 0], 
            Message[LDivision::badarg, "LDivision", "l x l"]; 
            Return[$Failed]]; 
If[Not[ExtPolQ[{a, b}, x]], 
            Message[General::pol]; Return[$Failed]]; 
If[
            Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a, 
                b} = {a, 
                  b} /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = True]; 
 (*Inizialization : 
              delta\ is\ the\ inverse\ of\ the\ leading\ coefficient\ matrix\ of\ matrix\ b*) 
{l, m} = 
            Dimensions[a]; 
u = 
            Table[0, {l}, {m}]; 
delta = 
            Inverse[Coefficient[b, x, 
                Max[Exponent[b, x]]]]; 
 (*Main\ \(loop : 
                go\ on\ until\ deg[a] < deg[b]\)*) While[
            Max[Exponent[a, x]] ≥ Max[Exponent[b, x]], 
            lambda = 
              delta . Coefficient[a, x, 
                  Max[Exponent[a, x]]]; 
esp = 
              Max[Exponent[a, x]] - 
                Max[Exponent[b, 
                    x]]; 
 (*Decrease\ the\ degree\ of\ matrix\ a*) a = Expand[a - b . \((lambda*x^esp)\), x]; 
u = 
              Expand[u + lambda*x^esp, x]]; 
If[
            ch, {u, a} = {u, a} /. {x^\((n : _)\) -> x^\(-n\), 
                  x -> x^\(-1\)}]; 
{u, 
            a}];\)\)

LQuotient[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, 
            x = xt}, If[# =!= $Failed, #[[1]], #] &[
            LDivision[a, b, x]]];

LRemainder[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, 
            x = xt}, If[# =!= $Failed, #[[2]], #] &[
            LDivision[a, b, x]]];

\(\(RDivision[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, x = xt, l, m, u, esp, lambda, delta, 
            ch = False}, 
          If[x === {} && Length[Variables[{a, b}]] > 1, 
            Message[General::badarg]; 
Return[$Failed], 
            If[x === {}, x = Variabile[{a, b}]]]; 
If[
            Not[MatrixQ[at] && 
                MatrixQ[
                  bt] && \(Dimensions[at]\)[\([2]\)] === \(Dimensions[
                      bt]\)[\([1]\)] === \(Dimensions[bt]\)[\([2]\)] && 
                Det[Coefficient[bt, x, Max[Exponent[bt, x]]]] =!= 0], 
            Message[LDivision::badarg, "RDivision", "m x m"]; 
            Return[$Failed]]; 
If[Not[ExtPolQ[{a, b}, x]], 
            Message[General::pol]; Return[$Failed]]; 
If[
            Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a, 
                b} = {a, 
                  b} /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = True]; 
 (*Some\ inizialization, 
            delta\ is\ the\ inverse\ of\ the\ leading\ coefficient\ matrix\ of\ matrix\ b*) {l, m} = Dimensions[a]; 
u = 
            Table[0, {l}, {m}]; 
delta = 
            Inverse[Coefficient[b, x, 
                Max[Exponent[b, x]]]]; 
 (*Main\ \(loop : 
                go\ on\ until\ deg[a] < deg[b]\)*) While[
            Max[Exponent[a, x]] ≥ Max[Exponent[b, x]], 
            lambda = 
              Coefficient[a, x, Max[Exponent[a, x]]] . 
                delta; 
esp = 
              Max[Exponent[a, x]] - 
                Max[Exponent[b, 
                    x]]; 
 (*Decrease\ the\ degree\ of\ matrix\ a*) a = Expand[a - \((lambda*x^esp)\) . b, x]; 
u = 
              Expand[u + lambda*x^esp, x]]; 
If[
            ch, {u, a} = {u, a} /. {x^\((n : _)\) -> x^\(-n\), 
                  x -> x^\(-1\)}]; 
{u, 
            a}];\)\)

RQuotient[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, 
            x = xt}, If[# =!= $Failed, #[[1]], #] &[
            RDivision[a, b, x]]];

RRemainder[at_, bt_, xt_: {}] := 
        Module[{a = at, b = bt, 
            x = xt}, If[# =!= $Failed, #[[2]], #] &[
            RDivision[a, b, x]]];

Rank[at_] := 
        Module[{a = at, x = xt, p, m}, {p, m} = 
              Dimensions[a]; 
a = 
              RowReduce[a]; 
Return[
              p - Count[a, n_ /; n === Table[0, {m}]]]] /; 
          MatrixQ[at] || Message[General::mtrx, "Rank"];

FullRowRankQ[at_] := 
        Module[{p = Dimensions[at][[1]]}, Return[Rank[at] === p]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "FullRowRankQ"];

\(\(PolynomialColumnReduce[at_, xt_: {}] := 
        Module[{a = at, x = xt, u, p, m, Ahc, v, deg, degmax, k, ch = False}, 
            If[x === {}, 
              If[Length[Variables[a]] > 1, 
                Message[General::badarg]; 
Return[$Failed], 
                x = Variabile[a]]]; 
If[Not[ExtPolQ[a, x]], 
              Message[General::pol]; Return[$Failed]]; 
If[
              Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]], 
              a = a /. {\((x^\((n : _)\))\) -> \((x^\(-n\))\)}; 
ch = 
                True]; 
If[\(! FullRowRankQ[
                  Transpose[a]]\), Message[PolynomialColumnReduce::badarg]; 
              Return[$Failed]]; 
{p, m} = 
              Dimensions[a]; 
u = 
              IdentityMatrix[
                m]; 
 (*Ahc\ is\ the\ matrix\ composed\ by\ the\ coefficients\ of\ the\ highest\ degree\ terms\ within\ each\ column*) 
              Ahc = Transpose[\(Table[
                      Coefficient[#[\([i]\)], x, 
                        Max[Exponent[#[\([i]\)], x]]], {i, 1, m}] &\)[
                  Transpose[
                    a]]]; 
 (*go\ on\ until\ Ahc\ has\ full\ column\ rank\ \((i . e . a\ is\ column\ reduced)\)*) 
              While[\(! FullRowRankQ[Transpose[Ahc]]\), 
              v = \(NullSpace[Ahc]\)[\([1]\)]; 
deg = 
                Table[Max[Exponent[\(Transpose[a]\)[\([i]\)], x]], {i, 1, 
                    m}]; 
k = \(Flatten[
                    Position[deg, 
                      degmax = 
                        Max[Table[
                            Exponent[\(Transpose[a]\)[\([i]\)]*v[\([i]\)], 
                              x], {i, 1, 
                              m}]]]]\)[\([1]\)]; 
 (*lower\ the\ degree\ of\ column\ k*) {u, a} = \({u . #, a . #} &\)[
                  Transpose[
                    ReplacePart[IdentityMatrix[m], 
                      Table[v[\([i]\)]*x^\((degmax - deg[\([i]\)])\), {i, 1, 
                          m}], k]]]; 
Ahc = 
                Transpose[\(Table[
                        Coefficient[#[\([i]\)], x, 
                          Max[Exponent[#[\([i]\)], x]]], {i, 1, m}] &\)[
                    Transpose[a]]]]; 
If[ch, 
              u = u /. {x^\((n : _)\) -> x^\(-n\), 
                    x -> x^\(-1\)}]; 
Return[u]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "PolynomialColumnReduce"];\)\)

\(\(PolynomialRowReduce[at_, xt_: {}] := 
        Module[{a = at, x = xt}, 
            If[x === {}, 
              If[Length[Variables[a]] > 1, 
                Message[General::badarg]; 
Return[$Failed], 
                x = Variabile[a]]]; 
If[Not[ExtPolQ[a, x]], 
              Message[General::pol]; 
              Return[$Failed]]; 
If[\(! FullRowRankQ[a]\), 
              Message[PolynomialRowReduce::badarg]; 
Return[$Failed]]; 
Return[
              Transpose[PolynomialColumnReduce[Transpose[a], x]]]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "PolynomialRowReduce"];\)\)

Test whether a diophantine equation has a solution

EquivQ[R1t_, R2t_, xt_: {}] := 
        Module[{R1 = R1t, R2 = R2t, x = xt, h1, h2}, 
          h1 = HermiteForm[R1, x]; 
h2 = 
            HermiteForm[R2, x]; 
If[h1 === h2, 
            Return[True], Return[False]]];

DiophantineSolve[at_, bt_, ct_, xt_: {}, dt_: {}] := 
        Module[{a = at, b = bt, c = ct, x = xt, Zsol, Ysol, t, i, k, sol, 
            deg, z, y, Z, Y, dummy = dt, ch = False}, 
          If[! MatrixQ[a] || ! MatrixQ[b] || ! MatrixQ[c] || 
              Not[Dimensions[a][[1]] == Dimensions[
                      b][[1]] == Dimensions[c][[1]]], 
            Message[DiophantineSolve::badarg]; 
            Return[$Failed]]; 
If[x === {}, 
            If[Length[Variables[{a, b, c}]] > 1, 
              Message[General::badarg]; 
Return[$Failed], 
              x = Variabile[{a, b, c}]]]; 
If[
            Not[ExtPolQ[{a, b, c}, x]], Message[General::pol]; 
            Return[$Failed]]; 
If[
            Not[Apply[
                And, PolynomialQ[#, x] & /@ Flatten[{a, b, c}]]], {a, b, 
                c} = {a, b, 
                  c} /. {(x^(n : _)) -> (x^-n)}; 
ch = True]; 
 (*Ysol, 
            Zsol are the polynomial matrix solutions of the equation a . Y + b . Z == c*) 
            Ysol = {{}}; 
Zsol = {{}}; 
 (*Test wether the equation has solutions or not*) If[
            Rank[Transpose[
                    Flatten[{Transpose[a], Transpose[b]}, 1]]] == 
                Rank[Transpose[
                    Flatten[{Transpose[a], Transpose[b], Transpose[c]}, 1]]] && 
              EquivQ[Flatten[{Transpose[a], Transpose[b], 
                    Table[0, {Dimensions[c][[2]]}, {Dimensions[
                            a][[1]]}]}, 1], 
                Flatten[{Transpose[a], Transpose[b], Transpose[c]}, 1], 
                x],  (*Main loop on t*) For[t = 1, 
              t ≤ Dimensions[c][[2]], t++, 
              sol = {}; 
deg = 
                0; 
 (*Find column t of the minimal column degree solutions Ysol and Zsol*) While[
                sol === {},  (*define column t of Ysol and Zsol, 
                  composed by polynomials with degree deg*) Y = 
                  Table[Sum[
                      y[i, k] x^k, {k, 0, 
                        deg}], {i, Dimensions[
                          a][[2]]}, {1}]; 
Z = 
                  Table[Sum[
                      z[i, k] x^k, {k, 0, 
                        deg}], {i, Dimensions[
                          
                          b][[2]]}, {1}]; 
 (*try to calculate column t of the solutions using polynomials with degree deg . If no solution exists (sol === {}) repeat the instructions with deg increased by 1*) sol = 
                  SolveAlways[
                    a . Y + b . Z == Transpose[c][[t]], 
                    Variables[{a, b, 
                        c}]]; 
deg++];  (*endwhile*)  (*Update the solutions with the new minimal degree columns just found*) Ysol = Transpose[
                  Flatten[{Transpose[Ysol], 
                      Transpose[(Y /. 
                              sol)[[1]]
 (*/.Map[# -> 0 &, Complement[Variables[{Z, Y}], {x}]]*) ]}, 
                    1]]; 
Zsol = 
                Transpose[
                  Flatten[{Transpose[Zsol], 
                      Transpose[(Z /. 
                              sol)[[1]]
 (*/.Map[# -> 0 &, Complement[Variables[{Z, Y}], {x}]]*) ]}, 
                    1]];];  (*endfor*) ];  (*endif*)  (*if dummy =!= {} rename the free variables of the result using the symbol specified in dummy . 
                  In this case you calculate all the minimal column degree solutions of the equation*) If[
            dummy =!= {}, For[i = 1, i ≤ Length[#], i++, 
                  Module[{t = dummy}, {Ysol, 
                        Zsol} = {Ysol, Zsol} /. #[[i]] -> t[i]]] &[
              Complement[Variables[{Ysol, Zsol}], 
                Variables[{a, b, c}]]], {Ysol, Zsol} = {Ysol, Zsol} /. 
                  Map[# -> 0 &, 
                    Complement[
                      Variables[{Ysol, 
                          Zsol}], {x}]];]; 
If[
            ch, {Ysol, 
                Zsol} = {Ysol, Zsol} /. {x^(n : _) -> x^-n, 
                  x -> x^-1}]; 
Return[{Ysol, 
              Zsol}]];

PolynomialDegree[at_List, xt_: {}] := 
        Module[{x = xt}, 
          If[x === {} && Length[Variables[at]] > 1, 
            Message[General::badarg]; 
Return[$Failed], 
            If[x === {}, x = Variabile[at]]]; 
Return[
            Max[Exponent[at, x]]]];

Epilog

End[]

EndPackage[]