##
##    Title: 	Maple function for computing resultants of homogeneous polynomials
##    Created:	june 2010
##    Author: 	Laurent Buse <Laurent.Buse@inria.fr>
## 	  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: