//SINGULAR Example 2.6.11

//================== we need here =======================
proc normalForm(matrix M)
{
  int n=nrows(M);
  int m=ncols(M);
  matrix N,K;
  matrix L[n][m];
  while(N!=M)
  {
     N=M;
     M=L;
     K=transpose(interred(transpose(interred(N))));     
     M[1..nrows(K),1..ncols(K)]=K;
  }
  return(N);
}

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

ring R=0,(x),(C,dp);

matrix M0[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,3;
matrix N[5][5]=1, 1, -1,  1,-1,
               2, 2,  1,  1, 0,
              -1, 2,  2,  1, 1,
              -2, 1,  1, -1, 0,
               1, 2, -2,  1, 1;

matrix M=lift(N,freemodule(nrows(N)))*M0*N-x*freemodule(5);
print(M);

N=normalForm(M);
print(N);

factorize(N[1,1]);

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);  
}


LIB"matrix.lib";
LIB"linalg.lib";

matrix M1[2][2]=x2+1,  0,
                0   ,x-1;
matrix N1[2][2]=1, 1,
                1, 2;
matrix N2[2][2]=0,-1,
                1, 1;
                
M=N1*M1*N2;
print(M);

list L=extendedNormalForm(M);
print(L[1]);

print(L[2]);

print(L[3]);