# Maple VII code for bivariate sparse resultants

# Author: Amit Khetan, University of Massachusetts at Amherst
# address@hidden, http://www.math.umass.edu/~khetan

# Created: 3/2002
# Last update: 10/2003

# Main function call:
# bires([poly1, poly2, poly3], var1, var2)
#  
# RETURNS unmixed sparse resultant matrix on the union of the supports


# Reference:
# Amit Khetan, "The resultant of an unmixed bivariate system", 
# Journal of Symbolic Computation 36 (2003) 425-442

with(linalg);
with(geometry);

################ utilities #########################################


supp_mat := proc (supp::list, vars::list)

#INPUT: list of monomials
#OUTPUT: matrix of exponent vectors of support

  local m, v, nmons, nvars, suppmat;
  
  nmons := nops(supp);
  nvars := nops(vars);
  nvars := nops (vars);
  suppmat := matrix(nmons, nvars);
  for m from 1 to nmons do
    for v from 1 to nvars do suppmat[m,v] := degree (supp[m], vars[v]); od;
  od;	
  evalm (suppmat);

end: # supp_mat


convhull := proc(suppmat, nmons)

#INPUT: support matrix
#OUTPUT: list of vertices in convex hull in counterclockwise order

local pts, m, hull, vertlist, v;

for m from 1 to nmons do
   pts[m] := point(p||m, convert(row(suppmat, m), list));
od;
hull := convexhull(convert(pts, list));

for v from 1 to nops(hull) do
   vertlist[v] := coordinates(hull[v]);
od;

eval(convert(vertlist, list));
end:



getpoints := proc(hull::list)

# The "polygon fill algorithm" for computing lattice points in a convex
# polygon

#INPUT: list of vertices in counterclockwise order
#OUTPUT: list of all lattice points in polygon

#RETURNS two lists: (points, interiorpoints);

local nvert, miny, minindex, maxy, maxindex, v, left, leftstart,
rightstart, leftend, rightend, leftx, rightx, pts, intpts, numpts,
numintpts, right, y, x, nextleft, nextright, leftslope, rightslope;

nvert := nops(hull); miny := 32767; maxy := -32767;
for v from 1 to nvert do
   if ((hull[v])[2] < miny) then
       miny := (hull[v])[2];
   minindex := v;
   fi;
   if ((hull[v])[2] > maxy) then
     maxy := (hull[v])[2];
     maxindex := v;
   fi;
od;


if  minindex = 1 then
  if (hull[nvert])[2] = miny then
      leftstart := nvert;
      rightstart := minindex;
  else
      if (hull[minindex+1])[2] = miny then
          leftstart := minindex;
          rightstart := minindex+1;
      else
      	leftstart := minindex;
      	rightstart := minindex;
      fi;
  fi;
else
  if (minindex = nvert) then
     leftstart := minindex;
     rightstart := minindex;
  else	
	if (hull[minindex+1])[2] = miny then
      		leftstart := minindex;	
      		rightstart := minindex+1;
	else
		leftstart := minindex;
      		rightstart := minindex;
	fi;
  fi;
fi;

v := rightstart;

while (hull[v][2] < maxy) do
  v := v + 1;
  if v > nvert then v := v - nvert; fi;
od;

rightend := v;

if (v < nvert) then
   if ((hull[v+1])[2] = maxy) then
      leftend := v+1;
   else
      leftend := v;
   fi;
else
   if ((hull[1])[2] = maxy) then
      leftend := 1;
   else
      leftend := v;
   fi;
fi;


v := leftstart;

left := leftstart;
nextleft := left - 1;
if (nextleft = 0) then nextleft := nvert; fi;
right := rightstart;
nextright := right + 1;
if (nextright > nvert) then nextright := 1; fi;

leftslope := (hull[nextleft][1] - hull[left][1])/(hull[nextleft][2] - hull[left][2]); 
rightslope := (hull[nextright][1] -hull[right][1])/(hull[nextright][2] - hull[right][2]);

numpts := 0;
numintpts := 0;
leftx := hull[leftstart][1];
rightx := hull[rightstart][1];
for y from miny to maxy do
  for x from ceil(leftx) to floor(rightx) do
      numpts := numpts + 1;
      pts[numpts] := [x, y];
      
      if ((x > leftx) and (x < rightx) and (y > miny) and (y < maxy)) then
	  numintpts := numintpts + 1;
          intpts[numintpts] := [x, y];
      fi;
  od;
  if (y < maxy) then
  	if (y = hull[nextleft][2]) then
  	   left := nextleft;
     	   nextleft := left - 1;
           if  (nextleft = 0) then nextleft := nvert; fi;
     	   leftslope := (hull[nextleft][1] - hull[left][1])/(hull[nextleft][2] - hull[left][2]);
  	fi;
 	if (y = hull[nextright][2]) then
    	   right := nextright;
    	   nextright := right + 1;
           if  (nextright > nvert) then nextright := 1; fi;
           rightslope := (hull[nextright][1] - hull[right][1])/(hull[nextright][2] - hull[right][2]);
        fi;
        leftx := leftx + leftslope;
        rightx := rightx + rightslope;
  fi;            
od;
if (numintpts = 0) then
   intpts := []; fi;
return[convert(pts, list), convert(intpts, list)];

end:


get_fan := proc(hull::list)

# Computes the normal fan of a polygon and corresponding linear functional

# INPUT: List of vertices in counterclockwise order
# OUTPUT:[[fan[1], ..., fan[s]], [a[1], ..., a[s]]] 
#        where fan[i] dot [x,y] >= -a[i] are the defining inequalities 


local nvert, fan, a, i, v, n, d;

nvert := nops(hull);
fan := array(1..nvert);
a := array(1..nvert);
for i from 1 to nvert-1 do
     v := (hull[i]-hull[i+1]);
     n := v[2];
     d := -v[1];
     n := n/gcd(n,d);
     d := d/gcd(v[2],d);
     fan[i] := [n,d];
     a[i] := -dotprod(vector(hull[i]), vector(fan[i]));
   od;
   v := (hull[nvert]-hull[1]);
   n := v[2]; 
   d := -v[1];
   n := n/gcd(n,d);
   d := d/gcd(v[2],d);
   fan[nvert] := [n,d];
   a[nvert] := -dotprod(vector(hull[nvert]), Vector(fan[nvert]));
   return(convert(eval(fan), list), convert(eval(a), list));
end:

get_homog := proc(pt, fan, a)

#Computes the "homogenous coordinates" (edge lattice distances) of point
#with respect to a polygon

#INPUT: The point, the normal fan, and the linear functional
#OUTPUT: The list [<pt, fan[1]> + a[1], ..., <pt, fan[s]> + a[s]]

local v1;

v1 := map(v->dotprod(pt, vector(v)), fan);
v1 := v1 + a;

return(convert(v1, list));
end: 

add_supports := proc(supp1, supp2)

# Pointwise sum of two supports

local size1, size2, i, j, S;

size1 := nops(supp1);
size2 := nops(supp2);
S := {};
for i from 1 to size1 do
    for j from 1 to size2 do
        S := S union {supp1[i] + supp2[j]};
    od;
od;
return convert(S, list);

end:

get_interior := proc(pts)

# INPUT: A list of points in homogeneous coordinates
# OUTPUT: The sublist of interior points (all coordinates at least one)

local numpts, size, interior, i, newpt, d;

interior := [];
numpts := nops(pts);
if (numpts = 0) then
  return interior
else
   size := nops(pts[1]);
   for i from 1 to size do
     d[i] := 1;
   od;
   d := convert(d, list);
   for i from 1 to numpts do
     if (verify(pts[i], d, 'list'('greater_equal'))) then
        interior := [op(interior), pts[i] - d];
     fi;
   od;
   return interior;
fi;
end:

###################################################################
###  Main Procedure 						###
###################################################################


bires := proc(f, s, t)

#INPUT: a list of 3 polynomials and two variables
#OUTPUT: resultant matrix

local f0mons, f1mons, f2mons, suppset, supp, C, i, j, supppts,
suppmat, nmons, hull,F, fan, a, P, pts, intpts, qpts, inta, intqpts,
twoqpts, int2qpts, maxang, ang, eta_1, eta_2, M, R_1, R_2, R_3, c, B,
e, m, n, alpha, w1, w, flag, v1, v,u, u1, beta, g, size, rows, cols;           ;

# input:
# x = f0, y = f2, z = f3 3 polynomials in two variables s,t

# Compute monomial supports


coeffs(collect(expand(f[1]), [s,t], distributed), [s, t], 'f0mons');
coeffs(collect(expand(f[2]), [s,t], distributed), [s, t], 'f1mons');
coeffs(collect(expand(f[3]), [s,t], distributed), [s, t], 'f2mons');
suppset := `union`({f0mons}, {f1mons},  {f2mons});
supp := convert(sort(suppset), list);

# Get coefficient matrix

C := matrix(3, nops(supp));
for i from 1 to 3 do
   for j from 1 to nops(supp) do
       C[i, j] := coeftayl(f[i], [s,t]=[0,0], [degree(supp[j], s), degree(supp[j], t)]);
   od;
od;


# Compute convex hull of support and lattice points/interior lattice points

supppts := [];
suppmat := supp_mat(supp, [s,t]);
nmons := nops(supp);
for i from 1 to nmons do supppts := [op(supppts), row(suppmat, i)];
od;
hull := convhull(suppmat, nmons);
F := get_fan(hull);
fan := F[1];
a := F[2];
P := getpoints(hull);
pts := P[1];
intpts := P[2];

# Convert all to homogeneous coordinates

qpts := map(m->get_homog(m, fan, a), pts);
inta := map(x->x-1, a);
intqpts := map(m->get_homog(m, fan, inta), intpts);
twoqpts := add_supports(qpts, qpts);
int2qpts := get_interior(twoqpts);
supppts := map(m->get_homog(m, fan, a), supppts);

# Partition the fan

maxang := 0;
   for i from 1 to nops(fan)-1 do
	ang := angle(fan[i], fan[i+1]);
	if (evalf(ang) > evalf(maxang)) then
	   maxang := ang;
	   eta_1 := i;
	   eta_2 := i+1;
	fi;
    od;
    ang := angle(fan[nops(fan)], fan[1]);
    if (evalf(ang) > evalf(maxang)) then
       maxang := ang;
       eta_1 := nops(fan);
       eta_2 := 1;
    fi;
	

   M := inverse(matrix([fan[eta_1], fan[eta_2]]));
   R_1 := {};
   R_2 := {};
   R_3 := {};
  
   for i from 1 to nops(fan) do
       c := multiply(vector(fan[i]), M);
       if (c[1] >= 0 and c[2] <= 0) then
	  R_1 := R_1 union {i}
	else
	  if (c[2] >= 0) then
	     R_2 := R_2 union {i};
	  else
	     R_3 := R_3 union {i};
	  fi;
	fi;
   od;
  if (R_3 = {}) then R_3 := {0} fi;

# Fill in Sylvester rows and columns

  B := matrix(nops(int2qpts) + 3, nops(qpts) + 3*nops(intqpts), 0);
  e := vector(nops(fan), 1);
  e := convert(e, list);
  for i from 1 to nmons do
     m := supppts[i];
     member(m, qpts, 'cl');
     for j from 1 to 3 do	
        B[nops(int2qpts) + j, cl] := C[j,i];
        for n in intqpts do
            member(n+m, int2qpts, 'rw');
            member(n, intqpts, 'cl1');
            B[rw, nops(qpts) + (j-1)*nops(intqpts) + cl1] := C[j,i];
        od;
     od;
  od; 

# Fill in Bezout matrix by finding all triples (w,v,u) satisfying conditions

  for alpha in qpts do
    member(alpha, qpts, 'cl');
    for w1 from 1 to nmons do
        w := supppts[w1];
        flag := false;
	for i in R_3 do
	if (i = 0) then
	   if ((w[1] + w[2]) >= (alpha[1] + alpha[2])) then
              flag := true;
              break;
           fi;
        else
           if (w[i] <= alpha[i]) then
               flag := true;
	       break;
	   fi;
        fi;
	od;
	    
	if (flag) then next; fi; 
	    for i in R_2 do
		if (w[i] <= alpha[i]) then
		    flag := true;
		    break;
		fi;
	    od;
	    if (not(flag)) then next; fi;
	
            for v1 from 1 to nmons do
 	    v := supppts[v1];
	    flag := false;
	    if (v1 = w1) then next; fi;
	    c := v + w;
	    for i in R_2 do
	        if (c[i] <= alpha[i]) then
		   flag := true;
		   break;
		fi;
	    od;
	    if (flag) then next; fi;
	    for i in R_1 do;
		if (c[i] <= alpha[i]) then 	
  	     	   flag := true;
		   break;
		fi;
	    od;
	if (not(flag)) then next; fi;
      	for u1 from 1 to nmons do	
            u := supppts[u1];
	    flag := false;
            if ((u1=v1) or (u1=w1)) then next; fi;
	    c := u + v + w;
	    for i in R_1 do;
		if (c[i] <= alpha[i]) then
		   flag := true;
		   break;
		fi;
	    od;
	    if (flag) then next; fi;
	    beta := c - alpha - e;
	    member(beta, int2qpts, 'rw');
	   
           
            g := det(submatrix(C, [1,2,3], [u1, v1, w1]));
	    	    B[rw,cl] := B[rw,cl] + g;
        od;
 od;
od;
od;

size := nops(int2qpts) + 3;
rows  := []:
cols := []:

for i from 1 to size do

if (i <= 3) then rows := [op(rows), i+nops(int2qpts)]: else
rows := [op(rows), i-3]: fi:
if (i <= 3*nops(intqpts)) then cols := [op(cols), i + nops(qpts)]: else
cols := [op(cols), i - 3*nops(intqpts)]: fi:
od:
M := submatrix(B, rows, cols):

return(evalm(M)):
 
  	
end: