(* ::Package:: *)

(************************************************************************)
(* 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.                                                         *)
(************************************************************************)



(*:Mathematica Version:6.0*)


(*:Package Version:1.1*)


(*:Name:*)


(*:Author:Adriano Pascoletti*)


(*:URL:*)


(*:Summary:*)


(*:History:
v1.0 (2005-11-22):written by Adriano Pascoletti for Mathematica 5.2 
v1.1 (2007-11-16):made installable under Mathematica 6.0
*)


(*:Keywords:polynomial matrices, integer matrices, invariant factors, Smith form*)


(*:Requirements:None.*)


(*:References:*)


BeginPackage["SmithFormV6`"]


IntegerSmithForm::usage="IntegerSmithForm[A, ExtendedForm \[Rule] False] computes the  Smith form of an integer matrix A.\nIntegerSmithForm[A] and IntegerSmithForm[A, ExtendedForm \[Rule] True] yield {s,{u,v}} where s is the Smith form of A, and u, v are integer matrices satisfying u.A.v=s.";


IntegerInvariantFactors::usage="IntegerInvariantFactors[A] yields the list of invariant factors of an integer matrix A.";


PolynomialSmithForm::usage="PolynomialSmithForm[A, x, ExtendedForm \[Rule] False] computes the  Smith form of a matrix A of polynomials in x.\nPolynomialSmithForm[A, x] and PolynomialSmithForm[A, x, ExtendedForm \[Rule] True] yield {s,{u,v}} where s is the Smith form of A, and u,v are polynomial matrices satisfying u.A.v=s.\nThe second argument must be a symbol.";


PolynomialInvariantFactors::usage="PolynomialInvariantFactors[A,x] yields the  invariant factors of a matrix A of polynomials in x (x must be a symbol).";


ExtendedForm::usage="ExtendedForm \[Rule] True|False is an option to IntegerSmithForm and PolynomialSmithForm that enables|disables the extended Smith form.";


Begin["`PrivateP`"]


Options[PolynomialSmithForm]={ExtendedForm->True};


PolynomialSmithForm::"mat"="The first argument to `1` must be a nonempty matrix of polynomials in `2`.";


PolynomialSmithForm::"var"="The second argument to `1` must be a symbol.";


PolynomialSmithForm::"ierr"="Unespected internal error in `1`. Please contact pascolet@dimi.uniud.it .";


PolynomialSmithForm[at_,xt_,opts___?OptionQ]:=(If[Head[xt]=!=Symbol,Message[PolynomialSmithForm::"var",PolynomialSmithForm];Abort[]];
If[\[Not](Min[Dimensions[at]]>0\[And]MatrixQ[at,PolynomialQ[#,xt]&]),Message[PolynomialSmithForm::"mat",PolynomialSmithForm,xt];Abort[]];
If[ExtendedForm/.{opts}/.Options[PolynomialSmithForm],PolynomialSmithForm0[at,xt],PolynomialSmithForm1[at,xt]]);


PolynomialInvariantFactors[at_,xt_]:=(If[Head[xt]=!=Symbol,Message[PolynomialSmithForm::"var",PolynomialInvariantFactors];Abort[]];If[\[Not](Min[Dimensions[at]]>0\[And]MatrixQ[at,PolynomialQ[#,xt]&]),Message[PolynomialSmithForm::"mat",PolynomialInvariantFactors,xt];Abort[]];
DeleteCases[Tr[#,List]&@PolynomialSmithForm1[at,xt],0]);


PolynomialSmithForm0[at_,xt_]:=Block[{a=at,x=xt,u,v,p,m,r,i,L},{p,m}=Dimensions[a];
u=IdentityMatrix[p];v=IdentityMatrix[m];
i=1;
While[(i<=Min[p,m])\[And](\[Not]MatchQ[a[[Range[i,p],Range[i,m]]],{{(0)..}..}|{}]),Which[i<Min[p,m],While[\[Not](MatchQ[a[[Range[i+1,p],i]],{(0)..}|{}]\[And]MatchQ[a[[i,Range[i+1,m]]],{(0)..}|{}]),rowops[i];colops[i]];++i,i===p===m,++i,
i===p,If[i<m,colops[i]];++i,
i===m,If[i<p,rowops[i]];++i,True,Message[PolynomialSmithForm::"ierr",PolynomialSmithForm0];Abort[]]];
r=--i;
If[r>0,
makemonic/@Range[r];L=Tr[a,List];
With[{o=Ordering[Abs[Exponent[#,x]]&/@L],rL=Range[Length[L]]},
u[[rL]]=u[[o]];
v[[All,rL]]=v[[All,o]];
a[[rL,rL]]=a[[o,o]];L=L[[o]]];
r=Length[L=DeleteCases[L,0]];
i=1;
While[i<r,
If[PolynomialRemainder[L[[i+1]],L[[i]],x]===0,++i,finalstep[i];i=Max[--i,1]]
];
(a[[#,#]]=L[[#]])&/@Range[r]];
{a,{u,v}}];


PolynomialSmithForm1[at_,xt_]:=Block[{a=at,x=xt,p,m,r,i,L},{p,m}=Dimensions[a];i=1;
While[(i<=Min[p,m])\[And](\[Not]MatchQ[a[[Range[i,p],Range[i,m]]],{{(0)..}..}|{}]),Which[i<Min[p,m],While[\[Not](MatchQ[a[[Range[i+1,p],i]],{(0)..}|{}]\[And]MatchQ[a[[i,Range[i+1,m]]],{(0)..}|{}]),rowops1[i];colops1[i]];++i,i===p===m,++i,
i===p,If[i<m,colops1[i]];++i,
i===m,If[i<p,rowops1[i]];++i,True,Message[PolynomialSmithForm::"ierr",PolynomialSmithForm1];Abort[]]];
r=--i;
If[r>0,
makemonic1/@Range[r];L=Tr[a,List];
With[{o=Ordering[Abs[Exponent[#,x]]&/@L],rL=Range[Length[L]]},
a[[rL,rL]]=a[[o,o]];L=L[[o]]];
r=Length[L=DeleteCases[L,0]];
i=1;
While[i<r,
If[PolynomialRemainder[L[[i+1]],L[[i]],x]===0,++i,finalstep1[i];i=Max[--i,1]]
];
(a[[#,#]]=L[[#]])&/@Range[r]];
a];


extpolyGCD[p01_,p02_,x_]:=Module[{U=IdentityMatrix[2],p1=p01,p2=p02},While[p2=!=0,U={U[[2]],Expand[U[[1]]-PolynomialQuotient[p1,p2,x]U[[2]]]};{p1,p2}={p2,PolynomialRemainder[p1,p2,x]}];(p1=Expand[p1/#];U[[1]]=Expand[U[[1]]/#])&[Coefficient[p1,x,Exponent[p1,x]]];{p1,U}];


makemonic=Function[i,If[a[[i,i]]=!=0,(a[[i,i]]=Expand[a[[i,i]]/#];u[[i]]=Expand[u[[i]]/#])&[Coefficient[#,x,Exponent[#,x]]&[a[[i,i]]]]]];


makemonic1=Function[i,If[a[[i,i]]=!=0,(a[[i,i]]=Expand[a[[i,i]]/#])&[Coefficient[#,x,Exponent[#,x]]&[a[[i,i]]]]]];


finalstep[i_]:=Module[{g,g0,g1,\[Gamma]0,\[Gamma]1},(g=#[[1]];{g0,g1}=#[[2,1]])&@extpolyGCD[L[[i]],L[[i+1]],x];{\[Gamma]0,\[Gamma]1}=PolynomialQuotient[#,g,x]&/@{L[[i]],L[[i+1]]};
u[[i+{0,1}]]=Expand[{{g0,1},{-g1  \[Gamma]1,\[Gamma]0}}.u[[i+{0,1}]]];v[[All,i+{0,1}]]=Expand[v[[All,i+{0,1}]].{{1,-\[Gamma]1},{g1 ,g0 \[Gamma]0}}];{L[[i]],L[[i+1]]}={g,Expand[\[Gamma]0  L[[i+1]] ]}];


finalstep1[i_]:=Module[{g,g0,g1,\[Gamma]0,\[Gamma]1},(g=#[[1]];{g0,g1}=#[[2,1]])&@extpolyGCD[L[[i]],L[[i+1]],x];{\[Gamma]0,\[Gamma]1}=PolynomialQuotient[#,g,x]&/@{L[[i]],L[[i+1]]};
{L[[i]],L[[i+1]]}={g,Expand[\[Gamma]0  L[[i+1]] ]}];


rowops=Function[i,NestWhile[(If[#[[1]]<Abs[Exponent[a[[i,i]],x]],rowswap[i,#[[2]]]];If[Exponent[a[[i,i]],x]>-\[Infinity],lowerbelow[i];rowmin[i]])&,rowmin[i],(#[[2]]=!=0 ) &]];


rowops1=Function[i,NestWhile[(If[#[[1]]<Abs[Exponent[a[[i,i]],x]],rowswap1[i,#[[2]]]];If[Exponent[a[[i,i]],x]>-\[Infinity],lowerbelow1[i];rowmin[i]])&,rowmin[i],(#[[2]]=!=0 ) &]];


rowmin=Function[i,Fold[If[-\[Infinity]<#2[[1]]<#1[[1]],#2,#1]&,{\[Infinity],0},MapIndexed[{Exponent[#,x],#2[[1]]+i}&,a[[Range[i+1,p],i]]]]];


rowswap=If[#1=!=#2,(a[[{#1,#2}]]=a[[{#2,#1}]];u[[{#1,#2}]]=u[[{#2,#1}]])]&;


rowswap1=If[#1=!=#2,a[[{#1,#2}]]=a[[{#2,#1}]]]&;


lowerbelow=Function[i,With[{q=PolynomialQuotient[a[[#,i]],a[[i,i]],x]},(a[[#]]=Expand[a[[#]]-q a[[i]]];u[[#]]=Expand[u[[#]]-q u[[i]]])]&/@  Range[i+1,p]];


lowerbelow1=Function[i,(a[[#]]=Expand[a[[#]]-PolynomialQuotient[a[[#,i]],a[[i,i]],x] a[[i]]])&/@  Range[i+1,p]];


colops=Function[i,NestWhile[(If[#[[1]]<Abs[Exponent[a[[i,i]],x]],colswap[i,#[[2]]]];If[Exponent[a[[i,i]],x]>-\[Infinity],lowerright[i];colmin[i]])&,colmin[i],(#[[2]]=!=0 ) &]];


colops1=Function[i,NestWhile[(If[#[[1]]<Abs[Exponent[a[[i,i]],x]],colswap1[i,#[[2]]]];If[Exponent[a[[i,i]],x]>-\[Infinity],lowerright1[i];colmin[i]])&,colmin[i],(#[[2]]=!=0 ) &]];


colmin=Function[i,Fold[If[-\[Infinity]<#2[[1]]<#1[[1]],#2,#1]&,{\[Infinity],0},MapIndexed[{Exponent[#,x],#2[[1]]+i}&,a[[i,Range[i+1,m]]]]]];


colswap=(If[#1=!=#2,(a[[All,{#1,#2}]]=a[[All,{#2,#1}]];v[[All,{#1,#2}]]=v[[All,{#2,#1}]])])&;


colswap1=(If[#1=!=#2,(a[[All,{#1,#2}]]=a[[All,{#2,#1}]])])&;


lowerright=Function[i,With[{q=PolynomialQuotient[a[[i,#]],a[[i,i]],x]},(a[[All,#]]=Expand[a[[All,#]]-q a[[All,i]]];v[[All,#]]=Expand[v[[All,#]]-q v[[All,i]]])]&/@  Range[i+1,m]];


lowerright1=Function[i,(a[[All,#]]=Expand[a[[All,#]]-PolynomialQuotient[a[[i,#]],a[[i,i]],x] a[[All,i]]])&/@  Range[i+1,m]];


End[]


Begin["`PrivateI`"]


Options[IntegerSmithForm]={ExtendedForm->True};


IntegerSmithForm::"mat"="The argument to `1` must be a nonempty integer matrix.";


IntegerSmithForm::"ierr"="Unespected internal error in `1`. Please contact pascolet@dimi.uniud.it .";


IntegerSmithForm[at_,opts___?OptionQ]:=(If[\[Not](Min[Dimensions[at]]>0\[And]MatrixQ[at,IntegerQ]),Message[IntegerSmithForm::"mat",IntegerSmithForm];Abort[]];
If[ExtendedForm/.{opts}/.Options[IntegerSmithForm],IntegerSmithForm0[at],IntegerSmithForm1[at]]);


IntegerInvariantFactors[at_]:=(If[\[Not](Min[Dimensions[at]]>0\[And]MatrixQ[at,IntegerQ]),Message[IntegerSmithForm::"mat",IntegerInvariantFactors];Abort[]];
DeleteCases[Tr[#,List],0]&@IntegerSmithForm1[at]);


IntegerSmithForm0[at_]:=Block[{a=at,u,v,p,m,r,i,L},
{p,m}=Dimensions[a];
u=IdentityMatrix[p];v=IdentityMatrix[m];
i=1;
While[(i<=Min[p,m])\[And](\[Not]MatchQ[a[[Range[i,p],Range[i,m]]],{{(0)..}..}|{}]),Which[i<Min[p,m],While[\[Not](MatchQ[a[[Range[i+1,p],i]],{(0)..}|{}]\[And]MatchQ[a[[i,Range[i+1,m]]],{(0)..}|{}]),rowops[i];colops[i]];++i,i===p===m,++i,
i===p,If[i<m,colops[i]];++i,
i===m,If[i<p,rowops[i]];++i,True,Message[IntegerSmithForm::"ierr",IntegerSmithForm0];Abort[]]];
r=--i;
If[r>0,
makepositive/@Range[r];L=Tr[a,List];
With[{o=Ordering[L/.{0->\[Infinity]}],rL=Range[Length[L]]},
u[[rL]]=u[[o]];
v[[All,rL]]=v[[All,o]];
a[[rL,rL]]=a[[o,o]];L=L[[o]]];
r=Length[L=DeleteCases[L,0]];
i=1;
While[i<r,
If[Mod[L[[i+1]],L[[i]]]===0,++i,finalstep[i];i=Max[--i,1]]
];
(a[[#,#]]=L[[#]])&/@Range[r]];
{a,{u,v}}];


IntegerSmithForm1[at_]:=Block[{a=at,p,m,r,i,L},
{p,m}=Dimensions[a];
i=1;
While[(i<=Min[p,m])\[And](\[Not]MatchQ[a[[Range[i,p],Range[i,m]]],{{(0)..}..}|{}]),Which[i<Min[p,m],While[\[Not](MatchQ[a[[Range[i+1,p],i]],{(0)..}|{}]\[And]MatchQ[a[[i,Range[i+1,m]]],{(0)..}|{}]),rowops1[i];colops1[i]];++i,i===p===m,++i,
i===p,If[i<m,colops1[i]];++i,
i===m,If[i<p,rowops1[i]];++i,True,Message[IntegerSmithForm::"ierr",IntegerSmithForm1];Abort[]]];
r=--i;
If[r>0,
makepositive1/@Range[r];L=Tr[a,List];
With[{o=Ordering[L/.{0->\[Infinity]}],rL=Range[Length[L]]},
a[[rL,rL]]=a[[o,o]];L=L[[o]]];
r=Length[L=DeleteCases[L,0]];
i=1;
While[i<r,
If[Mod[L[[i+1]],L[[i]]]===0,++i,finalstep1[i];i=Max[--i,1]]
];
(a[[#,#]]=L[[#]])&/@Range[r]];
a];


extGCD[p01_,p02_]:=Module[{U=IdentityMatrix[2],p1=p01,p2=p02},While[p2=!=0,U={U[[2]],Expand[U[[1]]-Quotient[p1,p2]U[[2]]]};{p1,p2}={p2,Mod[p1,p2]}];(p1=p1 #;U[[1]]=U[[1]] #)&[Sign[p1]];{p1,U}];


makepositive=Function[i,If[a[[i,i]]=!=0,(a[[i,i]]=a[[i,i]] #;u[[i]]=u[[i]] #)&[Sign[a[[i,i]]]]]];


makepositive1=Function[i,If[a[[i,i]]=!=0,(a[[i,i]]=a[[i,i]] #)&[Sign[a[[i,i]]]]]];


finalstep[i_]:=Module[{g,g0,g1,\[Gamma]0,\[Gamma]1},(g=#[[1]];{g0,g1}=#[[2,1]])&@extGCD[L[[i]],L[[i+1]]];{\[Gamma]0,\[Gamma]1}=Quotient[#,g]&/@{L[[i]],L[[i+1]]};
u[[i+{0,1}]]={{g0,1},{-g1  \[Gamma]1,\[Gamma]0}}.u[[i+{0,1}]];v[[All,i+{0,1}]]=v[[All,i+{0,1}]].{{1,-\[Gamma]1},{g1 ,g0 \[Gamma]0}};{L[[i]],L[[i+1]]}={g,\[Gamma]0  L[[i+1]] }];


finalstep1[i_]:=Module[{g,g0,g1,\[Gamma]0,\[Gamma]1},(g=#[[1]];{g0,g1}=#[[2,1]])&@extGCD[L[[i]],L[[i+1]]];{\[Gamma]0,\[Gamma]1}=Quotient[#,g]&/@{L[[i]],L[[i+1]]};
{L[[i]],L[[i+1]]}={g,\[Gamma]0  L[[i+1]] }];


rowops=Function[i,NestWhile[(If[Abs[#[[1]]]<(Abs[a[[i,i]]]/.{0->\[Infinity]}),rowswap[i,#[[2]]]];If[a[[i,i]]=!=0,lowerbelow[i];rowmin[i]])&,rowmin[i],(#[[2]]=!=0 ) &]];


rowops1=Function[i,NestWhile[(If[Abs[#[[1]]]<(Abs[a[[i,i]]]/.{0->\[Infinity]}),rowswap1[i,#[[2]]]];If[a[[i,i]]=!=0,lowerbelow1[i];rowmin[i]])&,rowmin[i],(#[[2]]=!=0 ) &]];


rowmin=Function[i,Fold[If[0<Abs[#2[[1]]]<Abs[#1[[1]]],#2,#1]&,{\[Infinity],0},MapIndexed[{Abs[#],#2[[1]]+i}&,a[[Range[i+1,p],i]]]]];


rowswap=If[#1=!=#2,(a[[{#1,#2}]]=a[[{#2,#1}]];u[[{#1,#2}]]=u[[{#2,#1}]])]&;


rowswap1=If[#1=!=#2,a[[{#1,#2}]]=a[[{#2,#1}]]]&;


lowerbelow=Function[i,With[{q=Quotient[a[[#,i]],a[[i,i]]]},(a[[#]]=a[[#]]-q a[[i]];u[[#]]=u[[#]]-q u[[i]])]&/@  Range[i+1,p]];


lowerbelow1=Function[i,(a[[#]]=a[[#]]-Quotient[a[[#,i]],a[[i,i]]] a[[i]])&/@  Range[i+1,p]];


colops=Function[i,NestWhile[(If[Abs[#[[1]]]<(Abs[a[[i,i]]]/.{0->\[Infinity]}),colswap[i,#[[2]]]];If[a[[i,i]]=!=0,lowerright[i];colmin[i]])&,colmin[i],(#[[2]]=!=0 ) &]];


colops1=Function[i,NestWhile[(If[Abs[#[[1]]]<(Abs[a[[i,i]]]/.{0->\[Infinity]}),colswap1[i,#[[2]]]];If[a[[i,i]]=!=0,lowerright1[i];colmin[i]])&,colmin[i],(#[[2]]=!=0 ) &]];


colmin=Function[i,Fold[If[0<Abs[#2[[1]]]<Abs[#1[[1]]],#2,#1]&,{\[Infinity],0},MapIndexed[{Abs[#],#2[[1]]+i}&,a[[i,Range[i+1,m]]]]]];


colswap=(If[#1=!=#2,(a[[All,{#1,#2}]]=a[[All,{#2,#1}]];v[[All,{#1,#2}]]=v[[All,{#2,#1}]])])&;


colswap1=(If[#1=!=#2,(a[[All,{#1,#2}]]=a[[All,{#2,#1}]])])&;


lowerright=Function[i,With[{q=Quotient[a[[i,#]],a[[i,i]]]},(a[[All,#]]=a[[All,#]]-q a[[All,i]];v[[All,#]]=v[[All,#]]-q v[[All,i]])]&/@  Range[i+1,m]];


lowerright1=Function[i,(a[[All,#]]=a[[All,#]]-Quotient[a[[i,#]],a[[i,i]]] a[[All,i]])&/@  Range[i+1,m]];


End[]


EndPackage[]
