(* first example *)

BmT = {{(135897/208970),(4316466/104485),(3266372/104485),(216107/208970),-(5001151/41794),(10166539/208970)},{-(702887/417940),-(3097393/104485),-(4332107/208970),-(593317/417940),(6870083/83588),-(15212989/417940)},{-(67353/20897),-(1225779/20897),-(929352/20897),-(49406/20897),(3487602/20897),-(1456703/20897)},{(503887/417940),(3266998/104485),(4362757/208970),(104697/417940),-(7253779/83588),(16157549/417940)},{-(33817/20897),-(709947/20897),-(496288/20897),-(37268/20897),(1950421/20897),-(849805/20897)},{-(602357/417940),-(2784933/104485),-(4022217/208970),-(555847/417940),(6408125/83588),-(14677399/417940)}};

(* second example *)

BmT2={{-(1113739/208970),(3397616/104485),(509822/104485),(951187/208970),-(8613471/208970),(2154167/208970)},{(1072219/417940),-(4889561/208970),-(208216/104485),-(1637497/417940),(11065201/417940),-(3831427/417940)},{(110059/20897),-(964879/20897),-(146652/20897),-(153766/20897),(1160374/20897),-(319179/20897)},{-(1382399/417940),(5147021/208970),(100916/104485),(1214277/417940),-(11525261/417940),(4063127/417940)},{(66789/20897),-(561997/20897),-(52438/20897),-(96448/20897),(630707/20897),-(204743/20897)},{(1077889/417940),-(4334391/208970),-(157896/104485),-(1544227/417940),(9999751/417940),-(3904057/417940)}};

(* let us compute the MAP *)

(* first let us introduce the canonical vectors *)

e1 = Table[0,{i,1,6}]; e1[[1]]=1;
e2 = Table[0,{i,1,6}]; e2[[2]]=1;
e3 = Table[0,{i,1,6}];e3[[3]]=1;
e4 = Table[0,{i,1,6}]; e4[[4]]=1;
e5 = Table[0,{i,1,6}]; e5[[5]]=1;
e6 = Table[0,{i,1,6}];e6[[6]]=1;

NullSpace[NestList[BmT.#&,e1,4]] //Length

(NestList[BmT.#&,e1,4] //NullSpace )87

genE1=NestList[BmT.#&, e1,4]

(* EvalPol replaces the variable la with the matrix A *)

EvalPol= Module[{A,pol,powA},
A = #1;
pol = #2;
powA=NestList[A.#&, IdentityMatrix[Length[A]],Length[CoefficientList[pol,la]]-1];
(CoefficientList[pol,la]).powA
]&;

carPol=Det[la IdentityMatrix[6]-BmT]

EvalPol[BmT,carPol]

Eigenvalues[BmT]

polTest1=PolynomialQuotient[carPol,la+2,la]

EvalPol[BmT,polTest1]

polTest2=PolynomialQuotient[polTest1,la+2,la]

(* by extracting another eigenvalue although the polynomial degree decreases, it fails to cancel the matrix *)
EvalPol[BmT,polTest2]

polTest3 = PolynomialQuotient[polTest1,la+3,la]

(* ... same with the extraction of the other eigenvalue *)EvalPol[BmT,polTest3]

polTest1

NullSpace[NestList[BmT.#&, e1,10]]

CoefficientList[polTest1,la].NestList[BmT.#&, e1,5]

(* Now we complement the cyclic space witha complementary space and there are many possible ways to do this, for example *)perp2=List[NullSpace[genE1][[1]]+genE1[[3]]+4genE1[[4]]]

PP=Join[genE1,perp2]

(* another possible choice is to stick with the null space of the cyclic susbpace as complementary subspace *)

PP0=Join[genE1,NullSpace[genE1]]

(* then let us change the coordinates to achieve an equivalent matrix (from the invariant polynomial perspective) where the canonical structure becomes more visible *)BmTPP=Inverse[Transpose[PP]].BmT.Transpose[PP]

BmTPP0=Inverse[Transpose[PP0]].BmT.Transpose[PP0]

BmTPP //MatrixForm

BmTPP0 //MatrixForm

BmTPP //Last//Last

(* Now let us perform the same analysis with BmT2 *)

BmT2

carPol2 = Det[la IdentityMatrix[6]-BmT2]

EvalPol[BmT2,carPol2]

Eigenvalues[BmT2]

poly21=PolynomialQuotient[carPol2, (la+2),la]

Eigenvalues[BmT2]

poly21=PolynomialQuotient[carPol2, (la+2),la]

EvalPol[BmT2,poly21]

poly22=PolynomialQuotient[poly21, (la+2),la]

EvalPol[BmT2,poly21]

poly23=PolynomialQuotient[poly22, (la+2),la]

EvalPol[BmT2,poly23]

poly24=PolynomialQuotient[poly23, (la+3),la]

poly25=PolynomialQuotient[poly23, (la+2),la]

EvalPol[BmT2,poly24]

EvalPol[BmT2,poly25]

(* the map is therefore poly23 *)

poly23

NestList[BmT2.#&, e1,10] //NullSpace

(* this is too large *)

NestList[BmT2.#&, e2,10] //NullSpace

NestList[BmT2.#&, e3,10] //NullSpace

NestList[BmT2.#&, e4,10] //NullSpace

NestList[BmT2.#&, e5,10] //NullSpace

NestList[BmT2.#&, e6,10] //NullSpace

gen2=NestList[BmT2.#&,e3,2]

PP2=Join[gen2,NullSpace[gen2]]

BmT2PP=Inverse[Transpose[PP2]].BmT2.Transpose[PP2]; BmT2PP //MatrixForm

BmQ=Drop[#,3]& /@ Drop[BmT2PP,3]

Solve[poly23 == 0,la]

Eigenvalues[BmQ]

Solve[carPol2==0,la]

(* let us use the classical library of functions to achieve the decomposition and computation of the invariant polynomials *)

Directory[]

(* ATTENTION: the following line should be changed to where you have put the file "SmithFormV6.m" *)

SetDirectory["~/Dropbox/Cours/EDOC-LinSyst/LinSystSoftware"]

<<SmithFormV6`

A1 = {{la-2,0,0},{0,la-2,0},{0,0,la-2}}

A2 = -{{la,1,0},{-4,-4+la,0},{0,0,la-2}}

A3 = la IdentityMatrix[3]-{{0,1,0},{0,0,1},{-8,-12,-6}}

{S1,{U1,V1}}=PolynomialSmithForm[A1,la]

{S2,{U2,V2}}=PolynomialSmithForm[A2,la];

{S3,{U3,V3}}=PolynomialSmithForm[A3,la];

{U1//MatrixForm, V1// MatrixForm}

S1//MatrixForm

S2//MatrixForm

S3//MatrixForm

cl=-CoefficientList[(s+3)^3 //Expand,s]//Most

c2=-CoefficientList[(s+2)^3 //Expand,s]//Most

Eigenvalues[{{0,1,0},{0,0,1},cl}]

Eigenvalues[{{0,1,0},{0,0,1},c2}]

(* let us compute the structure with the Smith form of the quotient matrix BmQ *)

{S1T,{U1T,V1T}} = PolynomialSmithForm[la IdentityMatrix[6]-BmT,la];

{S2T,{U2T,V2T}} = PolynomialSmithForm[la IdentityMatrix[6]-BmT2,la];

S1T//MatrixForm

U1T//MatrixForm

U1T.(la IdentityMatrix[6]-BmT).V1T//Simplify

S2T//MatrixForm

{Sq,{Uq,Vq}} = PolynomialSmithForm[la IdentityMatrix[3]-BmQ,la];

Sq//MatrixForm