include"mubase.mmx";

 firstcompanionmatrix(M: Matrix Generic)== {
   d:Int:=degree M;
   m:Int:=rows M;
   n:Int:=columns M;
   L:Matrix Rational:=coefficients(M,0);
   A:Matrix Rational:=transpose(L);
   for i in 1..d do {
     B:Matrix Rational:=transpose(coefficients(M,i));
     A:=horizontal_join(A,B);    
   }
    k:Int:=(d-1)*m;
   C:Matrix Rational:=identity_matrix(k); 
   D:Matrix Rational:=zero_matrix_rational(k,m);
   U:Matrix Rational:=horizontal_join(D,C); 
   if rows U=0 then return(A)
   else  Z:Matrix Rational:=vertical_join(U,A);
   return(Z);
   
 };
 secondcompanionmatrix(M: Matrix Generic)== {
   d:Int:=degree M;
   m:Int:=rows M;
   n:Int:=columns M;
   k:Int:=(d-1)*m;
   A:Matrix Rational:=(-1)*transpose(coefficients(M,d));
   if k=0 then return(A) else {
   B:Matrix Rational:= zero_matrix_rational(n,k);
   V:Matrix Rational:=horizontal_join(B,A);
   C:Matrix Rational:=identity_matrix(k);
   D:Matrix Rational:=zero_matrix_rational(k,m);
   U:Matrix Rational:=horizontal_join(C,D);
   Z:Matrix Rational:=vertical_join(U,V);}
   return(Z);
 };

remove_row(M: Matrix Rational,k: Int)=={
   
   if rows M = 1 then { Z:Matrix Rational:=[rational(0)| i in 0..0 || j in 0..0];return(Z);};
   else  Z:Matrix Rational:=[rational(0)|i in 0..columns M || j in 0..(rows M -1)];
   
    for j in 0..columns M do
    {  for i in 0..k do
         Z[i,j]:= M[i,j];
        for i in (k+1)..rows M do
          Z[i-1,j]:= M[i,j];
    }
 return(Z);
 };
 
remove_column(M: Matrix Rational,k: Int)==
{ if columns M=1 then  { Z:Matrix Rational:=[rational(0)| i in 0..0 || j in 0..0];return(Z);}
  else  Z:Matrix Rational:=[rational(0)| i in 0..(columns M-1)|| j in 0..rows M];
    for i in 0..rows M do
        {for j in 0..k do
                Z[i,j]:= M[i,j];
           for j in (k+1)..columns M do
              Z[i,j-1]:= M[i,j];}
  return(Z);
};
 

regular(M: Matrix Rational, N: Matrix Rational):Vector Generic=={
  A:Matrix Rational:=M; B: Matrix Rational:=N;
  m:Int:=rank B;
  if m = columns B then {A:Matrix Rational:= transpose A; B:Matrix Rational:=transpose B;}
   V:Vector Generic:=column_reduced_echelon_with_transform B;
   A1:Matrix Rational:=A*V[1]; 
   A2:Matrix Rational:=A1;
   for i in 0..m do 
     A1:=remove_column(A1,0);
      U:Vector Generic:=row_reduced_echelon_with_transform A1;
   n:Int:=rank A1;
   A2:Matrix Rational:=U[1]*A2;
   B2:Matrix Rational:=U[1]*V[0];
   for i in m..columns B do
   {A2:=remove_column(A2,m);
     B2:=remove_column(B2,m);}
   for i in 0..n do
   {A2:=remove_row(A2,0);
     B2:=remove_row(B2,0);}
   return[A2,B2];
     
 };

pencilregular(M: Matrix Rational,N: Matrix Rational): Vector Generic=={
  A:Matrix Rational:=M; B:Matrix Rational:=N;
  while rank B< max(rows B, columns B) do
    {C:Vector Generic:=regular(A,B);
  A:=C[0];
  B:=C[1];}
  return[A,B];
};


as_matrix_double (M : Matrix Rational) : Matrix Double == {[ M[i,j] :> 
  Double | i in 0..columns(M) || j in 0..rows(M) ]};

as_matrix_double (M: Matrix Double): Matrix Double =={return(M);};

generalized_eigenvalues(M: Matrix Generic)=={
  A:Matrix Rational:=firstcompanionmatrix M;
  B:Matrix Rational:=secondcompanionmatrix M;
  C: Vector Generic:=pencilregular(A,B);
  Z:Matrix Rational:=[rational(0)| i in 0..0|| j in 0..0];
  if C[1]= Z then return([]);
  else {
  D:Matrix Rational:=C[0]*invert C[1];
  listeigen:Vector Complex Double := eigenvalues( as_matrix_double D);}
  return(listeigen);
};