# Macaulay and rectangular Macaulay multiresultants for MapleV.4
# Use & distribute freely <Bill.Triggs@inrialpes.fr>  20-May-99 

#---------------------------------------------------------------
# List all monomials: 
# - of a (collection of) polynomial(s): monomials(poly,[vars])
# - of a given degree: monomials(vars,d)
# - between two degrees: monomials(vars,d1..d2)

monomials := (u,v) ->
   if nargs<2 then 
      procname(u,indets(u))
   elif v::integer then
      procname(expand((convert(u,`+`))^v),u)
   elif v::`..` then
      procname(expand((1+convert(u,`+`))^(op(2,v)-op(1,v))
      	       	      *(convert(u,`+`))^op(1,v)),u)
   elif u::{list,set,table} then
      op(sort(convert(map(procname,convert(u,set),v),`+`),v,'tdeg'))
   else # u,v = poly,vars
      op(sort(collect(u,v,'distributed',1),v,'tdeg'))
   fi:

#---------------------------------------------------------------------
# Given a list of polynomials and a list of monomials, build a matrix
# with one row per polynomial and columns labelled by the monomials,
# whose entries are the coefficients of the polynomials at the
# monomials. Used for building multiresultant matrices.

resultant_matrix :=
proc(Polys::list, Monos::list, Vars)
   local mtab,A,i,p,cofs,mons,vars;
   if nargs>=3 then vars:=Vars else vars:=indets(Monos) fi;
   mtab := table(['Monos[i]=i'$'i'=1..nops(Monos)]);
   A := matrix(nops(Polys),nops(Monos),0);
   for p to nops(Polys) do
      cofs := [coeffs(collect(Polys[p],vars,'distributed'),vars,'mons')];
      mons:=[mons];
      for i from 1 to nops(cofs) do
      	 A[p,mtab[mons[i]]] := cofs[i];
      od;
   od;
   op(A);
end:

#------------------------------------------------------------------------
# Classical Macaulay multiresultant of a list of n+1 dense inhomogeneous
# polynomials POLYS in n variables VARS (default: all variables in
# POLYS).  Returns a square matrix whose determinant vanishes: (i) if
# the polynomials have a common (projective) root; (ii) if they are
# insufficiently generic (e.g. if certain leading coefficients vanish).
# The matrix columns are labelled by monomials. Each row contains a
# monomial multiple of one of the input polynomials, entry
# 'column-monomial' being the coefficient of the multiple at that
# monomial, i.e. a coefficient of the original polynomial. Optional
# arguments MONOS,MULTS return the column monomials, row multiples.

macaulay := proc(Polys::list,Vars,Monos,Mults)
   local vars,degs,monos,mons, np,nv,nu,v,p, S,Sv, e_mults,s_mults;
   global `POLY`;
   if nargs>1 then vars:=Vars else vars:=sort([op(indets(Polys))]) fi;
   np := nops(Polys);
   nv := nops(vars);
   if np <> nv+1 then
      ERROR(`#polynomials must be #variables+1: #p = `.np.`, #v = `.nv);
   fi;

   # Macaulay degree calculation for dense polynomials

   degs := map((p,vars)->degree(p,vars),Polys,vars);
   nu := convert(degs,`+`) - nv;
   monos := [monomials(vars,0..nu)];

   # Macaulay row selection for a square resultant matrix. Arbitrarily
   # associate a polynomial p with each variable v, and work through the
   # variables.  Each variable `steals' all monomials of degree at least
   # (v)^degree(p) not yet stolen by an earlier one.  For purely
   # traditional reasons (the last variable is usually the one used for
   # homogenization), we work backwards and associate the first
   # polynomial with the last variable.

   S := {op(monos)};
   e_mults := []; 
   s_mults := [];
   for v from nv by -1 to 0 do
      p := np-v;
      if v>0 then
      	 Sv := {monomials(vars[v]^degs[p]*(1+convert(vars,`+`))^(nu-degs[p]),vars)}
	       intersect S;
      	 S := S minus Sv;
      	 mons := [op(sort(expand(convert(Sv,`+`)/vars[v]^degs[p]),vars,'tdeg'))];
      else
      	 mons := [op(sort(convert(S,`+`),vars,'tdeg'))];      
      fi;
      e_mults := [ op(e_mults), op(map(`*`,mons,Polys[p])) ];
      s_mults := [ op(s_mults), op(map(`*`,mons,`POLY`[p])) ];
   od;

   # Output the column monomials, row multiples and resultant matrix.

   if nargs>2 then Monos:=monos else lprint(nu,nops(monos),monos) fi;
   if nargs>3 then Mults:=s_mults else lprint(nops(s_mults),s_mults) fi;
   resultant_matrix(e_mults,monos,vars);
end:

#--------------------------------------------------------------------
# Rectangular Macaulay multiresultant matrix for dense inhomogeneous
# polynomials. The rows represent _all_ possible multiples of the given
# polynomials with terms of degree at most nu. For nu<=0, the routine
# chooses nu to be the traditional Macaulay multiresultant degree (or
# the degree required for the linear u-resultant if #polys = #vars). In
# this case, the traditional square Macaulay multiresultant matrix is a
# selection of the rows of this rectangular one, that is generically
# non-singular unless the system has a solution.  In any case, both
# matrices are rank deficient by at least the number of roots of the
# system.  But the rectangular matrix may have significantly better
# numerical condition than the square one, so it is recommended for
# numerical root calculations.

rect_macaulay := 
proc(Nu::integer,Polys::list,Vars,Monos,Mults)
   local nu,polys,vars,degs,monos,m, np,nv,p, e_mults,s_mults;
   global `POLY`;
   if nargs>2 then vars:=Vars else vars:=sort([op(indets(Polys))]) fi;
   polys := map((p,vars)->collect(p,vars,`distributed`),Polys,vars);
   degs := map((p,vars)->degree(p,vars),polys,vars);
   np := nops(polys);
   nv := nops(vars);
   if Nu>0 then 
      nu := Nu
   elif np < nv then 
      ERROR(`Multiresultant requires #polynomials >= #variables but #p=`.np.`, #v=`.nv);
   elif np = nv then
      nu := convert(degs,`+`) - nv + 1
   else
      nu := convert(sort(degs)[1..nv+1],`+`) - nv
   fi;
   monos := [monomials(vars,0..nu)];
   e_mults := []; 
   s_mults := [];
   for p from 1 to nops(polys) do
      if degs[p]>nu then 
         lprint(`Skipping polynomial `.p.` - degree `.(degs[p]).
	       ` exceeds hull degree `.nu)
      else
      	 for m in monomials(vars,0..nu-degs[p]) do
	    e_mults := [ op(e_mults), m * polys[p] ];
	    s_mults := [ op(s_mults), m * `POLY`[p] ];
      	 od;
      fi;
   od;
      
   if nargs>3 then Monos:=monos else lprint(nu,nops(monos),monos) fi;
   if nargs>4 then Mults:=s_mults else lprint(nops(s_mults),s_mults) fi;
   resultant_matrix(e_mults,monos,vars);
end: