polynomialmatrices.nb

单模多项式矩阵的分解算法

"Input",
  InitializationCell->True],

Cell[BoxData[
    \(\(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.\>";\)\)], "Input",
  InitializationCell->True],

Cell[BoxData[
    \(\(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].\>";\)\)], "Input",
  InitializationCell->True],

Cell[BoxData[
    \(\(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).\>";\)\)], "Input",
  InitializationCell->True],

Cell[BoxData[{
    \(\(\(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).\>";\)\(\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(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].\>";\)\)}], "Input",
  InitializationCell->True],

Cell[BoxData[
    \(\(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].\>";\)\)], "Input",
  InitializationCell->True],

Cell[BoxData[
    \(\(SmithForm::usage = "\<SmithForm[a,x] yields the Smith form of the \
matrix a of polynomials in x.\>";\)\)], "Input",
  InitializationCell->True]
}, Open  ]],

Cell[CellGroupData[{

Cell["Error messages", "Subsubsection"],

Cell[BoxData[{
    \(\(\(General::mtrx = "\<The argument of `1` must be a matrix.\>";\)\(\
\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(General::pol = "\<The elements of the argument(s) must be \
polynomials.\>";\)\(\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(PolyExtendedGCD::egcdz = "\<PolyExtendedGCD[{0,0}] has'nt a unique \
solution.\>";\)\(\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(General::badarg = "\<The argument has more than one variable. One \
more specifying an indeterminate is expected.\>";\)\(\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(General::gcd = "\<The dimensions of the two matrices a, b (`1`, `2` \
respectively), must satisfy the relation `3`.\>";\)\(\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(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.\>";\)\(\[IndentingNewLine]\
\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(PolynomialColumnReduce::badarg = "\<PolynomialColumnReduce requires \
a full column rank matrix.\>";\)\(\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(PolynomialRowReduce::badarg = "\<PolynomialRowReduce requires a \
full row rank matrix.\>";\)\(\[IndentingNewLine]\)
    \)\), "\[IndentingNewLine]", 
    \(\(\(DiophantineSolve::badarg = "\<DiophantineSolve requires 3 \
polynomial matrices with the same number of rows.\>"\)\(\[IndentingNewLine]\)
    \)\)}], "Input",
  InitializationCell->True]
}, Open  ]]
}, Open  ]],

Cell[CellGroupData[{

Cell["Implementation", "Section"],

Cell[BoxData[
    \(Begin["\<`Private`\>"]\)], "Input",
  InitializationCell->True]
}, Open  ]],

Cell[CellGroupData[{

Cell["code", "Subtitle"],

Cell["\<\
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\
\>", "SmallText"],

Cell[BoxData[
    \(\(Variabile[at_] := 
        Module[{a = at, xt}, 
          xt = Variables[a]; \[IndentingNewLine]If[xt === {}, d, 
            xt[\([1]\)]]];\)\)], "Input",
  InitializationCell->True],

Cell["\<\
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\
\>", "SmallText"],

Cell[BoxData[
    \(\(ExtPolQ[lt_List, xt_: {}] := 
        Module[{x = xt}, \[IndentingNewLine]If[
            x === {} && Length[Variables[lt]] > 1, 
            Message[General::badarg]; \[IndentingNewLine]Return[$Failed], 
            If[x === {}, 
              x = Variabile[
                  lt]]]; \[IndentingNewLine] (*If\ all\ the\ components\ of\ \
lt\ are\ polynomials\ in\ x*) 
            Apply[And, \(PolynomialQ[#, x] &\) /@ 
                Flatten[lt]] ||  (*or\ in\ x^\(-1\), 
              the\ function\ returns\ True, 
              False\ otherwise . \((Substituting\ all\ the\ symbols\ x^\(-n\) \
\[Rule] x^n\ we\ can\ use\ PolynomialQ\ also\ in\ this\ case)\)*) Apply[
              And, \(PolynomialQ[#, x] &\) /@ 
                Flatten[
                  Expand[lt] /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\), 
                      x \[Rule] x^\(-1\)}]]];\)\)], "Input",
  InitializationCell->True],

Cell["\<\
This function gives the extended polynomial GCD of a list of two \
given polynomials.It is the polynomial analogue of ExtendedGCD\
\>", \
"SmallText"],

Cell[BoxData[
    \(\(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]; \[IndentingNewLine]Return[$Failed], 
            If[x === {}, x = Variabile[f]]]; \[IndentingNewLine]If[
            f === {0, 0}, Message[PolyExtendedGCD::egcdz]; 
            Return[$Failed]]; \[IndentingNewLine]v = 
            IdentityMatrix[2]; \[IndentingNewLine]dega = 
            Exponent[f[\([1]\)], x]; \[IndentingNewLine]degb = 
            Exponent[f[\([2]\)], x]; \[IndentingNewLine]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]; \[IndentingNewLine]l = 
              PolynomialQuotient[f[\([max]\)], f[\([min]\)], 
                x]; \[IndentingNewLine]f = 
              ReplacePart[f, Expand[f[\([max]\)] - f[\([min]\)]*l], 
                max]; \[IndentingNewLine]v = 
              Transpose[\(ReplacePart[#, 
                      Expand[\((#[\([max]\)] - #[\([min]\)]*l)\)], max] &\)[
                  Transpose[v]]]; \[IndentingNewLine]dega = 
              Exponent[f[\([1]\)], x]; \[IndentingNewLine]degb = 
              Exponent[f[\([2]\)], x]]; \[IndentingNewLine]Return[
            Expand[\({f[\([#]\)], \(Transpose[v]\)[\([#]\)]}/
                    Coefficient[f[\([#]\)], x, Exponent[f[\([#]\)], x]] &\)[
                If[f[\([1]\)] === 0, 2, 1]], x]]];\)\)], "Input",
  InitializationCell->True],

Cell["\<\
The following function returns True if vector vt has exactly one \
nonzero component\
\>", "SmallText"],

Cell[BoxData[
    \(\(Onenonzero[vt_, xt_] := 
        Module[{v = vt, x = xt}, 
          Count[v, n_ /; \((n =!= 0)\)] \[Equal] 1];\)\)], "Input",
  InitializationCell->True],

Cell[BoxData[
    \(\(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]; \[IndentingNewLine]Return[$Failed], 
              If[x === {}, x = Variabile[a]]]; \[IndentingNewLine]If[
              Not[ExtPolQ[a, x]], Message[General::pol]; 
              Return[$Failed]]; \[IndentingNewLine]If[
              Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]], 
              a = a /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \
\[IndentingNewLine]ch = 
                True]; \[IndentingNewLine] (*Initialize\ the\ variables*) {p, 
                m} = Dimensions[a]; \[IndentingNewLine]u = 
              IdentityMatrix[p]; \[IndentingNewLine]v = 
              IdentityMatrix[m]; \[IndentingNewLine]k = 1; 
            t = False; \[IndentingNewLine] (*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 \[Equal] 
                  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]\)]; \[IndentingNewLine]{i, j} = {i, j} + 
                  k - 1; \[IndentingNewLine] (*bring\ it\ in\ position\ \((k, 
                    k)\), if\ not\ already\ there*) If[i \[NotEqual] k, 
                a = a /. {a[\([k]\)] \[Rule] a[\([i]\)], 
                      a[\([i]\)] \[Rule] a[\([k]\)]}; \[IndentingNewLine]u = 
                  u /. {u[\([k]\)] \[Rule] u[\([i]\)], 
                      u[\([i]\)] \[Rule] u[\([k]\)]}]; \[IndentingNewLine]If[
                j \[NotEqual] k, 
                a = Transpose[\(# /. {#[\([k]\)] \[Rule] #[\([j]\)], \
#[\([j]\)] \[Rule] #[\([k]\)]} &\)[Transpose[a]]]; \[IndentingNewLine]v = 
                  Transpose[\(# /. {#[\([k]\)] \[Rule] #[\([j]\)], #[\([j]\)] \
\[Rule] #[\([k]\)]} &\)[
                      Transpose[
                        v]]]]; \[IndentingNewLine] (*The\ current\ element\ \
is\ a[\([k, k]\)] . Make\ it\ monic*) deg = 
                Exponent[a[\([k, k]\)], x]; \[IndentingNewLine]coef = 
                Coefficient[a[\([k, k]\)], x, deg]; \[IndentingNewLine]If[
                coef =!= 1, 
                a = ReplacePart[a, Expand[a[\([k]\)]/coef], 
                    k]; \[IndentingNewLine]u = 
                  ReplacePart[u, Expand[u[\([k]\)]/coef], 
                    k]]; \[IndentingNewLine] (*Lower\ the\ degree\ of\ non - 
                  zero\ polynomials\ in\ column\ k\ below\ a[\([k, k]\)]*) If[
                Onenonzero[\(Transpose[a]\)[\([k]\)], x] \[Equal] False, 
                For[r = k + 1, r \[LessEqual] p, \(r++\), 
                  If[a[\([r, k]\)] =!= 0, 
                    lpc = PolynomialQuotient[a[\([r, k]\)], a[\([k, k]\)], 
                        x]; \[IndentingNewLine]a = 
                      ReplacePart[a, Expand[a[\([r]\)] - a[\([k]\)]*lpc], 
                        r]; \[IndentingNewLine]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 \[GreaterEqual] k)\)*)  (*else*) If[
                  Onenonzero[a[\([k]\)], x] \[Equal] False, 
                  For[r = k + 1, r \[LessEqual] m, \(r++\), 
                    If[a[\([k, r]\)] =!= 0, 
                      lpc = PolynomialQuotient[a[\([k, r]\)], a[\([k, k]\)], 
                          x]; \[IndentingNewLine]a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([r]\)] - #[\([k]\)]*lpc], r] &\)[
                            Transpose[a]]]; \[IndentingNewLine]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 \[NotEqual] 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]; \[IndentingNewLine]a = 
                        ReplacePart[a, 
                          Expand[a[\([k]\)] + a[\([k - 1]\)]*g1], 
                          k]; \[IndentingNewLine]u = 
                        ReplacePart[u, 
                          Expand[u[\([k]\)] + u[\([k - 1]\)]*g1], 
                          k]; \[IndentingNewLine]a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k - 1]\)] + #[\([k]\)]*g2],