### Procedures

 ## computes the list of exponents of the monomial mon
 get_explist := proc(mon, vars)
   local explist,i;
 explist := []:
 for i from 1 to nops(vars) do
   explist := [op(explist), degree(mon, vars[i])]:
 od:
 explist;
 end proc:

 ## computes the list of exponent vectors of the polynomial f
 get_supp := proc(f, vars)
    local mons, supp,i;
 coeffs(f, vars, 'mons'):
 mons := [mons]:
 supp := []:
 for i from 1 to nops(mons) do
   supp := [op(supp), get_explist(mons[i], vars)]:
 od:
 supp:
 end proc:

 ## computes the inner normals to the polygon given by points
 get_normals:=proc(points)
   local i,sled,l,v;
 v:=[]:
 for i from 1 to nops(points) do
   if i=nops(points) then sled:=1: else sled:=i+1: fi:
   l:=igcd(points[i][2]-points[sled][2],points[i][1]-points[sled][1]):
   v:=[op(v),[(points[i][2]-points[sled][2])/l,
       (points[sled][1]-points[i][1])/l]]:
   od:
 v:
 end proc:


### Write down f1 and f2

 f1:=0: f2:=0:
 for i from 1 to nops(exp_f1) do
  f1:=f1+(a[i])*x^(exp_f1[i][1])*y^(exp_f1[i][2]):
 od:
 for i from 1 to nops(exp_f2) do
  f2:=f2+(b[i])*x^(exp_f2[i][1])*y^(exp_f2[i][2]):
 od:


### Compute the Minkowski sum P and its normals

 supp:=get_supp(collect(f1*f2,[x,y],'distributed'),[x,y]):
 supp:=convexhull({op(supp)}):
 s:=nops(supp):
 normals:=transpose(matrix(get_normals(supp))):


### Compute the inequalities for P1, P2 and P

 v_b:=[]: v_b1:=[]: v_b2:=[]:
 for i from 1 to s do
  v_b:=[op(v_b),-min(seq(dotprod(supp[j],[normals[1,i],normals[2,i]]),j=1..s))]:
  v_b1:=[op(v_b1),-min(seq(dotprod(exp_f1[j],[normals[1,i],normals[2,i]]),j=1..nops(exp_f1)))]:
  v_b2:=[op(v_b2),-min(seq(dotprod(exp_f2[j],[normals[1,i],normals[2,i]]),j=1..nops(exp_f2)))]:
 od:


### Homogenize f1 and f2

 exp_F1:=[seq(evalm(exp_f1[j]&*normals+v_b1),j=1..nops(exp_f1))]:
 F1:=0:
 for i from 1 to nops(exp_f1) do
 F1:=F1+coeftayl(f1,[x,y]=[0,0],exp_f1[i])*product((x[j])^exp_F1[i][j],j=1..s):
 od:

 exp_F2:=[seq(evalm(exp_f2[j]&*normals+v_b2),j=1..nops(exp_f2))]:
 F2:=0:
 for i from 1 to nops(exp_f2) do
 F2:=F2+coeftayl(f2,[x,y]=[0,0],exp_f2[i])*product((x[j])^exp_F2[i][j],j=1..s):
 od:


### Homogenize the monomial h and compute F0

 exp_H:=convert(evalm(monom&*normals+v_b),list):
 H:=product((x[j])^exp_H[j],j=1..s):

 mult:=1:
 for i from 1 to nops(exp_H) do
   if exp_H[i]<0 then mult:=mult*x[i]^(-exp_H[i]): fi:
 od:

 F0:=mult*product(x[j],j=1..s):


### Homogenize the toric jacobian w.r.t. the polygon P

 JMat:=matrix([[diff(f1,x),diff(f1,y)],[diff(f2,x),diff(f2,y)]]):
 jac:=expand(x*y*det(JMat)):
 exp_jac:=get_supp(jac, [x,y]):
 exp_J:=[seq(evalm(exp_jac[j]&*normals+v_b),j=1..nops(exp_jac))]:
 J:=0:
 for i from 1 to nops(exp_jac) do
   term:=product((x[j])^exp_J[i][j],j=1..s): all:=product(x[j],j=1..s):
   if divide(term,all)=false then
     J:=J+coeftayl(jac,[x,y]=[0,0],exp_jac[i])*term:
   fi:
 od:


### Compute the normal forms of mult*J and mult*H and obtain the residue r

 T:=tdeg(seq(x[j],j=1..s)):
 G:=gbasis([F0,F1,F2],T): #seq(leadterm(G[i],T),i=1..nops(G)):

 cJ:=normalf(mult*J,G,T):

 cH:=normalf(mult*H,G,T):

 r:=factor(simplify(cH/cJ)):