(************** Content-type: application/mathematica ************** CreatedBy='Mathematica 5.2' Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 71201, 1496]*) (*NotebookOutlinePosition[ 71955, 1522]*) (* CellTagsIndexPosition[ 71911, 1518]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Polynomial Matrices", "Title"], Cell[CellGroupData[{ Cell["Author", "Subsubsection"], Cell["\<\ Adriano Pascoletti Dipartimento di Matematica e Informatica Universit\[AGrave] di Udine\ \>", "Text"], Cell["v. 2.1. december 1999", "SmallText"] }, Open ]], Cell[CellGroupData[{ Cell["Summary", "Subsubsection"], Cell[TextData[{ "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 \nthe form and the unimodular \ transformations leading to it. \nSeveral other functions are defined: \ computation of left and right GCD\[CloseCurlyQuote]s, \ lcm\[CloseCurlyQuote]s, quotients, remainders, row/column proper forms, full \ rank and coprimality tests, extended polynomial GCD, solution of diophantine \ equations. \n", "All functions work on matrices of polynomials in z or ", Cell[BoxData[ \(TraditionalForm\`z\^\(-1\)\)]], ".", "\n" }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Notebook version", "Subsubsection"], Cell["2.1", "Text"] }, Open ]], Cell[CellGroupData[{ Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], " Version" }], "Subsubsection"], Cell["3.x and 4.0", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["History", "Subsubsection"], Cell[TextData[{ "This package supersedes package 0207-751. \nVersion 2.1 removes some \ incompatibilities between ", StyleBox["Mathematica", FontSlant->"Italic"], " 2.2 and ", StyleBox["Mathematica", FontSlant->"Italic"], " 3.x and 4.0.\nVersion 2.0 is the result of a work with M. Anziutti." }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Keywords", "Subsubsection"], Cell["\<\ Polynomials, polynomial matrices, non commutative algebra, Smith \ form, Hermite Form, McMillan Form, diophantine equations, coprimality, \ properness, greatest common divisors, least common multiples.\ \>", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Interface", "Section"], Cell[CellGroupData[{ Cell["Package context", "Subsubsection"], Cell["BeginPackage[\"PolynomialMatrices`\"]", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False] }, Open ]], Cell[CellGroupData[{ Cell["Usage messages", "Subsubsection"], Cell[BoxData[ \(\(DiophantineSolve::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(EliminateDep::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(ExtendedHermiteForm::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(ExtendedMcMillanForm::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(ExtendedSmithForm::usage = "\";\)\)], \ "Input", InitializationCell->True], Cell[BoxData[ \(\(ExtPolQ::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(FullRowRankQ::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(HermiteForm::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LCoprime::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LDivision::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LGCD::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(Llcm::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LQuotient::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LRemainder::usage = "\";\)\)], \ "Input", InitializationCell->True], Cell[BoxData[ \(\(McMillanForm::usage = "\< McMillanForm[a,x] takes a matrix of \ rational forms in x and returns its McMillan form.\>";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(PolyExtendedGCD::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(PolynomialColumnReduce::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(PolynomialDegree::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(PolynomialRowReduce::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(Rank::usage = "\";\)\)], \ "Input", InitializationCell->True], Cell[BoxData[ \(\(RCoprime::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(RDivision::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(RGCD::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[{ \(\(\(Rlcm::usage = "\";\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(RQuotient::usage = "\";\)\)}], "Input", InitializationCell->True], Cell[BoxData[ \(\(RRemainder::usage = "\";\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(SmithForm::usage = "\";\)\)], "Input", InitializationCell->True] }, Open ]], Cell[CellGroupData[{ Cell["Error messages", "Subsubsection"], Cell[BoxData[{ \(\(\(General::mtrx = "\";\)\(\ \[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(General::pol = "\";\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(PolyExtendedGCD::egcdz = "\";\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(General::badarg = "\";\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(General::gcd = "\";\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(LDivision::badarg = "\";\)\(\[IndentingNewLine]\ \) \)\), "\[IndentingNewLine]", \(\(\(PolynomialColumnReduce::badarg = "\";\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(PolynomialRowReduce::badarg = "\";\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(DiophantineSolve::badarg = "\"\)\(\[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], 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, "\"];\)\)], "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, "\"];\)\)], "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, "\"];\)\)], "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]\)], a[\([k]\)] \[Rule] a[\([rmin]\)]}]; \[IndentingNewLine]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]]]]; (*Endwhile*) \[IndentingNewLine]If[ a[\([k, kcl]\)] =!= 0, For[i = 1, i < k, \(i++\), If[a[\([i, kcl]\)] =!= 0, While[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];] (*endwhile*) ]]]; \[IndentingNewLine]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]\(kcl++\)]; (*Endfor*) 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]]];] (*endif*) ]; (*endfor*) If[ch, a = a /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]Return[a]] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(ExtendedMcMillanForm[at_, xt_: {}] := Module[{a = at, x = xt, s, u, v, lcm, coef}, If[x === {} && Length[Variables[a]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[ a]]]; \[IndentingNewLine] (*Compute\ the\ monic\ lcm\ of\ \ the\ denominators*) lcm = Apply[PolynomialLCM, Flatten[Denominator[Together[a]]]]; \[IndentingNewLine]coef = Coefficient[lcm, x, Exponent[lcm, x]]; \[IndentingNewLine]If[ coef =!= 1, lcm = lcm/ coef]; \[IndentingNewLine] (*Compute\ the\ Smith\ form\ s\ \ of\ the\ polynomial\ matrix\ a*lcm*) {s, u, v} = ExtendedSmithForm[Apart[a*lcm], x]; \[IndentingNewLine] (*The\ McMillan\ form\ is\ s/ lcm*) {Together[s/lcm], u, v}] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], "Input",\ InitializationCell->True], Cell["Same as above without unimodular matrices", "SmallText"], Cell[BoxData[ \(\(McMillanForm[at_, xt_: {}] := Module[{a = at, x = xt, s, lcm, coef}, If[x === {} && Length[Variables[a]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[a]]]; \[IndentingNewLine]lcm = Apply[PolynomialLCM, Flatten[Denominator[Together[a]]]]; \[IndentingNewLine]coef = Coefficient[lcm, x, Exponent[lcm, x]]; \[IndentingNewLine]If[ coef =!= 1, lcm = lcm/coef]; \[IndentingNewLine]s = SmithForm[Apart[a*lcm], x]; \[IndentingNewLine]Together[ s/lcm]] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], "Input", InitializationCell->True], Cell["\<\ The following auxiliary function eliminates the dependent rows of \ the polynomial matrix at\ \>", "SmallText"], Cell[BoxData[ \(\(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]; \[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]u = IdentityMatrix[p]; \[IndentingNewLine]For[k = \(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]\)], a[\([k]\)] \[Rule] a[\([rmin]\)]}; \[IndentingNewLine]u = u /. {u[\([rmin]\)] \[Rule] u[\([k]\)], u[\([k]\)] \[Rule] u[\([rmin]\)]}]; \[IndentingNewLine]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, Expand[a[\([i]\)] - a[\([k]\)]*lpc, x], i];]]]; (*Endwhile*) \(kcl++\)]; (*endfor*) k = Flatten[Position[Inverse[u] . a, n_ /; n === Table[0, {m}]]]; \[IndentingNewLine]If[ch, b = b /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]Return[\(IdentityMatrix[ p]\)[\([Complement[Range[1, p], k]]\)] . b]] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], "Input", InitializationCell->True], Cell["\<\ The following auxiliary function computes a left GCD and a right \ lcm of matrices at and bt\ \>", "SmallText"], Cell[BoxData[ \(\(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]]; \[IndentingNewLine]If[ Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a, b} = {a, b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \ \[IndentingNewLine]ch = True]; \[IndentingNewLine] (*Initialize\ the\ variables\ *) \[IndentingNewLine]{l, m} = Dimensions[ a]; \[IndentingNewLine]q = \(Dimensions[ b]\)[\([2]\)]; \[IndentingNewLine]v = IdentityMatrix[m + q]; \[IndentingNewLine]f = Transpose[ Flatten[{Transpose[a], Transpose[b]}, 1]]; \[IndentingNewLine] (*Main\ loop\ on\ k, columns\ of\ f*) \[IndentingNewLine]For[k = \(krow = 1\), krow \[NotEqual] l + 1 && k \[NotEqual] m + q, \(k++\), done = False; \[IndentingNewLine]While[ f[\([krow, Range[k, m + q]]\)] === Table[0, {m + q - k + 1}] && krow < l, \(krow++\)]; \[IndentingNewLine] (*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]; \[IndentingNewLine]If[j \[NotEqual] k, f = Transpose[\(# /. {#[\([k]\)] \[Rule] #[\([j]\)], \ #[\([j]\)] \[Rule] #[\([k]\)]} &\)[Transpose[f]]]; \[IndentingNewLine]v = Transpose[\(# /. {#[\([k]\)] \[Rule] #[\([j]\)], \ #[\([j]\)] \[Rule] #[\([k]\)]} &\)[ Transpose[ v]]]]; \[IndentingNewLine] (*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 \[LessEqual] m + q, \(i++\), If[f[\([krow, i]\)] =!= 0, lambda = PolynomialQuotient[f[\([krow, i]\)], f[\([krow, k]\)], x]; \[IndentingNewLine]f = Transpose[\(ReplacePart[#, Expand[#[\([i]\)] - #[\([k]\)] \((lambda)\), x], i] &\)[ Transpose[f]]]; \[IndentingNewLine]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; \[IndentingNewLine]While[ j \[GreaterEqual] 1 && \(Transpose[f]\)[\([j]\)] === Table[0, {l}], j = j - 1]; \[IndentingNewLine] (*Return\ the\ left\ GCD\ and\ if\ \ Flatten[lambda] =!= {}\ return\ also\ a\ right\ lcm*) If[ ch, {a, v, f} = {a, v, f} /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]lambda = v[\([Range[1, m], Range[j + 1, q + m]]\)]; \[IndentingNewLine]{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], {}], {}]}];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LGCD[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, If[x === {} && Length[Variables[{a, b}]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]\(If[# =!= $Failed, \ #[\([1]\)], #] &\)[AuxLGCDlcm[a, b, x]]] /; MatrixQ[at] && MatrixQ[bt] && \(Dimensions[at]\)[\([1]\)] \[Equal] \(Dimensions[ bt]\)[\([1]\)];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(Rlcm[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, If[x === {} && Length[Variables[{a, b}]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]\(If[# =!= $Failed, \ #[\([2]\)], #] &\)[AuxLGCDlcm[a, b, x]]] /; MatrixQ[at] && MatrixQ[bt] && \(Dimensions[at]\)[\([1]\)] \[Equal] \(Dimensions[ bt]\)[\([1]\)];\)\)], "Input", InitializationCell->True], Cell["\<\ The following auxiliary function computes a right GCD and a left \ lcm of matrices at and bt\ \>", "SmallText"], Cell[BoxData[ \(\(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]]; \[IndentingNewLine]If[ Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a, b} = {a, b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \ \[IndentingNewLine]ch = True]; \[IndentingNewLine] (*Initialize\ the\ variables*) {l, m} = Dimensions[ a]; \[IndentingNewLine]p = \(Dimensions[ b]\)[\([1]\)]; \[IndentingNewLine]v = IdentityMatrix[l + p]; \[IndentingNewLine]f = Flatten[{a, b}, 1]; \[IndentingNewLine] (*Main\ loop\ on\ k, rows\ of\ f*) For[k = \(kcl = 1\), kcl \[NotEqual] m + 1 && k \[NotEqual] l + p, \(k++\), done = False; \[IndentingNewLine]While[ f[\([Range[k, l + p], kcl]\)] === Table[0, {l + p - k + 1}] && kcl < m, \(kcl++\)]; \[IndentingNewLine] (*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]; \[IndentingNewLine]If[j \[NotEqual] k, f = f /. {f[\([k]\)] \[Rule] f[\([j]\)], f[\([j]\)] \[Rule] f[\([k]\)]}; \[IndentingNewLine]v = v /. {v[\([k]\)] \[Rule] v[\([j]\)], v[\([j]\)] \[Rule] v[\([k]\)]}]; \[IndentingNewLine] (*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 \[LessEqual] l + p, \(i++\), If[f[\([i, kcl]\)] =!= 0, lambda = PolynomialQuotient[f[\([i, kcl]\)], f[\([k, kcl]\)], x]; \[IndentingNewLine]f = ReplacePart[f, Expand[f[\([i]\)] - f[\([k]\)] \((lambda)\), x], i]; \[IndentingNewLine]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; \[IndentingNewLine]While[ j \[GreaterEqual] 1 && f[\([j]\)] === Table[0, {m}], j = j - 1]; \[IndentingNewLine]If[ ch, {a, v, f} = {a, v, f} /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]lambda = v[\([Range[j + 1, l + p], Range[1, l]]\)]; \[IndentingNewLine] (*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], {}], {}]}];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(\(RGCD[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, If[x === {} && Length[Variables[{a, b}]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]\(If[# =!= $Failed, \ #[\([1]\)], #] &\)[AuxRGCDlcm[a, b, x]]] /; MatrixQ[at] && MatrixQ[bt] && \(Dimensions[at]\)[\([2]\)] \[Equal] \(Dimensions[ bt]\)[\([2]\)];\)\(\ \)\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(Llcm[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, If[x === {} && Length[Variables[{a, b}]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]\(If[# =!= $Failed, \ #[\([2]\)], #] &\)[AuxRGCDlcm[a, b, x]]] /; MatrixQ[at] && MatrixQ[bt] && \(Dimensions[at]\)[\([2]\)] \[Equal] \(Dimensions[ bt]\)[\([2]\)];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(RCoprime[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt, d}, If[x === {} && Length[Variables[{a, b}]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]If[ Not[ExtPolQ[{a, b}, x]], Message[General::pol]; Return[$Failed]]; \[IndentingNewLine]If[ Not[Apply[ And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], \({a, b} = {a, b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)};\)]; \ \[IndentingNewLine]d = Det[RGCD[a, b]]; \[IndentingNewLine]\((FreeQ[d, x] && d \[NotEqual] 0)\)] /; MatrixQ[at] && MatrixQ[bt] && \(Dimensions[at]\)[\([2]\)] \[Equal] \(Dimensions[ bt]\)[\([2]\)];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LCoprime[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt, d}, If[x === {} && Length[Variables[{a, b}]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]If[ Not[ExtPolQ[{a, b}, x]], Message[General::pol]; Return[$Failed]]; \[IndentingNewLine]If[ Not[Apply[ And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], \({a, b} = {a, b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)};\)]; \ \[IndentingNewLine]d = Det[LGCD[a, b]]; \[IndentingNewLine]\((FreeQ[d, x] && d \[NotEqual] 0)\)] /; MatrixQ[at] && MatrixQ[bt] && \(Dimensions[at]\)[\([1]\)] \[Equal] \(Dimensions[ bt]\)[\([1]\)];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(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]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]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, "\", "\"]; Return[$Failed]]; \[IndentingNewLine]If[Not[ExtPolQ[{a, b}, x]], Message[General::pol]; Return[$Failed]]; \[IndentingNewLine]If[ Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a, b} = {a, b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \ \[IndentingNewLine]ch = True]; \[IndentingNewLine] (*Inizialization : delta\ is\ the\ inverse\ of\ the\ leading\ coefficient\ matrix\ \ of\ matrix\ b*) \[IndentingNewLine]{l, m} = Dimensions[a]; \[IndentingNewLine]u = Table[0, {l}, {m}]; \[IndentingNewLine]delta = Inverse[Coefficient[b, x, Max[Exponent[b, x]]]]; \[IndentingNewLine] (*Main\ \(loop : go\ on\ until\ deg[a] < deg[b]\)*) While[ Max[Exponent[a, x]] \[GreaterEqual] Max[Exponent[b, x]], lambda = delta . Coefficient[a, x, Max[Exponent[a, x]]]; \[IndentingNewLine]esp = Max[Exponent[a, x]] - Max[Exponent[b, x]]; \[IndentingNewLine] (*Decrease\ the\ degree\ of\ \ matrix\ a*) a = Expand[a - b . \((lambda*x^esp)\), x]; \[IndentingNewLine]u = Expand[u + lambda*x^esp, x]]; \[IndentingNewLine]If[ ch, {u, a} = {u, a} /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]{u, a}];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LQuotient[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, \(If[# =!= $Failed, #[\([1]\)], #] &\)[ LDivision[a, b, x]]];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(LRemainder[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, \(If[# =!= $Failed, #[\([2]\)], #] &\)[ LDivision[a, b, x]]];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(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]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[{a, b}]]]; \[IndentingNewLine]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, "\", "\"]; Return[$Failed]]; \[IndentingNewLine]If[Not[ExtPolQ[{a, b}, x]], Message[General::pol]; Return[$Failed]]; \[IndentingNewLine]If[ Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a, b} = {a, b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \ \[IndentingNewLine]ch = True]; \[IndentingNewLine] (*Some\ inizialization, delta\ is\ the\ inverse\ of\ the\ leading\ coefficient\ matrix\ \ of\ matrix\ b*) {l, m} = Dimensions[a]; \[IndentingNewLine]u = Table[0, {l}, {m}]; \[IndentingNewLine]delta = Inverse[Coefficient[b, x, Max[Exponent[b, x]]]]; \[IndentingNewLine] (*Main\ \(loop : go\ on\ until\ deg[a] < deg[b]\)*) While[ Max[Exponent[a, x]] \[GreaterEqual] Max[Exponent[b, x]], lambda = Coefficient[a, x, Max[Exponent[a, x]]] . delta; \[IndentingNewLine]esp = Max[Exponent[a, x]] - Max[Exponent[b, x]]; \[IndentingNewLine] (*Decrease\ the\ degree\ of\ \ matrix\ a*) a = Expand[a - \((lambda*x^esp)\) . b, x]; \[IndentingNewLine]u = Expand[u + lambda*x^esp, x]]; \[IndentingNewLine]If[ ch, {u, a} = {u, a} /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]{u, a}];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(RQuotient[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, \(If[# =!= $Failed, #[\([1]\)], #] &\)[ RDivision[a, b, x]]];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(RRemainder[at_, bt_, xt_: {}] := Module[{a = at, b = bt, x = xt}, \(If[# =!= $Failed, #[\([2]\)], #] &\)[ RDivision[a, b, x]]];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(Rank[at_] := Module[{a = at, x = xt, p, m}, {p, m} = Dimensions[a]; \[IndentingNewLine]a = RowReduce[a]; \[IndentingNewLine]Return[ p - Count[a, n_ /; n === Table[0, {m}]]]] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(FullRowRankQ[at_] := Module[{p = \(Dimensions[at]\)[\([1]\)]}, Return[Rank[at] === p]] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(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]; \[IndentingNewLine]Return[$Failed], 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]If[\(! FullRowRankQ[ Transpose[a]]\), Message[PolynomialColumnReduce::badarg]; Return[$Failed]]; \[IndentingNewLine]{p, m} = Dimensions[a]; \[IndentingNewLine]u = IdentityMatrix[ m]; \[IndentingNewLine] (*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]]]; \[IndentingNewLine] (*go\ on\ until\ Ahc\ has\ full\ \ column\ rank\ \((i . e . a\ is\ column\ reduced)\)*) While[\(! FullRowRankQ[Transpose[Ahc]]\), v = \(NullSpace[Ahc]\)[\([1]\)]; \[IndentingNewLine]deg = Table[Max[Exponent[\(Transpose[a]\)[\([i]\)], x]], {i, 1, m}]; \[IndentingNewLine]k = \(Flatten[ Position[deg, degmax = Max[Table[ Exponent[\(Transpose[a]\)[\([i]\)]*v[\([i]\)], x], {i, 1, m}]]]]\)[\([1]\)]; \[IndentingNewLine] (*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]]]; \[IndentingNewLine]Ahc = Transpose[\(Table[ Coefficient[#[\([i]\)], x, Max[Exponent[#[\([i]\)], x]]], {i, 1, m}] &\)[ Transpose[a]]]]; \[IndentingNewLine]If[ch, u = u /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]Return[u]] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], \ "Input", InitializationCell->True], Cell[BoxData[ \(\(PolynomialRowReduce[at_, xt_: {}] := Module[{a = at, x = xt}, If[x === {}, If[Length[Variables[a]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], x = Variabile[a]]]; \[IndentingNewLine]If[Not[ExtPolQ[a, x]], Message[General::pol]; Return[$Failed]]; \[IndentingNewLine]If[\(! FullRowRankQ[a]\), Message[PolynomialRowReduce::badarg]; \ \[IndentingNewLine]Return[$Failed]]; \[IndentingNewLine]Return[ Transpose[PolynomialColumnReduce[Transpose[a], x]]]] /; MatrixQ[at] || Message[General::mtrx, "\"];\)\)], "Input", InitializationCell->True], Cell["Test whether a diophantine equation has a solution", "SmallText"], Cell[BoxData[ \(\(EquivQ[R1t_, R2t_, xt_: {}] := Module[{R1 = R1t, R2 = R2t, x = xt, h1, h2}, h1 = HermiteForm[R1, x]; \[IndentingNewLine]h2 = HermiteForm[R2, x]; \[IndentingNewLine]If[h1 === h2, Return[True], Return[False]]];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(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]\)] \[Equal] \(Dimensions[ b]\)[\([1]\)] \[Equal] \(Dimensions[c]\)[\([1]\)]], Message[DiophantineSolve::badarg]; Return[$Failed]]; \[IndentingNewLine]If[x === {}, If[Length[Variables[{a, b, c}]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], x = Variabile[{a, b, c}]]]; \[IndentingNewLine]If[ Not[ExtPolQ[{a, b, c}, x]], Message[General::pol]; Return[$Failed]]; \[IndentingNewLine]If[ Not[Apply[ And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b, c}]]], {a, b, c} = {a, b, c} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \ \[IndentingNewLine]ch = True]; \[IndentingNewLine] (*Ysol, Zsol\ are\ the\ polynomial\ matrix\ solutions\ of\ the\ equation\ \ a . Y + b . Z \[Equal] c*) Ysol = {{}}; \[IndentingNewLine]Zsol = {{}}; \[IndentingNewLine] \ (*Test\ wether\ the\ equation\ has\ solutions\ or\ not*) If[ Rank[Transpose[ Flatten[{Transpose[a], Transpose[b]}, 1]]] \[Equal] 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 \[LessEqual] \(Dimensions[c]\)[\([2]\)], \(t++\), sol = {}; \[IndentingNewLine]deg = 0; \[IndentingNewLine] (*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}]; \[IndentingNewLine]Z = Table[Sum[ z[i, k]\ x^k, {k, 0, deg}], {i, \(Dimensions[ b]\)[\([2]\)]}, {1}]; \[IndentingNewLine] (*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 \[Equal] \(Transpose[c]\)[\([t]\)], Variables[{a, b, c}]]; \[IndentingNewLine]\(deg++\)]; (*endwhile*) \ (*Update\ the\ solutions\ with\ the\ new\ minimal\ degree\ columns\ just\ \ found*) Ysol = Transpose[ Flatten[{Transpose[Ysol], Transpose[\((Y /. sol)\)[\([1]\)]\[IndentingNewLine] \ (*\(/.\)\(Map[# \[Rule] 0 &, Complement[Variables[{Z, Y}], {x}]]\)*) ]}, 1]]; \[IndentingNewLine]Zsol = Transpose[ Flatten[{Transpose[Zsol], Transpose[\((Z /. sol)\)[\([1]\)]\[IndentingNewLine] \ (*\(/.\)\(Map[# \[Rule] 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 \[LessEqual] Length[#], \(i++\), Module[{t = dummy}, {Ysol, Zsol} = {Ysol, Zsol} /. #[\([i]\)] \[Rule] t[i]]] &\)[ Complement[Variables[{Ysol, Zsol}], Variables[{a, b, c}]]], \({Ysol, Zsol} = {Ysol, Zsol} /. Map[# \[Rule] 0 &, Complement[ Variables[{Ysol, Zsol}], {x}]];\)]; \[IndentingNewLine]If[ ch, {Ysol, Zsol} = {Ysol, Zsol} /. {x^\((n : _)\) \[Rule] x^\(-n\), x \[Rule] x^\(-1\)}]; \[IndentingNewLine]Return[{Ysol, Zsol}]];\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(PolynomialDegree[at_List, xt_: {}] := Module[{x = xt}, If[x === {} && Length[Variables[at]] > 1, Message[General::badarg]; \[IndentingNewLine]Return[$Failed], If[x === {}, x = Variabile[at]]]; \[IndentingNewLine]Return[ Max[Exponent[at, x]]]];\)\)], "Input", InitializationCell->True], Cell[CellGroupData[{ Cell["Epilog", "Section"], Cell[BoxData[ \(\(\(\[IndentingNewLine]\)\(End[]\)\)\)], "Input", InitializationCell->True], Cell[BoxData[ \(\(\(\[IndentingNewLine]\)\(EndPackage[]\)\)\)], "Input", InitializationCell->True] }, Open ]] }, Open ]] }, Open ]] }, FrontEndVersion->"5.2 for Macintosh", ScreenRectangle->{{0, 1400}, {0, 878}}, AutoGeneratedPackage->Automatic, WindowToolbars->"RulerBar", PageWidth->PaperWidth, WindowSize->{758, 555}, WindowMargins->{{-2, Automatic}, {Automatic, 1}}, StyleDefinitions -> "Default.nb" ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1776, 53, 36, 0, 88, "Title"], Cell[CellGroupData[{ Cell[1837, 57, 31, 0, 28, "Subsubsection"], Cell[1871, 59, 111, 4, 62, "Text"], Cell[1985, 65, 42, 0, 26, "SmallText"] }, Open ]], Cell[CellGroupData[{ Cell[2064, 70, 32, 0, 28, "Subsubsection"], Cell[2099, 72, 850, 16, 190, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[2986, 93, 41, 0, 28, "Subsubsection"], Cell[3030, 95, 19, 0, 30, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[3086, 100, 101, 4, 28, "Subsubsection"], Cell[3190, 106, 27, 0, 30, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[3254, 111, 32, 0, 28, "Subsubsection"], Cell[3289, 113, 323, 9, 82, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[3649, 127, 33, 0, 28, "Subsubsection"], Cell[3685, 129, 225, 4, 62, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[3947, 138, 28, 0, 69, "Section"], Cell[CellGroupData[{ Cell[4000, 142, 40, 0, 28, "Subsubsection"], Cell[4043, 144, 137, 3, 27, "Input", InitializationCell->True] }, Open ]], Cell[CellGroupData[{ Cell[4217, 152, 39, 0, 28, "Subsubsection"], Cell[4259, 154, 527, 8, 219, "Input", InitializationCell->True], Cell[4789, 164, 222, 4, 91, "Input", InitializationCell->True], Cell[5014, 170, 239, 4, 91, "Input", InitializationCell->True], Cell[5256, 176, 266, 4, 107, "Input", InitializationCell->True], Cell[5525, 182, 280, 5, 107, "Input", InitializationCell->True], Cell[5808, 189, 188, 3, 75, "Input", InitializationCell->True], Cell[5999, 194, 174, 3, 75, "Input", InitializationCell->True], Cell[6176, 199, 170, 3, 59, "Input", InitializationCell->True], Cell[6349, 204, 221, 4, 91, "Input", InitializationCell->True], Cell[6573, 210, 260, 4, 107, "Input", InitializationCell->True], Cell[6836, 216, 228, 4, 91, "Input", InitializationCell->True], Cell[7067, 222, 244, 4, 107, "Input", InitializationCell->True], Cell[7314, 228, 265, 4, 107, "Input", InitializationCell->True], Cell[7582, 234, 283, 5, 123, "Input", InitializationCell->True], Cell[7868, 241, 184, 3, 75, "Input", InitializationCell->True], Cell[8055, 246, 259, 4, 107, "Input", InitializationCell->True], Cell[8317, 252, 272, 4, 107, "Input", InitializationCell->True], Cell[8592, 258, 175, 3, 75, "Input", InitializationCell->True], Cell[8770, 263, 252, 4, 107, "Input", InitializationCell->True], Cell[9025, 269, 127, 3, 43, "Input", InitializationCell->True], Cell[9155, 274, 212, 4, 91, "Input", InitializationCell->True], Cell[9370, 280, 294, 5, 123, "Input", InitializationCell->True], Cell[9667, 287, 213, 4, 91, "Input", InitializationCell->True], Cell[9883, 293, 483, 8, 203, "Input", InitializationCell->True], Cell[10369, 303, 271, 4, 123, "Input", InitializationCell->True], Cell[10643, 309, 166, 3, 59, "Input", InitializationCell->True] }, Open ]], Cell[CellGroupData[{ Cell[10846, 317, 39, 0, 28, "Subsubsection"], Cell[10888, 319, 1617, 30, 651, "Input", InitializationCell->True] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[12554, 355, 33, 0, 69, "Section"], Cell[12590, 357, 83, 2, 27, "Input", InitializationCell->True] }, Open ]], Cell[CellGroupData[{ Cell[12710, 364, 24, 0, 52, "Subtitle"], Cell[12737, 366, 294, 5, 54, "SmallText"], Cell[13034, 373, 206, 5, 43, "Input", InitializationCell->True], Cell[13243, 380, 270, 5, 54, "SmallText"], Cell[13516, 387, 933, 18, 171, "Input", InitializationCell->True], Cell[14452, 407, 160, 4, 26, "SmallText"], Cell[14615, 413, 1615, 27, 331, "Input", InitializationCell->True], Cell[16233, 442, 113, 3, 26, "SmallText"], Cell[16349, 447, 176, 4, 27, "Input", InitializationCell->True], Cell[16528, 453, 8042, 133, 1163, "Input", InitializationCell->True], Cell[24573, 588, 136, 3, 26, "SmallText"], Cell[24712, 593, 5532, 93, 795, "Input", InitializationCell->True], Cell[30247, 688, 103, 3, 26, "SmallText"], Cell[30353, 693, 474, 9, 107, "Input", InitializationCell->True], Cell[30830, 704, 4832, 81, 715, "Input", InitializationCell->True], Cell[35665, 787, 114, 3, 70, "SmallText"], Cell[35782, 792, 3310, 56, 70, "Input", InitializationCell->True], Cell[39095, 850, 1088, 22, 70, "Input", InitializationCell->True], Cell[40186, 874, 62, 0, 70, "SmallText"], Cell[40251, 876, 784, 14, 70, "Input", InitializationCell->True], Cell[41038, 892, 121, 3, 70, "SmallText"], Cell[41162, 897, 2439, 40, 70, "Input", InitializationCell->True], Cell[43604, 939, 121, 3, 70, "SmallText"], Cell[43728, 944, 4098, 72, 70, "Input", InitializationCell->True], Cell[47829, 1018, 586, 12, 70, "Input", InitializationCell->True], Cell[48418, 1032, 586, 12, 70, "Input", InitializationCell->True], Cell[49007, 1046, 121, 3, 70, "SmallText"], Cell[49131, 1051, 3774, 66, 70, "Input", InitializationCell->True], Cell[52908, 1119, 596, 12, 70, "Input", InitializationCell->True], Cell[53507, 1133, 586, 12, 70, "Input", InitializationCell->True], Cell[54096, 1147, 952, 18, 70, "Input", InitializationCell->True], Cell[55051, 1167, 952, 18, 70, "Input", InitializationCell->True], Cell[56006, 1187, 2191, 38, 70, "Input", InitializationCell->True], Cell[58200, 1227, 225, 5, 70, "Input", InitializationCell->True], Cell[58428, 1234, 226, 5, 70, "Input", InitializationCell->True], Cell[58657, 1241, 2162, 37, 70, "Input", InitializationCell->True], Cell[60822, 1280, 225, 5, 70, "Input", InitializationCell->True], Cell[61050, 1287, 226, 5, 70, "Input", InitializationCell->True], Cell[61279, 1294, 356, 7, 70, "Input", InitializationCell->True], Cell[61638, 1303, 246, 5, 70, "Input", InitializationCell->True], Cell[61887, 1310, 2718, 48, 70, "Input", InitializationCell->True], Cell[64608, 1360, 758, 14, 70, "Input", InitializationCell->True], Cell[65369, 1376, 71, 0, 70, "SmallText"], Cell[65443, 1378, 319, 6, 70, "Input", InitializationCell->True], Cell[65765, 1386, 4770, 84, 70, "Input", InitializationCell->True], Cell[70538, 1472, 366, 7, 70, "Input", InitializationCell->True], Cell[CellGroupData[{ Cell[70929, 1483, 25, 0, 70, "Section"], Cell[70957, 1485, 97, 2, 70, "Input", InitializationCell->True], Cell[71057, 1489, 104, 2, 70, "Input", InitializationCell->True] }, Open ]] }, Open ]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)