userdet:=proc(a::Matrix,n::integer,userzero,variablelist)::polynomial;
   local sig0,sig0s,u,v,b1,b2:
   local j,k,l,m,c,small,new_small,ones,r:
   local x,y,q,p,p0,ct,cs:

   if n=1 then
     return(a[1,1]):
   end if:

   ct:=Matrix(n,n+1):
   x:=Vector(n):
   y:=Vector(n):
   q:=Vector(n):

   ct[1,1]:=1:
   ct[1,2]:=1:
   for l from 2 to n do
     ct[l,1]:=1:
     for j from 2 to l do 
       ct[l,j]:=ct[l-1,j-1]+ct[l-1,j]:
     end do:
     ct[l,l+1]:=1:
   end do:

   p:=Vector(n):
   for l from 1 to n do
     p[n-l+1]:=a[1,l]:
   end do:

   sig0s:=-1:

   for l from 2 to n do

     for k from 1 to l do
       x[k]:=1:
     end do:

     for k from l+1 to n do
       x[k]:=0:
     end do:

     for k from 1 to n do
       y[k]:=0:
     end do:

     r:=ct[n,l+1]:
     p0:=Vector(r): 

     k:=l:
     m:=1:
     while k<>0 do

       k:=1:
       for j from 1 to n do
         if x[j]=1 then
           q[k]:=j:
           k:=k+1:
         end if:
       end do:

       sig0:=sig0s:
       p0[r-m+1]:=userzero:
 
       for j from 1 to l do
         u:=1:
         v:=l-1:
         b1:=n-1:

         for k from l by -1 to j+1 do

            for b2 from b1 by -1 to q[k] do
               u:=u+ct[b2,v]:
            end do:
            b1:=b2-1:
            v:=v-1:

         end do:

         for k from j-1 by -1 to 1 do

            for b2 from b1 by -1 to q[k] do
               u:=u+ct[b2,v]:
            end do:
            b1:=b2-1:
            v:=v-1:

         end do:

         if sig0=-1 then
	    p0[r-m+1]:=SDMPolynom(p0[r-m+1]-a[l,q[j]]*p[u],variablelist):
         else
	    p0[r-m+1]:=SDMPolynom(p0[r-m+1]+a[l,q[j]]*p[u],variablelist):
         end if:

         sig0:=-sig0:

       end do:

       m:=m+1:

       for k from 1 to n do
         if x[k]<>0 then
           break:
         end if:
       end do:
       small:=k:

       c:=1:
       k:=small:
       while c<>0 and k<=n do
         if x[k]=1 and c=1 then
           x[k]:=0:
         else
           if x[k]=0 and c=1 then
              x[k]:=1:
              c:=0:
           end if:
         end if:
         k:=k+1:
       end do:

       for k from 1 to n do
         if x[k]<>0 then
            break:
         end if:
       end do:
       new_small:=k:

       if new_small<=n then
          ones:=new_small-small-1:

          for k from 1 to ones do
             if x[k]=0 then
                x[k]:=1:
             end if:
          end do:
       end if:

       for k from n by -1 to 1 do
         if x[k]<>y[k] then
           break:
         end if:
       end do:

     end do:

     p:=p0:
     sig0s:=-sig0s:
   end do:

   return(p[1]):

end proc:

with(LinearAlgebra):
with(MTM):

p:=x^12+x^11+a10*x^10+a9*x^9+a8*x^8+a7*x^7+a6*x^6+a5*x^5+a4*x^4+a3*x^3+a2*x^2+a1*x+a0:

q:=diff(p,x):

m12:=BezoutMatrix(p,q,x):

m12:=map(expand,m12):

n12:=max(degree(p,x),degree(q,x)):

vlist:=indets(m12):
nv:=nops(vlist):

exp0:=[seq(vlist[i],i=1..nv)];

m12sdmp:=map(SDMPolynom,m12,exp0):
userzero:=SDMPolynom(0,exp0):

p:=userdet(m12sdmp,n12,userzero,exp0):

q:=convert(p,'polynom'):

nops(q);

quit: