## ## Title: Maple function for computing resultants of homogeneous polynomials ## Created: june 2010 ## Author: Laurent Buse ## Description: Sylvester, Bezout and Hybrid matrices for the Sylvester resultant ## Macaulay formula for the resultant of several homogeneous polnyomials ## with(LinearAlgebra): with(algcurves): #### Sylvester Matrix############################################## ## This function returns the Sylvester matrix of two polynomials ## f and g in the variable var seen of degree m and n respectively. ####################################################################### SylvMat:= proc(f,g,m,n,var) local i,j,M; M:=Matrix(m+n,m+n); for i from 1 to n do for j from 1 to m+n do M[j,i]:=coeftayl(var^(n-i)*f,var=0,m+n-j); od; od; for i from 1 to m do for j from 1 to m+n do M[j,i+n]:=coeftayl(var^(m-i)*g,var=0,m+n-j); od; od; return(M); end: ## This function resturns the resultant (the determinant of the previous matrix) ## It allows print out the "true" degrees of f and g SylvRes:= proc(f,g,m,n,var) print(degree(f,var),degree(g,var)); return(Determinant(SylvMat(f,g,m,n,var))); end: ## This function returns the discriminant Disc:= proc(f,m,var) print(degree(f,var)); return(simplify(expand(Determinant(SylvMat(f,diff(f,var),m,m-1,var)))/coeftayl(f,var=0,m))); end: #### Bezout Matrix################################################# ## This function returns the Bezout matrix of two polynomials ## f and g in the variable var seen as degree n polynomials. ####################################################################### BezMat:= proc(f,g,n,var) local i,j,b,M,_var; M:=Matrix(n,n); b:=simplify((f*subs(var=_var,g)-g*subs(var=_var,f))/(var-_var)); for i from 1 to n do for j from 1 to n do M[i,j]:=coeftayl(b,[var,_var]=[0,0],[n-i,n-j]); od; od; return(M); end: #### Hybrid Sylvester/Bezout Matrix ############################################# ## This function returns the hybrid Sylvester/Bezout matrix of two polynomials ## f and g in the variable var seen as degree d and e polynomials respectively. ################################################################################ HBezMat:= proc(f,g,d,e,var) local P,Q,m,n,i,j,M,A,B,C,D,K; if d=e then return(BezMat(f,g,d,var)); end if; if d>e then P:=f; Q:=g; m:=d; n:=e; else P:=g; Q:=f; m:=e; n:=d; end if; M:=Matrix(m,m); for i from 1 to m-n do ## colonnes du type Sylvester for j from 1 to m do M[j,i]:=coeftayl(var^(m-n-i)*Q,var=0,m-j); od; od; for i from 1 to n do ## colonnes du type Bezout A:=0; for j from 0 to n-i do A:=A+coeftayl(P,var=0,m-n+i+j)*var^(j); od; B:=0; for j from 0 to i-1 do B:=B+coeftayl(Q,var=0,j)*var^(m-n+j); od; C:=0; for j from 0 to m-n+i-1 do C:=C+coeftayl(P,var=0,j)*var^(j); od; D:=0; for j from 0 to n-i do D:=D+coeftayl(Q,var=0,j+i)*var^(j); od; K:=A*B-C*D; for j from 1 to m do M[j,m-n+i]:=coeftayl(K,var=0,m-j); od; od; return(M); end: #### Macaulay's Formula ############################################# ## This function gives the Macaulay formula for the resultant of ## several of homogeneous polynomials ###################################################################### #return the coeff. of a polynomial p in the variables var of the monomial m: coeffof := proc(m,p,var) local lm,lc,k; lc := [coeffs(p,var,'lm')]; lm := [lm]; if member(m,lm,'k') then lc[k] else 0 fi; end: ##lp is a list of homogeneous polynomials and var is a list of variables MacMatrix:=proc (lp::list,var::list) local i,j,k,listmon,N,n,ld,nu,t,dodu,Base,listpol,tt,M,s; n:=nops(lp); #Test:nombre de polynome=nombre de variables if(nops(var)<>n) then ERROR(`the number of polynomials in the list must be the number of variables +1 `); fi; #Calcul des degres, liste ld, puis de nu=sum(ld)-n+1 ld := NULL; for i to n do ld := ld,degree(lp[i],var); od; ld := [ld]; nu := convert(ld,`+`)-n+1; #List de tous les monomes unassign('s'); coeffs(expand((convert(var,`+`))^(nu)),var,'s'); listmon:={s}; dodu:=[NULL]; ##liste d'entiers qui donne la sous-matrice Base:=[NULL]; ##liste des monomes qui indexent la matrice listpol:=[NULL]; ##la liste des polynomes de la matrice for i to n-1 do #i correspond a la variable var[i] unassign('s'); coeffs(expand((convert(var,`+`))^(nu-ld[i])),var,'s'); t:=[s]; tt:=NULL; for j to nops(t) do tt:=tt,(t[j]*var[i]^(ld[i])); od; tt:=[op({tt} intersect listmon)]; for j to nops(tt) do listpol:=[op(listpol),expand((lp[i])*tt[j]/(var[i]^(ld[i])))]; for k to n do if ( (k<>i) and (degree(tt[j],var[k])>=ld[k]) ) then dodu:=[op(dodu),j+nops(Base)]; fi; od; od; listmon:=listmon minus {op(tt)}; Base:=[op(Base),op(tt)]; od; tt:=[op(listmon)]; for j from 1 to nops(tt) do listpol:=[op(listpol),expand((lp[n])*tt[j]/var[n]^(ld[n]))]; od; Base:=[op(Base),op(tt)]; dodu:=[op({op(dodu)})]; #On construit la matrice de Macaulay: N:=nops(listpol); M:=Matrix(N); for i to N do for j to N do M[j,i]:=coeffof(Base[j],listpol[i],var); od; od; #print(Base); #print(dodu); return(M,SubMatrix(M,dodu,dodu)); end: ##Here is an example ##M,N:=MacMatrix([a0*x^2+a1*x*y+a2*x*z+a3*y^2+a4*y*z+a5*z^2, ##b0*x^2+b1*x*y+b2*x*z+b3*y^2+b4*y*z+b5*z^2, ##c0*x^2+c1*x*y+c2*x*z+c3*y^2+c4*y*z+c5*z^2],[x,y,z]); ##factor(Determinant(M)); ##factor(Determinant(N)); ##lp is a list of homogeneous polynomials and var is a list of variables, nu (the "critical" degree in which the Macaulay formula is built) is now given as a parameter MacMatrixGen:=proc (lp::list,var::list,nu) local i,j,k,listmon,N,n,ld,t,dodu,Base,listpol,tt,M,s; n:=nops(lp); #Test:nombre de polynome=nombre de variables if(nops(var)<>n) then ERROR(`the number of polynomials in the list must be the number of variables +1 `); fi; #Calcul des degres, liste ld, puis de nu=sum(ld)-n+1 ld := NULL; for i to n do ld := ld,degree(lp[i],var); od; ld := [ld]; # nu := convert(ld,`+`)-n+1; #List de tous les monomes unassign('s'); coeffs(expand((convert(var,`+`))^(nu)),var,'s'); listmon:={s}; dodu:=[NULL]; ##liste d'entiers qui donne la sous-matrice Base:=[NULL]; ##liste des monomes qui indexent la matrice listpol:=[NULL]; ##la liste des polynomes de la matrice for i to n-1 do #i correspond a la variable var[i] unassign('s'); coeffs(expand((convert(var,`+`))^(nu-ld[i])),var,'s'); t:=[s]; tt:=NULL; for j to nops(t) do tt:=tt,(t[j]*var[i]^(ld[i])); od; tt:=[op({tt} intersect listmon)]; for j to nops(tt) do listpol:=[op(listpol),expand((lp[i])*tt[j]/(var[i]^(ld[i])))]; for k to n do if ( (k<>i) and (degree(tt[j],var[k])>=ld[k]) ) then dodu:=[op(dodu),j+nops(Base)]; fi; od; od; listmon:=listmon minus {op(tt)}; Base:=[op(Base),op(tt)]; od; tt:=[op(listmon)]; for j from 1 to nops(tt) do listpol:=[op(listpol),expand((lp[n])*tt[j]/var[n]^(ld[n]))]; od; Base:=[op(Base),op(tt)]; dodu:=[op({op(dodu)})]; #On construit la matrice de Macaulay: N:=nops(listpol); M:=Matrix(N); for i to N do for j to N do M[j,i]:=coeffof(Base[j],listpol[i],var); od; od; #print(Base); #print(dodu); return(M,SubMatrix(M,dodu,dodu)); end: ## Get the resultant via determinant computation MacRes:=proc (lp::list,var::list) local M,N; M,N:=MacMatrix(lp,var); return(simplify(Determinant(M)/Determinant(N))); end: ## A simple call is the case of three polynomials given as affine polynomials ## in two variables. The homogeneization is done in the projective plane. BiMacRes:=proc (f1,f2,f3,var1,var2) local F1,F2,F3,_hvar; F1:=homogeneous(f1,var1,var2,_hvar); F2:=homogeneous(f2,var1,var2,_hvar); F3:=homogeneous(f3,var1,var2,_hvar); return(MacRes([F1,F2,F3],[var1,var2,_hvar])); end: