Package $contrib/matrixnormalform

Alias MatNF := $contrib/matrixnormalform;

Define About()
  PrintLn "
    Topic   : Matrix normal forms
    Keywords: matrix, normal form, Smith, rational
    Authors : S. Defrancisci, A. Bigatti
    Version : CoCoA 4.3
    Date    : 1 July 2003
  ";
EndDefine; -- About


Define Man()
  PrintLn "
Suggested alias for this package:
    Alias MatNF := $contrib/matrixnormalform;

SYNTAX
    MatNF.Smith(M: MAT): MAT
    MatNF.SmithFactor(M: MAT): RECORD of MAT

DESCRIPTION 
    MatNF.Smith(M)
  Returns the Smith normal form of M, a matrix over the integers or over
  the univariate polynomials
    MatNF.SmithFactor(M)
  Returns Record[Smith, U, V] where U * M * V = Smith (the Smith normal form)

>EXAMPLE<

M  := Mat([
        [-x + 2, 0, 0, 0],
        [-1, -x + 1, 0, 0],
        [0, -1, -x, -1],
        [1, 1, 1, -x + 2]
        ]);
MatNF.Smith(M);
Mat([
  [1, 0, 0, 0],
  [0, 1, 0, 0],
  [0, 0, x - 1, 0],
  [0, 0, 0, x^3 - 4x^2 + 5x - 2]
])
-------------------------------
SFs := MatNF.SmithFactor(M);
SFs.U * M * SFs.V = SFs.Smith;
TRUE
-------------------------------

M := Mat([
         [2, 0, 0, 0, 0, 0, 0],
         [0, 2, 0, 0, 0, 0, 0],
         [0, 0, 3, 0, 0, 0, 0],
         [0, 0, 0, 3, 0, 0, 0],
         [0, 0, 0, 0, 5, 0, 0],
         [0, 0, 0, 0, 0, 5, 0],
         [0, 0, 0, 0, 0, 0, 5]
         ]);
MatNF.Smith(M);
Mat([
  [1, 0, 0, 0, 0, 0, 0],
  [0, 1, 0, 0, 0, 0, 0],
  [0, 0, 1, 0, 0, 0, 0],
  [0, 0, 0, 1, 0, 0, 0],
  [0, 0, 0, 0, 5, 0, 0],
  [0, 0, 0, 0, 0, 30, 0],
  [0, 0, 0, 0, 0, 0, 30]
])
-------------------------------
SFs := MatNF.SmithFactor(M);
SFs.U * M * SFs.V = SFs.Smith;
TRUE
-------------------------------
";
EndDefine; -- Man


----------------------------------------------------------------------
/*                RIDUZIONE A FORMA NORMALE DI SMITH                */
----------------------------------------------------------------------
/* INDICE 

Funzioni principali:
- Smith (presente in questa sezione)
- Verifica (presente in questa sezione)
- RecDiag (presente in questa sezione)
- LSuperDiag (presente in questa sezione)

Funzioni ausiliarie 1:
- FirstNonZero (presente in questa sezione)
- SeIntOPol (presente in questa sezione)
- PosMinAbs (presente in questa sezione)
- PosMinDegLMinInII (presente in questa sezione)

Funzioni ausiliarie 2:
- MatDiScambioR (presente in questa sezione)
- MatDiScambioC (presente in questa sezione)
- MatDiSommaR (presente in questa sezione)
- MatDiSommaC (presente in questa sezione)
- PDiv (presente in questa sezione)
- PMod (presente in questa sezione)
- Sgn (presente in questa sezione)  
*/


  
Define ManDecom1()
  PrintLn "
-- EXAMPLE 1:

M := [[2,0,0,0],[-1,1,0,0],[0,-1,0,-1],[1,1,1,2]];
Decom1(M);

/* Interpretazione dei risultati: V si decompone in somma diretta di
due sottomoduli ciclici w1 e w2 con w1 = e2-e3, w2 = -e1+e3. */

-- EXAMPLE 2:

M := [[0,0,-1],[1,1,1],[1,0,2]];
Decom1(M);

/* Interpretazione dei risultati: V si decompone in somma diretta di
due sottomoduli ciclici w1 e w2 con w1 = e2, w2 = -e2+e3. */
";
EndDefine; -- ManDecom1

  
Define ManDecom2()
  PrintLn "
-- EXAMPLE 1:

M := [[2,0,0,0],[-1,1,0,0],[0,-1,0,-1],[1,1,1,2]];
Decom2(M);

/* Interpretazione dei risultati: V si decompone in somma diretta di
tre sottomoduli ciclici w1 e z1 e z2 con w1 = e2-e3, z1 = e2-2e3, z2 =
-e1+e2-e4.  Confrontandolo con l'Esempio 1 precedente si osserva che
il modulo ciclico <w1> rimane immutato, mentre il modulo ciclico <w2>
si e' spezzato nella somma diretta dei moduli ciclici irriducibili
<z1> e <z2>. */ 

-- EXAMPLE 2:

M := [[0,0,-1],[1,1,1],[1,0,2]];
Decom2(M);

/* Interpretazione dei risultati: V si decompone in somma diretta di
due sottomoduli ciclici w1 e w2 con w1 = e2-e3, w2 = -e1+e3.
Confrontandolo con l'Esempio 2 precedente si osserva che la
decomposizione rimane immutata, perche' i moduli ciclici erano gia'
irriducibili.*/
";
EndDefine; -- ManDecom2


Define ManFormaCanonicaRazionale()
  PrintLn "
-- EXAMPLE 1:

M := [[0,0,-1],[1,1,1],[1,0,2]];
FormaCanonicaRazionale(M);
/*
 Mat[
  [1, 0, 0],
  [0, 0, -1],
  [0, 1, 2]
]
 */

-- EXAMPLE 2:

M := [[2,0,0,0],[-1,1,0,0],[0,-1,0,-1],[1,1,1,2]];
FormaCanonicaRazionale(M);
/*
 Mat[
  [1, 0, 0, 0],
  [0, 0, 0, 2],
  [0, 1, 0, -5],
  [0, 0, 1, 4]
]
 */

-- EXAMPLE 3:

M := [[0,0,0,4],[1,0,0,0],[0,1,0,0],[0,0,1,0]];
FormaCanonicaRazionale(M);
/*
 Mat[
  [0, 0, 0, 4],
  [1, 0, 0, 0],
  [0, 1, 0, 0],
  [0, 0, 1, 0]
]
 */"
EndDefine; -- ManFormaCanonicaRazionale


Define SmithFactor(A)
  --trova la Forma Normale di Smith di M
  NR := Len(A);  NC := Len(A[1]);
  IdNR := Identity(NR);  IdNC := Identity(NC);
  A1 := Mat(A);
  N := $.SeIntOPol(A1,NR,NC);  A1 := N[1];  Pol := N[2];
  If NR <= NC
    Then L := $.RecDiag([IdNR, A1, IdNC], 1, Pol);
  Else
    LT := $.RecDiag([IdNC,Transposed(A1),IdNR],1, Pol); 
    L := [Transposed(LT[3]),Transposed(LT[2]),Transposed(LT[1])];
  EndIf;
  Return Record(U=L[1], V=L[3], Smith=L[2]);
EndDefine; -- SmithFactor


Define Smith(A)
  SF := $.SmithFactor(A);
  Return SF.Smith;
EndDefine; -- Smith


Define Verify(M,L)
  U := L.U;
  MM := Mat(M);
  V := L.V;
  D := L.Smith;
  For I := 1 To Len(D) Do
    For J := 1 To Len(D[1]) Do
      If I<>J And D[I,J]<>0 Then
        Error("non diagonal: ", D);
      EndIf;
    EndFor;
  EndFor;
  For I := 2 To Min(Len(D), Len(D[1])) Do
    If D[I-1,I-1]<>0 Then 
      If $.PMod(D[I,I], D[I-1,I-1])<>0 Then Error('non SmithNF', D);
      EndIf;
    EndIf;
  EndFor;
  If U*MM*V=D Then Print 'OK : U*M*V = SmithNF'
  Else Print 'ERROR : U*M*V <> SmithNF';
  EndIf;
EndDefine; -- Verifica

-----------------------------------------------------------------------
/* La funzione RecDiag e' costruita in base ai passi 2 e 3 della
dimostrazione del Teorema 3.1.3.

Usando solo operazioni razionali trasforma la matrice M in una
diagonale equivalente con entrate non ancora ridotte .  */

Define RecDiag(L,I,Pol)
  -- assumiamo che la matrice abbia NR <= NC
  U := L[1];  M := L[2];  V := L[3];
  NR := Len(M);  NC := Len(M[1]);
  If I>NR Then Return $.LSuperDiag([U,M,V],1,Pol);  EndIf;
  Mlist := $.LMinInII(M,I,Pol);
  U := Mlist[1]*U;   M := Mlist[2];   V := V*Mlist[3];
  If M[I,I]=0 Then Return $.LSuperDiag([U,M,V],1,Pol); EndIf;
  -- lavora sulla i-ma riga
  For J := I+1 To NC Do
    While M[I,J]<>0 Do
      Q := $.PDiv(M[I,J],M[I,I]);
      ES := $.MatDiSommaC([I,J],-Q,NC);
      M := M*ES;  V := V*ES;
      -- cerca il nuovo minimo
      Mlist := $.LMinInII(M,I,Pol); 
      U := Mlist[1]*U;   M := Mlist[2];   V := V*Mlist[3];
    EndWhile;
  EndFor;
  -- lavora sulla i-ma colonna
  IR := I+1;
  OK := TRUE;
  While IR<=NR And OK Do
    If M[IR,I]<>0 Then
      Q1 := $.PDiv(M[IR,I],M[I,I]);
      ES1 := $.MatDiSommaR([IR,I],-Q1,NR);
      M := ES1*M;  U := ES1*U;
    EndIf;
    If M[IR,I]<>0 Then OK := FALSE; EndIf;
    IR := IR+1;
  EndWhile;
  -- chiamata ricorsiva
  If OK Then Return $.RecDiag([U,M,V], I+1, Pol);
  Else Return $.RecDiag([U,M,V], I, Pol);
  EndIf;
EndDefine; -- RecDiag

----------------------------------------------------------------------
/* La funzione LSuperDiag raffina la diagonalizzazione ottenuta da
RecDiag in base all'Osservazione 3.1.5. */

Define LSuperDiag(LD, I, Pol) 
  -- raffinamento in base all'osservazione 3.1.5
  U := LD[1];  M := LD[2];  V := LD[3];
  NR := Len(M);  NC := Len(M[1]);
  If I+1>NR Then Return [U,M,V] EndIf;
  S := $.MatDiSommaR([I,I+1],1,NR);
  M := S*M;  U := S*U;
  Mlist := $.LMinInII(M,I,Pol);
  U := Mlist[1]*U;   M := Mlist[2];   V := V*Mlist[3];
  If M[I,I]=0 Then Return [U,M,V] EndIf;
  --rende monico M[I,I] o cambia segno a colonna I
  If Pol Then X := 1/LC(M[I,I]); Else X := $.Sgn(M[I,I]); EndIf;
  If X <> 1 Then 
    E := $.MatDiSommaC([I,I],X,NC);
    M := M*E;  V := V*E;
  EndIf;
  -- lavora sulla i-ma riga
  While M[I,I+1]<>0 Do
    --rende monico M[I,I+1] o cambia segno a colonna I+1
    If Pol Then X := 1/LC(M[I,I+1]); Else X := $.Sgn(M[I,I+1]); EndIf;
    If X <> 1 Then 
      E := $.MatDiSommaC([I,I+1],X,NC);
      M := M*E;  V := V*E;
    EndIf; 
    Q := $.PDiv(M[I,I+1],M[I,I]);  
    ES := $.MatDiSommaC([I,I+1],-Q,NC);
    M := M*ES;  V := V*ES; 
    Mlist := $.LMinInII(M,I,Pol); --cerca il nuovo minimo
    U := Mlist[1]*U;   M := Mlist[2];   V := V*Mlist[3];
    --rende monico M[I,I] o cambia segno a colonna I
    If Pol Then X := 1/LC(M[I,I]); Else X := $.Sgn(M[I,I]); EndIf;
    If X <> 1 Then 
      E := $.MatDiSommaC([I,I],X,NC);
      M := M*E;  V := V*E;
    EndIf;  
  EndWhile; 
  -- lavora sull'elemento diagonale utilizzando l'i-ma colonna
  IR := I+1;
  OK := TRUE;
  While IR<=NR And OK Do
    If M[IR,I]<>0  Then 
      --rende monico M[IR,I] oppure cambia segno a M[IR,I]
      If Pol Then X := 1/LC(M[IR,I]); Else X := $.Sgn(M[IR,I]); EndIf;
      If X <> 1 Then 
	E := $.MatDiSommaR([IR,IR],X,NR);
        M := E*M;  U := E*U;
      EndIf;
      Q1 := $.PDiv(M[IR,I],M[I,I]); 
      ES1 := $.MatDiSommaR([IR,I],-Q1,NR);
      M := ES1*M;   U := ES1*U;
    EndIf;  
    If M[IR,I]<>0 Then OK := FALSE; EndIf;
    If Pol Then X := 1/LC(M[IR,IR]); Else X := $.Sgn(M[IR,IR]); EndIf;
    If X <> 1 Then 
      E := $.MatDiSommaC([IR,IR],X,NC);
      M := M*E;  V := V*E;
    EndIf;
    IR := IR+1;
  EndWhile;
  -- chiamata ricorsiva
  If OK Then 
    S := $.LSuperDiag([U,M,V],I+1,Pol);
    If $.PMod(S[2][I+1,I+1],S[2][I,I])=0 Then Return S;
    Else  Return $.LSuperDiag(S,I,Pol);  EndIf; 
  Else Return $.LSuperDiag([U,M,V], I,Pol);
  EndIf;
EndDefine; --LSuperDiag

----------------------------------------------------------------------

-- Funzioni ausiliarie 1

/* Le seguenti funzioni ausiliarie individuano l'elemento minimo ( in
valore assoluto o in grado) e lo spostano in (I,I). In particolare
FirstNonZero prende in input la matrice M, le dimensioni NR = numero
righe di M, NC = numero colonne di M e l'intero I che indica la
sottomatrice che sto considerando ad un generico passo e restituisce
l'elemento piu' piccolo diverso da zero.
 
SeIntOPol legge la matrice e controlla se vi e' almeno un'entrata
polinomiale.  Se la trova trasforma anche gli eventuali interi in
polinomi e Pol = True.  Se non la trova dichiara che la matrice ha
entrate intere e Pol = False.  Se Pol = False entra in gioco PosMinAbs
che prende gli stessi input della funzione FirstNonZero e restituisce
la posizione dell'elemento minimo in valore assoluto.
 
Se Pol = True entra in gioco PosMinDeg che prende gli stessi input
della funzione FirstNonZero e restituisce la posizione dell'elemento
minimo in grado.  Infine LMinInII prende in input la matrice M,
l'intero I e il booleano Pol e restituisce la lista [E1,Mat(M),E2]
dove E1 e' la matrice delle operazioni fatte sulle righe per spostare
l'elemento minimo nel posto (I,I), M e' la matrice di prima in cui
pero' l'elemento minimo si trova in (I,I) e E2 e' la matrice delle
operazioni fatte sulle colonne. */

Define FirstNonZero(M,NR,NC,I)
  --trova il primo elemento della matrice diverso da zero 
  For IP := I To NR Do
    For JP := I To NC Do
      If M[IP, JP]<>0 Then Return M[IP, JP]; EndIf;
    EndFor;
  EndFor;
  Return 0;
EndDefine; -- FirstNonZero


Define PosMinAbs(M,NR,NC,I)
  --posizione del minimo per matrici di interi
  MinM := Abs($.FirstNonZero(M, NR, NC, I))+1;
  If MinM=1 Then Return [I,I]; EndIf; 
  MinI := -1;  MinJ := -1;
  For IP := I To NR Do
    For JP := I To NC Do
      If M[IP,JP]<>0 And Abs(M[IP,JP])<MinM Then
        MinM := Abs(M[IP,JP]);
        MinI := IP;  MinJ := JP;
      EndIf;
    EndFor;
  EndFor;
  Return [MinI,MinJ];
EndDefine;  -- PosMinAbs


Define PosMinDeg(M,NR,NC,I)
  --posizione del minimo per matrici di polinomi
  MinM := Deg($.FirstNonZero(M, NR, NC, I))+1;
  If MinM=1 Then Return [I,I]; EndIf; 
  MinI := -1;  MinJ := -1;
  For IP := I To NR Do
    For JP := I To NC Do
      If M[IP,JP]<>0 And  Deg(M[IP,JP])<MinM Then 
          MinM := Deg(M[IP,JP]);
          MinI := IP;   MinJ := JP;          
          EndIf;
    EndFor;
  EndFor;
  Return [MinI,MinJ];
EndDefine;  --PosMinDeg


Define LMinInII(M,I,Pol)
  --sposta il minimo in (I,I) 
  NR := Len(M);  NC := Len(M[1]);
  If Pol Then L := $.PosMinDeg(M,NR,NC,I);
  Else L := $.PosMinAbs(M,NR,NC,I);
  EndIf;
  IL := L[1];  JL := L[2];
  If I<=NR And I<=NC Then
    If IL<>I Then 
      E1 := $.MatDiScambioR([IL,I],NR);
      M := E1*M;
    Else E1 := Identity(NR);       
    EndIf;
    If JL<>I Then 
      E2 := $.MatDiScambioC([I,JL],NC);
      M := M*E2;
    Else E2 := Identity(NC);              
    EndIf;
    Return [E1,Mat(M),E2];
  Else -- se le  matrici sono rettangolari
    Return 'Indice I troppo grande';
  EndIf;
EndDefine;  -- LMinInII

----------------------------------------------------------------------

-- Funzioni ausiliarie 2 

/* La funzione MatDiScambioR prende in input la lista L = [i,j] degli
indici e la lunghezza NR delle righe della matrice per costruire la
matrice elementare Eij che scambia la riga i con riga j. Poiche'
l'algoritmo funziona per matrici rettangolari la funzione
MatDiScambioR e' costruita per essere moltiplicata a sinistra e quindi
per scambiare le righe, mentre MatDiScambioC e' costruita per essere
moltiplicata a destra e quindi per scambiare le colonne.
 
La funzione MatDiSommaR prende in input la lista L = [i,j] degli
indici, l'intero A e la lunghezza NR delle righe della matrice per
costruire la matrice elementare E(i,j,A) che aggiunge alla riga i la
riga j moltiplicata per A. Come sopra la funzione MatDiSommaR e'
costruita per essere moltiplicata a sinistra e quindi per operare
sulle righe, mentre MatDiSommaC e' costruita per essere moltiplicata a
destra e quindi per operare sulle colonne. */

Define MatDiScambioR(L,NR)
  I := L[1];
  J := L[2];
  E := Identity(NR);
  E[I,I] := 0;
  E[J,J] := 0;
  E[I,J] := 1;
  E[J,I] := 1;
  Return E;
EndDefine;  -- MatDiScambioR

Define MatDiScambioC(L,NC)
  I := L[1];
  J := L[2];
  E := Identity(NC);
  E[I,I] := 0;
  E[J,J] := 0;
  E[I,J] := 1;
  E[J,I] := 1;
  Return E;
EndDefine;  -- MatDiScambioC

Define MatDiSommaR(L,A,NR)
  E := Identity(NR);
  I := L[1];
  J := L[2];
  E[I,J] := A;
  Return E ;
EndDefine;  -- MatDiSommaR

Define MatDiSommaC(L,A,NC)
  E := Identity(NC);
  I := L[1];
  J := L[2];
  E[I,J] := A;
  Return E ;
EndDefine;  -- MatDiSommaC


/* La funzione PDiv prende in input due elementi che possono essere
sia interi sia polinomi e restituisce il risultato della divisione tra
interi nel primo caso e il risultato della divisione tra polinomi nel
secondo.
             
La funzione PMod prende in input due elementi che possono essere sia
interi sia polinomi e restituisce il resto della divisione tra interi
nel primo caso e il resto della divisione tra polinomi nel secondo.

La funzione Sgn prende in input un intero e restituisce 1 se e'
positivo e -1 se e' negativo. */

Define PDiv(F,G)
  If Type(F)=INT And Type(G)=INT Then Return Div(F,G);
  Else
   E := DivAlg(F,[G]);
   Return Head(E.Quotients);
  EndIf;
EndDefine;  -- PDiv 

Define PMod(F,G)
  If Type(F)=INT And Type(G)=INT Then Return Mod(F,G);
  Else
   E := DivAlg(F,[G]);
   Return E.Remainder;
  EndIf;
EndDefine;  -- PMod

Define Sgn(X)
  If X>0 Then Return 1; EndIf;
  If X<0 Then Return -1; EndIf;
  Return 0;
EndDefine; -- Sgn



----------------------------------------------------------------------------
/*        DECOMPOSIZIONE DI V COME SOMMA DIRETTA DI MODULI CICLICI        */
----------------------------------------------------------------------------
/* INDICE 

Funzione principale:
- Decom1 (presente in questa sezione)

Funzioni ausiliarie 1:
- SmithOfxI_M (presente in questa sezione)
- Rho (presente in questa sezione)

Funzioni ausiliarie 2:
- Smith (in sezione "Riduzione a Forma Normale di Smith")
*/

-- Funzione principale

/* La funzione Decom1 prende in input la matrice di interi M e restituisce la l
ista dei vettori wi da leggere come coefficienti della base canonica. */

Define Decom1(M)
  ML := $.SmithOFxI_M(M);
  U  := ML.U;
  MD := ML.Smith;
  NZ := Len([ I In 1..Len(M) | Deg(MD[I,I]) <> 0 ]);
  Return Rho(NZ,Inverse(U),M);
EndDefine; -- Decom1

--Funzioni ausiliarie 1

/* In particolare la Funzione principale chiama per prima SmithOFxI_M(M) che pr
ende in input la matrice M e restituisce la Forma Normale di Smith della matric
e (XI-M), richiamando su (XI-M) la funzione Smith descritta nella sezione "Ridu
zione a Forma Normale di Smith". La seconda funzione richiamata e' Rho che pren
de in input NZ = non zero,ossia il numero di elementi della diagonale di Smith 
con grado maggiore di zero, UI = inversa della matrice U,dove U e' la matrice d
i operazioni fatte sulle righe per trasformare (XI-M) in Forma di Smith e M la 
matrice di interi di partenza.
La funzione Rho in base all'Osservazione 4.2.6 calcola prima la base (v1,...,vn
)=(e1,...,en)U^(-1) e poi usa le formule w1=rho(v_(n-r+1)),...,wr=rho(vn). */

Define SmithOFxI_M(M)
  Id := Identity(Len(M));
  A := x*Id-Mat(M);
  Return $.SmithFactor(A);
EndDefine; -- SmithOFxI_M

Define Rho(NZ,UI,M) 
  TM := Transposed(Mat(M));
  E := Transposed(Mat(UI));
  NR := Len(M);
  Id := Identity(NR);
  W := NewMat(NZ,NR,0);
  C := NewList(NR,0);
  For K := 1 To NZ Do
    For J := 1 To NR Do
      C[J] := Coefficients(Poly(E[NR-NZ+K,J]));
      If  C[J] <> [0] Then 
            D := Len(C[J]); 
            For I := 1 To D Do
              If I-1 = 0 Then TT := Id; Else TT := TM^(I-1);
              EndIf;  
              W[K] := W[K]+(C[J][D-I+1])*TT[J];
            EndFor;
      EndIf;
    EndFor;
  EndFor;
  Return List(W);
EndDefine; -- Rho


---------------------------------------------------------------------------
/* DECOMPOSIZIONE DI V COME SOMMA DIRETTA DI MODULI CICLICI IRRIDUCIBILI */
---------------------------------------------------------------------------
/* INDICE 

Funzione principale:
- Decom2 (presente in questa sezione)

Funzioni ausiliarie 1:
- Fatt (presente in questa sezione)
- SmithOfxI_M (in sezione "Decomposizione di V come somma diretta di moduli cic
lici")
- Rho (in sezione "Decomposizione di V come somma diretta di moduli ciclici")

Funzioni ausiliarie 2:
- Smith (in sezione "Riduzione a Forma Normale di Smith")
*/

-- Funzione principale

/* La funzione Decom2 prende in input la matrice di interi M e restituisce in b
ase al Teorema 3.1.31 l'ulteriore fattorizzazione di V. Ricordiamo, infatti, ch
e se 0:m =  (d) con d = p1^r1....ps^rs dove p1,...,ps sono elementi primi disti
nti e r1,...,rs sono interi positivi,se si pone di = d|(pi^ri) allora si ha che
 <m> e' uguale alla somma diretta dei <d1 m>....<ds m>. In particolare Decom2 c
ostruisce la decomposizione di V con lo stesso procedimento di Decom1, ma usand
o al posto del fattore invariante un fattore alla volta della sua scomposizione
. Osserviamo, inoltre, che usando l'algoritmo di fattorizzazione Fatt si puo' p
rocedere a questa decomposizione solo ammenttendo l'esistenza di un tale algori
tmo. Le altre funzioni ausiliarie usate sono SmithOFxI_M e Rho descritte nella 
sezione "Decomposizione di V come somma diretta di moduli ciclici",  e Fatt des
critta qui di seguito. */


Define Decom2(M)
  ML := $.SmithOFxI_M(M);
  U  := ML.U;
  MD := ML.Smith;
  NR := Len(M);
  LL := Fatt(NR,MD);
  F := LL[1];
  NZ := LL[2];
  Id := Identity(NR);
  W := Rho(NZ,Inverse(U),M);
  V := NewList(NZ,0);
  TM := Transposed(Mat(M));
  For I := 1 To NZ Do
    S := Len(F[I]); 
    If S=1 Then
      V[I] := W[I];
    Else
      V1 := NewMat(S,NR,0);
      P := NewList(S,0);
      C := NewList(S,0);
      For J := 1 To S Do
        P[J] := (F[I,J][1])^(F[I,J][2]); 
        C[J] := Coefficients(P[J]);
        DJ := Len(C[J]); 
        For H := 1 To DJ Do
          If H=1 Then 
            If Deg(P[J])>=1 And DJ=1 Then TT := TM;
            Else TT := Id;
            EndIf;
          Else TT := TM^(H-1);
          EndIf;
          For T := 1 To NR Do
            V1[J] := V1[J]+(W[I,T]*C[J][DJ-H+1])*TT[T];
          EndFor;
        EndFor;
      EndFor;
      Append(V,List(V1)); 
    EndIf;  
  EndFor; 
  VV := [];
  For I := 1 To Len(V) Do
    If V[I]<>0 Then Append(VV,V[I]); EndIf;
  EndFor;
  Return VV;
EndDefine; -- Decom2


-- Funzioni ausiliarie 1

/* La funzione Fatt prende in input NR = numero di righe, MD = Smith(XI-M) e re
stituisce la lista di F = fattorizzazione dei fattori invarianti,NZ = numero de
gli elementi di grado maggiore di zero,cioe' non costanti,presenti sulla diagon
ale della matrice MD. Per fattorizzare i fattori invarianti usa Factor, funzion
e gia' esistente nella versione utilizzata di "CoCoA language." */ 


Define Fatt(NR,MD)
  NZ := Len([ I In 1..Len(MD) | Deg(MD[I,I]) <> 0 ]); 
  F := [ MD[I,I] | I In 1..Len(MD) And MD[I,I] <> 1 ];
  For I := 1 To NZ Do
    F[I] := Factor(F[I]);   
  EndFor; 
  Return [F,NZ];
EndDefine; -- Fatt






--------------------------------------------------------------------------
/*           FORMA CANONICA RAZIONALE DI UNA MATRICE QUADRATA           */
--------------------------------------------------------------------------
/* INDICE 

Funzione principale:
- FormaCanonicaRazionale (presente in questa sezione)

Funzioni ausiliarie 1:
- FattoriInvarianti (presente in questa sezione)
- SmithOfxI_M (in sezione "Decomposizione di V come somma diretta di moduli cic
lici")
- BloccoFattInv (presente in questa sezione)

Funzioni ausiliarie 2:
- Smith (in sezione "Riduzione a Forma Normale di Smith")
*/

-- Funzione principale

/* La funzione FormaCanonicaRazionale calcola la Forma Canonica
Razionale di una matrice quadrata nel modo seguente. Trova i fattori
invarianti di XI-M tramite la funzione FattoriInvarianti (descritta
nelle funzioni ausiliarie). Per ogni fattore invariante costruisce il
blocco razionale elementare corrispondente chiamando la funzione
BlocoFattInv (descritta nelle funzioni ausiliarie) e infine
restituisce la matrice cercata. Osserviamo che usa solo operazioni
razionali.
*/

Define FormaCanonicaRazionale(M)
  LF := FattoriInvarianti(M);
  N := Len(LF);
  MB := NewMat(N,N,0);
  For I := 1 To N Do
    MB[I,I] := BloccoFattInv(LF[I]);
  EndFor;
  Return BlockMatrix(MB);
EndDefine; -- FormaCanonicaRazionale
  
  
-- Funzioni ausiliarie 1: 
 
/* La funzione FattoriInvarianti prende in input una matrice quadrata
M e restituisce la lista dei fattori invarianti di (XI-M). La funzione
SmithOFxI_M e' la stessa descritta nella sezione "Decomposizione di V
come somma diretta di moduli ciclici", mentre la funzione
BloccoFattInv prende in input un polinomio di grado n (nel nostro caso
un fattore invariante alla volta) e restituisce la corrispondente
matrice razionale elementare (o companion matrix) che ha 1 sulla sotto
diagonale e i coefficienti del polinomio nell'ultima colonna a partire
dal termine noto cambiato di segno fino al coefficiente di grado n-1
sempre cambiato di segno. */
  
  
Define FattoriInvarianti(M)
  --trova i fattori invarianti di (XI-M)
  ML := $.SmithOFxI_M(M);
  MD := ML.Smith;
  Return [ MD[I,I] | I In 1..Len(MD) And MD[I,I]<>1 ];
EndDefine; -- FattoriInvarianti

  
Define BloccoFattInv(P)
  --matrice razionale elementare di P
  If P=0 Then Return 0; EndIf;
  CoefPol := Tail(Coefficients(P,x));
  N := Len(CoefPol);
  B := NewMat(N,N,0);
  For I := 1 To N Do
    If I >= 2 Then B[I,I-1] := 1;EndIf;
    B[I,N] := -CoefPol[N-I+1];
  EndFor;
  Return B;
EndDefine; -- BloccoFattInv




---------------------------------------------------------------------------
/*             FORMA QUASI DI JORDAN DI UNA MATRICE QUADRATA             */
---------------------------------------------------------------------------
/* INDICE 

Funzione principale:
- QuasiJordan (presente in questa sezione)

Funzioni ausiliarie 1:
- TrovaQM (presente in questa sezione)
- TrovaCM (presente in questa sezione)
- FattoriInvarianti (in sezione "Forma Canonica Razionale di una matrice quadra
ta")

Funzioni ausiliarie 2:
- TrovaQ (presente in questa sezione)
- TrovaC (presente in questa sezione)
- BloccoFattInv (in sezione "Forma Canonica Razionale di una matrice quadrata")

Funzioni ausiliarie 3:
- TrovaQAux (presente in questa sezione)
- Pillar (presente in questa sezione)
- PDiv (in sezione "Riduzione a Forma Normale di Smith") 

*/

-- Funzione principale

/*La funzione QuasiJordan costruisce, in base al Teorema 6.1.7, le matrici QM e
 CM tali che QM^-1*CM*QM = QuasiJordan. La matrice QM viene costruita dalla fun
zione ausiliaria TrovaQM tramite il calcolo dei fattori invarianti di M e la co
struzione della Q corrispondente. La matrice CM viene costruita sempre tramite 
il calcolo dei fattori invarianti di M concludendo, pero', con la costruzione d
ella corrispondente matrice C.
Una descrizione piu' dettagliata si trova di seguito in Funzioni ausiliarie.  
Osserviamo che se i fattori invarianti si spezzano in polinomi lineari irriduci
bili, allora QuasiJordan restituisce la Forma Canonica di Jordan; se i fattori 
invarianti si spezzano in polinomi non lineari irriducibili, allora QuasiJordan
 restituisce la Forma Quasi di Jordan della matrice M. */


Define QuasiJordan(M)
  Q := $.TrovaQM(M); 
  Q1 := Inverse(Q);
  C := $.TrovaCM(M); 
  QJ := Q1*C*Q;
  Return QJ;
EndDefine; -- QuasiJordan

  
  
-- Funzioni ausiliarie 1: 

/*La funzione TrovaQM prende in input una matrice quadrata e restituisce la mat
rice QM che contiene sulla diagonale principale i blocchi Q relativi ad ogni fa
ttore invariante di M. I blocchi Q sono costruiti dalla funzione TrovaQ, descri
tta in ausiliarie 2. La funzione TrovaCM prende in input una matrice quadrata e
 restituisce la matrice CM che contiene sulla diagonale principale i blocchi C 
relativi ad ogni fattore invariante di M. I blocchi C sono costruiti dalla funz
ione TrovaC, descritta in ausiliarie 2. Entrambe utilizzano la funzione Fattori
Invarianti gia' descritta in sezione "Forma Canonica Razionale di una matrice q
uadrata." */


Define TrovaQM(M)
  LFI := $.FattoriInvarianti(M);
  N := Len(LFI);
  QM := NewMat(N,N,0);
  For I := 1 To N Do
    QM[I,I] := $.TrovaQ(LFI[I]);
  EndFor;
  Return BlockMatrix(QM);
EndDefine; -- TrovaQM
  
Define TrovaCM(M)
  LFC := $.FattoriInvarianti(M); 
  N := Len(LFC);
  CM := NewMat(N,N,0);
  For I := 1 To N Do
    CM[I,I] := $.TrovaC(LFC[I]);
  EndFor;
  Return BlockMatrix(CM);
EndDefine; -- TrovaCM
 
 
-- Funzioni ausiliarie 2: 
 
/*

 La funzione TrovaQ prende in input il polinomio P, che rappresenta
ogni fattore invariante di M, lo fattorizza con la funzione Factor e
per ogni fattore di P trova il blocco corrispondente chiamando
TrovaQAux. Infine restituisce la ma trice Q = [Q1,...,Qn].

 La funzione TrovaC prende in input il polinomio P, lo fa ttorizza con
la funzione Factor e per ogni fattore di P trova il blocco corrisp
ondente chiamando BloccoFattInv. Supponiamo che P := F1^n1...FS^ns,
allora Trov aC contiene sulla diagonale principale le matrici
razionali elementari di ogni Fi^ni.

 Infine ricordiamo che la funzione BloccoFattInv e' gia' stata
descritta nella sezione "Forma Canonica Razionale di una matrice
quadrata."

 */


Define TrovaQ(P)
  LP := Factor(P);
  R := Len(LP);
  AA := NewMat(R,R,0);
  For K := 1 To R Do
    PP := LP[K,1];
    N := LP[K,2];
    AA[K,K] := $.TrovaQAux(PP,N);
  EndFor;
  Return BlockMatrix(AA);
EndDefine; -- TrovaQ

Define TrovaC(P)
  LP := Factor(P); 
  R := Len(LP); 
  AA := NewMat(R,R,0);
  For K := 1 To R Do
    PP := LP[K,1]; 
    N := LP[K,2];
    AA[K,K] := $.BloccoFattInv(PP^N); 
  EndFor;
  Return BlockMatrix(AA);
EndDefine; -- TrovaC


-- Funzioni ausiliarie 3:

/*

 La funzione TrovaQAux prende in input il fattore del polinomio e il
rispettivo grado e ne calcola la corrispondente
(0,p^(n-j+1),m,p^n)-pillar, chiamando la funzione Pillar. La funzione
Pillar prende in input la costante T che nel nostro caso e' sempre
uguale a 0, P che e' il fattore che stiamo considerando del polinomio
F, M che e' il grado del fattore P e F che e' il polinomio non ancora
 fattorizzato di grado n.

 La funzione Pillar calcola la matrice di tipo n x M l a cui k-esima
colonna contiene i coefficienti, eventualmente completati con 0, del
polinomo (x)^(k-1)(F/P), come voluto dalla definizione 6.1.1. In
conclusion e ricordiamo che la funzione PDiv e' gia' stata descritta
nella sezione "Riduzi one a Forma Normale di Smith."

 */


Define TrovaQAux(PP,N)
--trova Qi dalle pillar
  P := PP^N;
  M := Deg(PP);
  Q1 :=  NewMat(N*M,M,0);
  For J := 1 To M Do
    Q1[J,J] := 1
  EndFor; 
  QJ := [ $.Pillar(0,PP^(N-(J-1)),M,P) | J In 2..N ];
  A := Concat([Q1],QJ);
  Return BlockMatrix([A]);
EndDefine; -- TrovaQAux

Define Pillar(T,P,M,F)
--trova le (t,p,m,f)-pillar
  N := Deg(F);
  Q := NewMat(N,M,0);
  G := $.PDiv(F,P);
  For I := 1 To N Do
    Q[I,1] := CoeffOfTerm(x^(I-1),G);
  EndFor;
  For K := 2 To M Do
    G := (x-T)*G;
    For J := 1 To N Do
      Q[J,K] := CoeffOfTerm(x^(J-1),G);
    EndFor;
  EndFor;
  Return Q;
EndDefine; -- Pillar

----------------------------------------------------------------------
-- Check Input
----------------------------------------------------------------------

  Define SeIntOPol(M,NR,NC)
  --trova se entrate di interi o polinomi
  P := 0;
  For I := 1 To NR Do
    For J := 1 To NC Do
      If M[I,J] <> 0 Then 
        If Deg(M[I,J]) <> 0 Then P := P+1 EndIf;
      EndIf; 
    EndFor;
  EndFor;
  If P = 0 Then Pol := False;
  Else For I := 1 To NR Do
    For J := 1 To NC Do
      M[I,J] := Poly(M[I,J]);
    EndFor;
  EndFor;
    Pol := True;
  EndIf;
  Return [M,Pol];
EndDefine; -- SeIntOPol



  
EndPackage;


/*


----------------------------------------------------------------------
-- EXAMPLE 1:


/* Osserviamo che in questo caso troviamo
la Forma Quasi di Jordan di M. */
M := Mat([
	 [0,0,0,0,0,-125],
	 [1,0,0,0,0,-150],
	 [0,1,0,0,0,-135],
	 [0,0,1,0,0,-68],
	 [0,0,0,1,0,-27],
	 [0,0,0,0,1,-6]]);
$contrib/matrixnormalform.QuasiJordan(M);


/*
 Mat[
  [0, -5, 0, 0, 0, 0],
  [1, -2, 0, 0, 0, 0],
  [0, 1, 0, -5, 0, 0],
  [0, 0, 1, -2, 0, 0],
  [0, 0, 0, 1, 0, -5],
  [0, 0, 0, 0, 1, -2]
]
-------------------------------
 */


----------------------------------------------------------------------
-- EXAMPLE 2:
/* Osserviamo che in questo caso troviamo
la Forma Quasi di Jordan di M. */

M := [[0,0,0,4],[1,0,0,0],[0,1,0,0],[0,0,1,0]];
$contrib/matrixnormalform.QuasiJordan(M);
/*
 Mat[
  [0, 2, 0, 0],
  [1, 0, 0, 0],
  [0, 0, 0, -2],
  [0, 0, 1, 0]
]
-------------------------------
 */


----------------------------------------------------------------------
-- EXAMPLE 3:
/* Osserviamo che in questo caso troviamo proprio
la Forma Canonica di Jordan di M. */

M := [[2,0,0,0],[-1,1,0,0],[0,-1,0,-1],[1,1,1,2]];
$contrib/matrixnormalform.QuasiJordan(M);
/*
 Mat[
  [1, 0, 0, 0],
  [0, 2, 0, 0],
  [0, 0, 1, 0],
  [0, 0, 1, 1]
]
-------------------------------
 */


 */