//SINGULAR Example 2.6.15

//================== we need here =======================

proc extendedNormalForm(matrix M)
{
   int n=nrows(M);
   int m=ncols(M);
   intvec v=1..n; 
   intvec w=n+1..n+m; 
   intvec u=1..m; 
   intvec x=m+1..n+m; 
   matrix E=unitmat(n);
   matrix B=unitmat(m);
   matrix N=M;            //to keep M for the test
   matrix D,K;             
   
   while(D!=N)
   {
      D=N;
      K=transpose(interred(transpose(concat(E,D))));  
      E=submat(K,v,v);
      N=submat(K,v,w);
      K=interred(transpose(concat(transpose(B),transpose(N))));
      K=simplify(K,1);
      B=submat(K,u,u);
      N=submat(K,x,u);
   }
   matrix C=inverse(E);
   if(M*B!=C*D){ERROR("something went wrong");}  //test
   list Re=B,C,D;
   return(Re);  
}

//=======================================================



option(redSB);
ring R=0,(x),(C,dp);
LIB"matrix.lib";
LIB"linalg.lib";

matrix M[5][5]=1, 1,0,0,0,
              -2,-1,0,0,0,
               0, 0,2,1,0,
               0, 0,0,2,0,
               0, 0,0,0,2;

matrix N[5][5]=1,2, 2, 2,-1,
               1,1, 2, 1, 1,
              -1,1, 2,-1, 2,
              -1,1, 1,-1, 2,
               1,2,-1, 2, 1;
M=lift(N,freemodule(nrows(N)))*M*N;
print(M);

matrix A=M-x*freemodule(5);

list L= extendedNormalForm(A);

print(L[2]);    //the new basis

print(L[3]);   //the diagonal form


matrix V1[5][4]=concat(L[1][1],M*L[1][1],M*M*L[1][1],
                                      M*M*M*L[1][1]);
matrix V2[5][1]=L[1][2];  //the 2 invariant subspaces

list F=factorize(L[2][1,1]);
F;

proc polyOfEndo(matrix B,poly p)
{
  int i;
  int d=nrows(B);
  matrix A=coeffs(p,var(1));
  matrix E[d][d]=freemodule(d);
  matrix C[d][d]=A[1,1]*E;
  for(i=2;i<=nrows(A);i++)
  {
     E=E*B;
     C=C+A[i,1]*E;
  }
  return(C);
}

matrix S=polyOfEndo(M,F[1][3]^2);//the decomposition of V1
matrix V11=std(S*V1);
print(V11);

S=polyOfEndo(M,F[1][2]); 
matrix V12=std(S*V1);
print(V12);

matrix B=concat(V11,V12,V2);
det(B);                 //we obtained a basis
reduce(M*V11,std(V11)); //the subspaces are invariant
reduce(M*V12,std(V12));
reduce(M*V2,std(V2));

matrix C=lift(B,M*B); //the matrix with respect to the
print(C);             //basis B

matrix v[5][1]=V12[1];//special basis for Jordan form
B=concat(V11,M*v-2*v,v,V2);
C=lift(B,M*B);        //the matrix with respect to the
print(C);             //basis B