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