userddet:=proc(a::Matrix,n::integer)::integer;

   local r,s,t,u,v,S0,S1;
   local x,y,z;

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

   x:=Vector(2):
   z:=Vector(n):

   x[2]:=a[1,1]:

   for r from 1 to n-1 do

     y:=Vector(r+2):

     z[1]:=a[r+1,r+1]:

     S0:=Vector(r):
     for t from 1 to r do
       S0[t]:=a[t,r+1]:
     end do:
 
     if r<n-1 then
       y[2]:=max(z[1],x[2]):
     end if:

     for s from 2 to r+1 do

       for t from s-1 by -1 to 1 do
         z[t+1]:=z[t]: 
       end do:

       z[1]:=0:
       for t from 1 to r do
         z[1]:=max(z[1],a[r+1,t]+S0[t]):       
       end do:

       if s<>r+1 then
         S1:=Vector(r):
         for u from 1 to r do
           S1[u]:=0:
           for v from 1 to r do
             S1[u]:=max(S1[u],a[u,v]+S0[v]):
           end do:
         end do:
         S0:=S1:
       end if:

       if r<n-1 then
         if s<r+1 then
           y[s+1]:=z[1]:
           for t from 2 to s do 
             y[s+1]:=max(y[s+1],z[t]+x[t]):
           end do:
	   y[s+1]:=max(y[s+1],x[s+1]):
         else
           y[s+1]:=z[1]:
           for t from 2 to s do 
             y[s+1]:=max(y[s+1],z[t]+x[t]):
           end do:
         end if:
       else
         if s=r+1 then
           for t from 2 to s do 
             z[1]:=max(z[1],z[t]+x[t]):
           end do:
         end if:
       end if:

     end do:

     x:=y:
   end do:

   return(z[1]):

end proc:

det_maxdegree:=proc(a::Matrix,n::integer)::integer;
   local vlist0,n0,d0,r0,row0,colm0,b1,b2:
   local j,k,l:

   b1:=Matrix(n,n):
   b2:=Matrix(n,n):
   vlist0:=indets(a):
   n0:=nops(vlist0):
   r0:=0:

   for l from 1 to n0 do
     for j from 1 to n do
        for k from 1 to n do
           d0:=degree(a[j,k],vlist0[l]):
           if d0>0 then
             b1[j,k]:=d0:
             b2[k,j]:=d0:
           else
             b1[j,k]:=0:
             b2[k,j]:=0:
           end if:
        end do:
     end do:

     row0:=userddet(b1,n):
     colm0:=userddet(b2,n):

     if row0<colm0 then
        d0:=row0:
     else
        d0:=colm0:
     end if:

     if r0<d0 then
        r0:=d0:
     end if:

   end do:

   return(r0):
end proc:

userdet:=proc(a::Matrix,n::integer,userzero,usermone)::polynomial;

   local r,s,t,u,v,S0,S1;
   local x,y,z,w;

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

   x:=Vector(2):
   z:=Vector(n):

   x[2]:=sdmp:-Multiply(usermone,a[1,1]):

   for r from 1 to n-1 do

     y:=Vector(r+2):

     z[1]:=sdmp:-Multiply(usermone,a[r+1,r+1]):

     S0:=Vector(r):
     for t from 1 to r do
       S0[t]:=a[t,r+1]:
     end do:

     if r<n-1 then
       y[2]:=sdmp:-Add(z[1],x[2]):
     end if:

     for s from 2 to r+1 do

       for t from s-1 by -1 to 1 do
         z[t+1]:=z[t]: 
       end do:

       w:=userzero:
       for t from 1 to r do
         w:=w,sdmp:-Multiply(a[r+1,t],S0[t]):       
       end do:
       w:=sdmp:-Add(w):
       z[1]:=sdmp:-Multiply(usermone,w):

       if s<>r+1 then
         S1:=Vector(r):
         for u from 1 to r do
           w:=userzero:
           for v from 1 to r do
             w:=w,sdmp:-Multiply(a[u,v],S0[v]):
           end do:
           S1[u]:=sdmp:-Add(w):
         end do:
         S0:=S1:
       end if:

       if r<n-1 then
         if s<r+1 then
           w:=z[1]:
           for t from 2 to s do 
             w:=w,sdmp:-Multiply(z[t],x[t]):
           end do:
           w:=w,x[s+1];
           y[s+1]:=sdmp:-Add(w):
         else
           w:=z[1]:
           for t from 2 to s do 
             w:=w,sdmp:-Multiply(z[t],x[t]):
           end do:
           y[s+1]:=sdmp:-Add(w):
         end if:
       else
         if s=r+1 then
           for t from 2 to s do 
             z[1]:=z[1],sdmp:-Multiply(z[t],x[t]):
           end do:
           z[1]:=sdmp:-Add(z[1]):
         end if:
       end if:

     end do:

     x:=y:
   end do:

   if (n mod 2)=1 then
     z[1]:=sdmp:-Multiply(usermone,z[1]):
     return(z[1]):
   else
     return(z[1]):
   end if:

end proc:

with(LinearAlgebra):
with(MTM):

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:

q:=diff(p,x):

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

m9:=map(expand,m9):

n9:=9:

d0:=det_maxdegree(m9,n9);
d1:=ceil(log2(d0+1.0));
d2:=floor(64/d1);

if d2=0 then
  print("error"):
  quit:
end if:

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

if d2>nv then
  d2:=nv:
end if:

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

m9sdmp:=map(sdmp:-Import,m9,plex(exp0),pack=d2):

userzero:=sdmp:-Import(0,plex(exp0),pack=d2):
usermone:=sdmp:-Import(-1,plex(exp0),pack=d2):

p:=userdet(m9sdmp,n9,userzero,usermone);

quit: