polynomialmatrices.nb

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

                                k - 1] &\)[
                            Transpose[a]]]; \[IndentingNewLine]v = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k - 1]\)] + #[\([k]\)]*g2], 
                                k - 1] &\)[
                            Transpose[v]]]; \[IndentingNewLine]a = 
                        a /. {a[\([k]\)] \[Rule] a[\([k - 1]\)], 
                            a[\([k - 1]\)] \[Rule] 
                              a[\([k]\)]}; \[IndentingNewLine]u = 
                        u /. {u[\([k]\)] \[Rule] u[\([k - 1]\)], 
                            u[\([k - 1]\)] \[Rule] 
                              u[\([k]\)]}; \[IndentingNewLine]lpc = 
                        PolynomialQuotient[a[\([k, k - 1]\)], g, 
                          x]; \[IndentingNewLine]a = 
                        ReplacePart[a, 
                          Expand[a[\([k]\)] - a[\([k - 1]\)]*lpc], 
                          k]; \[IndentingNewLine]u = 
                        ReplacePart[u, 
                          Expand[u[\([k]\)] - u[\([k - 1]\)]*lpc], 
                          k]; \[IndentingNewLine]lpc = 
                        PolynomialQuotient[a[\([k - 1, k]\)], g, 
                          x]; \[IndentingNewLine]a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k]\)] - #[\([k - 1]\)]*lpc], 
                                k] &\)[Transpose[a]]]; \[IndentingNewLine]v = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k]\)] - #[\([k - 1]\)]*lpc], 
                                k] &\)[
                            Transpose[v]]]; \[IndentingNewLine]coef = 
                        Coefficient[a[\([k, k]\)], x, 
                          Exponent[a[\([k, k]\)], x]]; \[IndentingNewLine]If[
                        coef =!= 1, 
                        a = ReplacePart[a, Expand[a[\([k]\)]/coef], 
                            k]; \[IndentingNewLine]u = 
                          ReplacePart[u, Expand[u[\([k]\)]/coef], k]]; 
                      k = k - 2]]; \[IndentingNewLine]If[
                    k \[Equal] p || k \[Equal] m, 
                    t = True, \(k++\)]]]\[IndentingNewLine] (*endwhile*) ]; \
\[IndentingNewLine]If[
              ch, {a, u, 
                  v} = {a, u, v} /. {x^\((n : _)\) \[Rule] x^\(-n\), 
                    x \[Rule] x^\(-1\)}]; \[IndentingNewLine]Return[{a, u, 
                v}]] /; MatrixQ[at] || 
            Message[General::mtrx, "\<ExtendedSmithForm\>"];\)\)], "Input",
  InitializationCell->True],

Cell["\<\
Same as above with the only difference that in this case the \
unimodular matrices u and v are not computed\
\>", "SmallText"],

Cell[BoxData[
    \(\(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]; \[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]{p, m} = 
              Dimensions[a]; \[IndentingNewLine]k = 1; 
            t = False; \[IndentingNewLine]While[\(! MatrixQ[
                    a[\([Range[k, p], Range[k, m]]\)], # === 0 &]\) && 
                t \[Equal] 
                  False,  (*\(\(\(=!=\)\(Table[
                        0, {p - k + 1}, {m - k + 1}]\ For[r = k; 
                        i = \(j = k\); deg = \(+Infinity\), 
                        r \[LessEqual] p, \(r++\), 
                        For[q = k, q \[LessEqual] m, \(q++\), 
                          esp = Exponent[a[\([r, q]\)], x]; 
                          If[esp \[NotEqual] \(-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]\)]; \[IndentingNewLine]{i, j} = {i, j} + 
                  k - 1; \[IndentingNewLine]If[i \[NotEqual] k, 
                a = a /. {a[\([k]\)] \[Rule] a[\([i]\)], 
                      a[\([i]\)] \[Rule] a[\([k]\)]}]; \[IndentingNewLine]If[
                j \[NotEqual] k, 
                a = Transpose[\(# /. {#[\([k]\)] \[Rule] #[\([j]\)], \
#[\([j]\)] \[Rule] #[\([k]\)]} &\)[Transpose[a]]]]; \[IndentingNewLine]deg = 
                Exponent[a[\([k, k]\)], x]; \[IndentingNewLine]coef = 
                Coefficient[a[\([k, k]\)], x, deg]; \[IndentingNewLine]If[
                coef =!= 1, 
                a = ReplacePart[a, Together[a[\([k]\)]/coef], 
                    k]]; \[IndentingNewLine]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, Together[a[\([r]\)] - a[\([k]\)]*lpc], 
                        r]]],  (*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[#, 
                                Together[#[\([r]\)] - #[\([k]\)]*lpc], r] &\)[
                            Transpose[a]]]]],  (*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]a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k - 1]\)] + #[\([k]\)]*g2], 
                                k - 1] &\)[
                            Transpose[a]]]; \[IndentingNewLine]a = 
                        a /. {a[\([k]\)] \[Rule] a[\([k - 1]\)], 
                            a[\([k - 1]\)] \[Rule] 
                              a[\([k]\)]}; \[IndentingNewLine]lpc = 
                        PolynomialQuotient[a[\([k, k - 1]\)], g, 
                          x]; \[IndentingNewLine]a = 
                        ReplacePart[a, 
                          Expand[a[\([k]\)] - a[\([k - 1]\)]*lpc], 
                          k]; \[IndentingNewLine]lpc = 
                        PolynomialQuotient[a[\([k - 1, k]\)], g, 
                          x]; \[IndentingNewLine]a = 
                        Transpose[\(ReplacePart[#, 
                                Expand[#[\([k]\)] - #[\([k - 1]\)]*lpc], 
                                k] &\)[
                            Transpose[a]]]; \[IndentingNewLine]coef = 
                        Coefficient[a[\([k, k]\)], x, 
                          Exponent[a[\([k, k]\)], x]]; \[IndentingNewLine]If[
                        coef =!= 1, 
                        a = ReplacePart[a, Expand[a[\([k]\)]/coef], k]]; 
                      k = k - 2]]; \[IndentingNewLine]If[
                    k \[Equal] p || k \[Equal] m, 
                    t = True, \(k++\)]]]\[IndentingNewLine] (*endwhile*) ]; \
\[IndentingNewLine]If[ch, 
              a = a /. {x^\((n : _)\) \[Rule] x^\(-n\), 
                    x \[Rule] x^\(-1\)}]; \[IndentingNewLine]Return[
              Expand[a]]] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "\<SmithForm\>"];\)\)], "Input",
  InitializationCell->True],

Cell["\<\
Compute the least degree element in vector at[[j]]/; (j\
\[GreaterEqual]k)\
\>", "SmallText"],

Cell[BoxData[
    \(\(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 \[LessEqual] p, \(i++\), 
            If[a[\([i]\)] =!= 0 && Exponent[a[\([i]\)], x] < degmin, 
              rgmin = i; \[IndentingNewLine]degmin = 
                Exponent[a[\([i]\)], 
                  x]]]; \[IndentingNewLine]rgmin];\)\)], "Input",
  InitializationCell->True],

Cell[BoxData[
    \(\(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]; \[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] (*Main\ loop\ on\ k\ \((and\ kcl)\)*) 
              For[k = 1; kcl = 1, 
              k \[LessEqual] p && 
                kcl \[LessEqual] 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++\)]; \[IndentingNewLine] (*With\ this\ loop\
\ we\ eliminate\ all\ the\ elements\ in\ column\ kcl\ below\ position\ k*) 
                While[If[k \[Equal] 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]; \[IndentingNewLine]If[rmin \[NotEqual] k, 
                  a = a /. {a[\([rmin]\)] \[Rule] a[\([k]\)], 
                        a[\([k]\)] \[Rule] 
                          a[\([rmin]\)]}; \[IndentingNewLine]u = 
                    u /. {u[\([rmin]\)] \[Rule] u[\([k]\)], 
                        u[\([k]\)] \[Rule] 
                          u[\([rmin]\)]}]; \[IndentingNewLine] (*Lower\ the\ \
degree\ of\ non - zero\ polynomials\ in\ column\ kcl\ below\ a[\([k, kcl]\)]*) 
                  For[i = k + 1, i \[LessEqual] p, \(i++\), 
                  If[a[\([i, kcl]\)] =!= 0, 
                    lpc = PolynomialQuotient[a[\([i, kcl]\)], 
                        a[\([k, kcl]\)], x]; \[IndentingNewLine]a = 
                      ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc], 
                        i]; \[IndentingNewLine]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] \[GreaterEqual] 
                        Exponent[a[\([k, kcl]\)], x], 
                      lpc = PolynomialQuotient[a[\([i, kcl]\)], 
                          a[\([k, kcl]\)], x]; \[IndentingNewLine]a = 
                        ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc], 
                          i]; \[IndentingNewLine]u = 
                        ReplacePart[u, Expand[u[\([i]\)] - u[\([k]\)]*lpc], 
                          i]]]]]; \[IndentingNewLine] (*Make\ a[\([k, 
                      kcl]\)]\ monic*) If[a[\([k, kcl]\)] =!= 0, 
                lpc = Coefficient[a[\([k, kcl]\)], x, 
                    Exponent[a[\([k, kcl]\)], x]]; \[IndentingNewLine]If[
                  lpc =!= 1, 
                  a = ReplacePart[a, Expand[a[\([k]\)]/lpc], 
                      k]; \[IndentingNewLine]u = 
                    ReplacePart[u, Expand[u[\([k]\)]/lpc], 
                      k]]]; \[IndentingNewLine]\(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; \[IndentingNewLine]While[
                  a[\([i]\)] === Table[0, {m}] && i > k, 
                  i = i - 1]; \[IndentingNewLine]If[i > k, 
                  For[j = k, j \[LessEqual] i, \(j++\), 
                    a = ReplacePart[a, a[\([j + 1]\)], 
                        j]; \[IndentingNewLine]u = 
                      ReplacePart[u, u[\([j + 1]\)], 
                        j]]];] (*endif*) ];  (*endfor*) If[
              ch, {a, u} = {a, u} /. {x^\((n : _)\) \[Rule] x^\(-n\), 
                    x \[Rule] x^\(-1\)}]; \[IndentingNewLine]{a, u}] /; 
          MatrixQ[at] || 
            Message[General::mtrx, "\<ExtendedHermiteForm\>"];\)\)], "Input",
  InitializationCell->True],

Cell["\<\
Same as above with the only difference that the unimodular matrix u \
is not computed\
\>", "SmallText"],

Cell[BoxData[
    \(\(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]; \[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]{p, m} = 
              Dimensions[a]; \[IndentingNewLine]For[k = 1; kcl = 1, 
              k \[LessEqual] p && kcl \[LessEqual] m, \(k++\), 
              While[a[\([Range[k, p], kcl]\)] === Table[0, {p - k + 1}] && 
                  kcl < m, \(kcl++\)]; \[IndentingNewLine]While[
                If[k \[Equal] p, False, 
                  a[\([Range[k + 1, p], kcl]\)] =!= Table[0, {p - k}]], 
                rmin = FindMin[\(Transpose[a]\)[\([kcl]\)], k, x, 
                    p]; \[IndentingNewLine]If[rmin \[NotEqual] k, 
                  a = a /. {a[\([rmin]\)] \[Rule] a[\([k]\)],