with(linalg): 

idealsource:=proc(f,part)

local U,s,p,y1,nf,y2,i,j,k,da,dv,t,H,L,T,n,V,A,h,fz,varf,vf,zz,g,su;

global W;

print(´f´=f);

nf:=nops(f); varf:=indets(f); vf:=nops(varf); fz:=f; g:=f[nf];

if nf <> vf then

   ERROR("Dimensions of source and target are different");

fi:

for i to (nf-1) do

  fz[nf]:=subs(f[i]=0,fz[nf]):

   if degree(fz[nf]) <=0 then

     ERROR("f is not finite,verify input!");

   fi;

od;

su:=0;

  for i to nops(part) do

    su:=part[i]+su;

  od:

if su > vf then

   ERROR(`Use partition for integer <= source dimension.`);

fi:

 zz:=indets(fz[nf])[];

 print(Partition = part);

 y1:=0;

 y2:=0;

  for j to nops(part) do

      for k from 0 to part[j] do

         y1:=y1+1;

         y2:=y2+1:

         H[y1]:=z[j,k];

         t[y2]:=z[j,k+1]=z[j,0];            

           if k = part[j] then  

              y2:=y2-1;

            fi;

      od; 

    od;

  L:=[seq(H[n], n=1..y1)]:

  T:=[seq(t[n], n=1..y2)];

  y2:=0:

  y1:=0:

  s:=0;

    for n from 2 to nops(L) do

      V:=vandermonde(L):

        dv:=factor(det(V)): 

          for i to nops(L) do  

             V[i,n]:=subs(zz=V[i,2],g):

          od:

              A:=convert(V,matrix):

             da:=factor(det(A)):

              s:=s+1;

        U[n-1]:=subs(T,simplify(da/dv)):

        print(h[s]=U[n-1]);

        W[s]:=U[n-1]:

    od:

end: