with(PolynomialIdeals):with(linalg): with(LinearAlgebra): with(Groebner):
#This finction computes the minimal poly of matrices over Q by MinPoly1 and MinimalPolynomial.
mp := proc (M, x)
local MM,n, N1, N2, i, j, G, f, g, h, t1, t2, b1, b2;
b1, t1 := kernelopts(bytesused, cputime);
n := op(M)[1];
N1 := M;
N2 := NULL;
MM:=M:
for i from 2 to n do
    MM:=Multiply(M,MM):
    N1 := N1, MM:
end do;
for i to n do
    for j to n do
        if i <> j then
            N2 := N2, sum(a[k]*[N1][k][i, j], k = 1 .. n):
        else
            N2 := N2, sum(a[k]*[N1][k][i, j], k = 1 .. n)+c:
        end if:
    end do:
end do;
G := Basis([N2], plex(seq(a[n+1-k], k = 1 .. n), c));
f := NormalForm(sum(a[k]*x^k, k = 1 .. nops(G))+c, subs(seq(a[i] = 0, i = nops(G)+1 .. n), G), plex(seq(a[n+1-k], k = 1 .. nops(G)), c, x));
g := collect(simplify(f/LeadingCoefficient(f,plex(x))), x);
b2, t2 := kernelopts(bytesused, cputime);
#print(g);
print(t2-t1, b2-b1);
print("$$$$$$$$$$$$$$$$$");
b1, t1 := kernelopts(bytesused, cputime);
h := collect(MinimalPolynomial(M, x), x);
b2, t2 := kernelopts(bytesused, cputime);
print(t2-t1, b2-b1);
print("$$$$$$$$$$$$$$$$$");
print(evalb(g = h)):
end proc:

with(linalg): with(LinearAlgebra): with(Groebner): with(PolynomialIdeals):
##This procedure computes the inverse of a monomial in algebraic number field.
inv:=proc(m,J,T)
local h,f,i,u,plc,sum,g,S,L,G:
#global ns,nn;
#option trace;
G:=Basis([op(J),z*m-1],plex(z,op(T))):
RETURN(z-G[-1]/LeadingCoefficient(G[-1],plex(z,op(T))));
end:
#This finction computes the minimal poly of matrices over Alg field num by MinPoly2.
mp4 := proc (M,J,V,x)
local HC,GG,F,IHC, N1, N2, i, j, G, f, g, h, t1, t2, b1, b2;
#global n,ns,nn;
#option trace:
b1, t1 := kernelopts(bytesused, cputime);
ns:=NormalSet(J,plex(op(V)))[1];
nn:=nops(ns);
n := op(M)[1];
N1 := M;
N2 := NULL;
MM:=M:
for i from 2 to n do
    MM:=MM.M;
    N1 := N1, MM:
end do;
for i to n do
    for j to n do
        if i <> j then
            N2 := N2, sum(a[k]*[N1][k][i, j], k = 1 .. n):
        else
            N2 := N2, sum(a[k]*[N1][k][i, j], k = 1 .. n)+c:
        end if:
    end do:
end do;
G := Basis([N2,op(J)], plex(seq(a[n+1-k], k = 1 .. n), c,op(V)));
GG:=subs(seq(a[i] = 0, i = nops(G)-nops(J)+1 .. n), G):
F:=[seq(NormalForm(inv(LeadingCoefficient(GG[i],plex(seq(a[n+1-k], k = 1 .. n),c)),J,V)*GG[i],J,plex(op(V))) , i=2..nops(GG))]:F:=[op({op(F)} minus {0})]:
f := NormalForm(sum(a[k]*x^k, k = 1 .. nops(G)-nops(J))+c, F, plex(seq(a[n+1-k], k = 1 .. nops(G)), c,op(V)));
f:=simplify(f/LeadingCoefficient(f,plex(x,op(V)))):
HC:=LeadingCoefficient(f,plex(x));
IHC:=inv(HC,J,V);
f:=NormalForm(IHC*f,J,plex(op(V)));
g := collect(simplify(f/LeadingCoefficient(f,plex(x))),x);
b2, t2 := kernelopts(bytesused, cputime);
#print(g);
print(t2-t1, b2-b1);
print("$$$$$$$$$$$$$$$$$");
end proc:

# Following function computes, the minimal polynomial of a matrix over alg. extension of finite field.
pinv:=proc(m,J,T,p)
local h,f,i,u,plc,sum,g,S,L,G:
#global ns,nn;
#option trace;
G:=Basis([op(J),z*m-1],plex(z,op(T)),characteristic=p):
RETURN(z-G[-1]/LeadingCoefficient(G[-1],plex(z,op(T))));
end:
pmp4 := proc (M,J,V,x,p)
local HC,GG,F,IHC, N1, N2, i, j, G, f, g, h, t1, t2, b1, b2;
#global n,ns,nn;
b1, t1 := kernelopts(bytesused, cputime);
ns:=NormalSet(J,plex(op(V)),characteristic=p)[1];
nn:=nops(ns);
n := op(M)[1];
N1 := M;
N2 := NULL;
MM:=M:
for i from 2 to n do
    MM:=MM.M;
    N1 := N1, MM:
end do;
for i to n do
    for j to n do
        if i <> j then
            N2 := N2, sum(a[k]*[N1][k][i, j], k = 1 .. n):
        else
            N2 := N2, sum(a[k]*[N1][k][i, j], k = 1 .. n)+c:
        end if:
    end do:
end do;
G := Basis([N2,op(J)], plex(seq(a[n+1-k], k = 1 .. n), c,op(V)),characteristic=p);
GG:=subs(seq(a[i] = 0, i = nops(G)-nops(J)+1 .. n), G):
F:=[seq(NormalForm(pinv(LeadingCoefficient(GG[i],plex(seq(a[n+1-k], k = 1 .. n),c)),J,V,p)*GG[i],J,plex(op(V)),characteristic=p) , i=2..nops(GG))]:F:=[op({op(F)} minus {0})]:
f := NormalForm(sum(a[k]*x^k, k = 1 .. nops(G)-nops(J))+c, F, plex(seq(a[n+1-k], k = 1 .. nops(G)), c,op(V)),characteristic=p);
f:=simplify(f/LeadingCoefficient(f,plex(x,op(V)))):
HC:=LeadingCoefficient(f,plex(x));
IHC:=pinv(HC,J,V,p);
f:=NormalForm(IHC*f,J,plex(op(V)),characteristic=p);
g := collect(simplify(f/LeadingCoefficient(f,plex(x))),x);
b2, t2 := kernelopts(bytesused, cputime);
print(g);
print(t2-t1, b2-b1);
end proc:


for i in {5,10,15,20,25,30} do
    p:='p':q:='q':
    print(i);
    M:=RandomMatrix(i,i):
    M[1,1]:=p;M[2,2]:=q:
    mp4(M,[9*p^2-2,8*q^3-3],[p,q],x);
    p:=sqrt(2)/3:q:=3^(1/3)/2:
    b1, t1 := kernelopts(bytesused, cputime);
    h := collect(MinimalPolynomial(M, x), x):
    b2, t2 := kernelopts(bytesused, cputime);
    print(t2-t1, b2-b1);
    print("**********");
    print("$$$000$$$ooo$$$000$$$ooo$$$000$$$ooo$$$000$$$");
od: