(******************************************************************* This file was generated automatically by the Mathematica front end. It contains Initialization cells from a Notebook file, which typically will have the same name as this file except ending in ".nb" instead of ".m". This file is intended to be loaded into the Mathematica kernel using the package loading commands Get or Needs. Doing so is equivalent to using the Evaluate Initialization Cells menu command in the front end. DO NOT EDIT THIS FILE. This entire file is regenerated automatically each time the parent Notebook file is saved in the Mathematica front end. Any changes you make to this file will be overwritten. ***********************************************************************) BeginPackage["PolynomialMatrices`"] 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]1,Message[General::badarg]; Return[$Failed],If[x==={},x=Variabile[lt]]]; (*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}]]]; 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[dega1,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:_))\[Rule](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\[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]]; {i,j}={i,j}+k-1; (*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]]}; u=u/.{u[[k]]\[Rule]u[[i]],u[[i]]\[Rule]u[[k]]}]; If[j\[NotEqual]k, a=Transpose[#/.{#[[k]]\[Rule]#[[j]],#[[j]]\[Rule]#[[k]]}&[ Transpose[a]]]; v=Transpose[#/.{#[[k]]\[Rule]#[[j]],#[[j]]\[Rule]#[[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]\[Equal]False, For[r=k+1,r\[LessEqual]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\[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]; 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\[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]; 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]]\[Rule]a[[k-1]],a[[k-1]]\[Rule]a[[k]]}; u=u/.{u[[k]]\[Rule]u[[k-1]],u[[k-1]]\[Rule]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\[Equal]p||k\[Equal]m,t=True,k++]]] (*endwhile*)]; If[ch,{a,u,v}={a,u,v}/.{x^(n:_)\[Rule]x^-n,x\[Rule]x^-1}]; Return[{a,u,v}]]/; MatrixQ[at]||Message[General::mtrx,"ExtendedSmithForm"]; 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:_))\[Rule](x^-n)}; ch=True]; {p,m}=Dimensions[a]; k=1;t=False; 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&&esp1,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:_))\[Rule](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\[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}]&&kclk,i=i-1]; If[i>k,For[j=k,j\[LessEqual]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:_)\[Rule]x^-n,x\[Rule]x^-1}]; {a,u}]/;MatrixQ[at]||Message[General::mtrx,"ExtendedHermiteForm"]; 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:_))\[Rule](x^-n)}; ch=True]; {p,m}=Dimensions[a]; For[k=1;kcl=1,k\[LessEqual]p&&kcl\[LessEqual]m,k++, While[a[[Range[k,p],kcl]]===Table[0,{p-k+1}]&&kclk,i=i-1]; 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}]; 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]]]; (*Compute the monic lcm of the denominators*)lcm= Apply[PolynomialLCM,Flatten[Denominator[Together[a]]]]; coef=Coefficient[lcm,x,Exponent[lcm,x]]; If[coef=!=1,lcm=lcm/coef]; (*Compute the Smith form s of the polynomial matrix a*lcm*){s,u,v}= ExtendedSmithForm[Apart[a*lcm],x]; (*The McMillan form is s/lcm*){Together[s/lcm],u,v}]/; MatrixQ[at]||Message[General::mtrx,"ExtendedMcMillanForm"]; 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"]; 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:_))\[Rule](x^-n)}; ch=True]; {p,m}=Dimensions[a]; u=IdentityMatrix[p]; For[k=kcl=1,k\[LessEqual]p&&kcl\[LessEqual]m,k++, While[a[[Range[k,p],kcl]]===Table[0,{p-k+1}]&&kcl(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[l1,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]]\[Equal]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]]\[Equal]Dimensions[bt][[1]]; 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:_))\[Rule](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\[NotEqual]m+1&&k\[NotEqual]l+p,k++,done=False; While[f[[Range[k,l+p],kcl]]===Table[0,{l+p-k+1}]&&kcl(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[m1,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]]\[Equal]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]]\[Equal]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:_))\[Rule](x^-n)};]; d=Det[RGCD[a,b]]; (FreeQ[d,x]&&d\[NotEqual]0)]/; MatrixQ[at]&&MatrixQ[bt]&& Dimensions[at][[2]]\[Equal]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:_))\[Rule](x^-n)};]; d=Det[LGCD[a,b]]; (FreeQ[d,x]&&d\[NotEqual]0)]/; MatrixQ[at]&&MatrixQ[bt]&& Dimensions[at][[1]]\[Equal]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:_))\[Rule](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]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:_))\[Rule](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]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:_))\[Rule](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:_)\[Rule]x^-n,x\[Rule]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"]; 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]]\[Equal]Dimensions[b][[1]]\[Equal] 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:_))\[Rule](x^-n)}; ch=True]; (*Ysol, Zsol are the polynomial matrix solutions of the equation a.Y+ b.Z\[Equal]c*)Ysol={{}}; Zsol={{}}; (*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={}; 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\[Equal]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[#\[Rule]0&,Complement[Variables[{Z,Y}],{x}]]*)]}, 1]]; Zsol=Transpose[Flatten[{Transpose[Zsol],Transpose[(Z/.sol)[[1]] (*/.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}]];]; If[ch,{Ysol,Zsol}={Ysol,Zsol}/.{x^(n:_)\[Rule]x^-n,x\[Rule]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]]]]; End[] EndPackage[]