def det_berkowitz(A,N){

   if (N==1){
     return(A[0][0])$
   }

   X=newvect(1)$
   Z=newvect(N)$

   X[0]=-A[0][0]$

   for(Cr=1;Cr<=N-1;Cr++){

     Z[0]=-A[Cr][Cr]$

     Y=newvect(Cr+1)$
     if (Cr<N-1){
       Y[0]=Z[0]+X[0]$
     }

     for(Cs=1;Cs<=Cr;Cs++){

       if(Cs==1){
         Z[Cs]=0$
         for(Ct=0;Ct<Cr;Ct++){
           Z[Cs]=Z[Cs]-A[Cr][Ct]*A[Ct][Cr]$
         }

         if (Cs!=Cr){
           S0=newvect(Cr)$
           for(Cu=0;Cu<Cr;Cu++){
             S0[Cu]=0$
             for(Cv=0;Cv<Cr;Cv++){
               S0[Cu]=S0[Cu]+A[Cu][Cv]*A[Cv][Cr]$
             }
           }
         }
       }
       else{
         Z[Cs]=0$
         for(Ct=0;Ct<Cr;Ct++){
           Z[Cs]=Z[Cs]-A[Cr][Ct]*S0[Ct]$
         }

         if (Cs!=Cr){
           S1=newvect(Cr)$
           for(Cu=0;Cu<Cr;Cu++){
             S1[Cu]=0$
             for(Cv=0;Cv<Cr;Cv++){
               S1[Cu]=S1[Cu]+A[Cu][Cv]*S0[Cv]$
             }
           }
           S0=S1$
         }
       }

       if(Cr<N-1){
         if(Cs<Cr){
           Y[Cs]=Z[Cs]$
           for(Ct=0;Ct<Cs;Ct++){
              Y[Cs]=Y[Cs]+Z[Cs-Ct-1]*X[Ct]$
           }
           Y[Cs]=Y[Cs]+X[Cs]$
         }
         else{
           Y[Cs]=Z[Cs]$
           for(Ct=0;Ct<Cs;Ct++){
              Y[Cs]=Y[Cs]+Z[Cs-Ct-1]*X[Ct]$
           }
         }
       }else{
         if(Cs==Cr){
           for(Ct=0;Ct<Cs;Ct++){
             Z[N-1]=Z[N-1]+Z[Cs-Ct-1]*X[Ct]$
           }
           if (irem(N,2)==1)
             return(-Z[N-1])$
           else
             return(Z[N-1])$
          }
       }

     }

     X=Y$
   }

}

def seki_res(C,A,B){

	M=deg(A,C)$
	N=deg(B,C)$

	if(M>=N){
		L=M$
		S=N$
	}
	else{
		L=N$
		S=M$
		T=A$
		A=B$
		B=T$
	}

	MAT=newmat(L,L)$

	for(I=0;I<(L-S);I++)
		for(J=0;J<=S;J++)
			MAT[I][J+I]=coef(B,J,C)$

	B=C^(L-S)*B$

	F=srem(A,C^L)$
	G=srem(B,C^L)$
	H=coef(A,L,C)*G-coef(B,L,C)*F$
	for(J=0;J<L;J++){
		MAT[2*L-S-L][J]=coef(H,J,C)$
	}

	for(I=L-1;I>L-S;I--){
		H=C*srem(H,C^(L-1))+coef(A,I,C)*G-coef(B,I,C)*F$
		for(J=0;J<L;J++){
			MAT[2*L-S-I][J]=coef(H,J,C)$
		}
	}

	return MAT$
}

P=a0+a1*x+a2*x^2+a3*x^3+a4*x^4+a5*x^5+a6*x^6+a7*x^7+x^8+x^9$
MAT=seki_res(x,P,diff(P,x))$
Q=det_berkowitz(MAT,9)$

nmono(Q);

end$