#--------------------------------------------------------------------------------------------------------
# author : Nan Li
# date   : Mar. 2nd 2010
#--------------------------------------------------------------------------------------------------------

#####   Differential Function acts on a Polynomial System at Zero  ######################################
##   
##	Calling sequence:  DifOnPol(lsys,ord,dualsupport)
##      
##	Input:  An polynomial system, an order and a dualsupport
##
##	Output:  The polynomial after differential operator acts on the polynomial system
##
#########################################################################################################


DifOnPol:=proc(lsys,ord,dualsupport)

	local L,v,dualsupp,coefficient,i,dif,j,k,u;

	v:=Vector(nops(lsys));
	dualsupp:=dualsupport+1/ord[1];
	if nops(dualsupp)>2 then
		for i from 1 to nops(dualsupport) do
			L:=[];
			coefficient:=1;
			dif:=Conv(op(dualsupport)[i],ord);
			for j from 1 to nops(ord) do
				for k from 1 to dif[j] do
					L:=[op(L),ord[j]];
				od;
				coefficient:=coefficient*dif[j]!;
			od;
			for u from 1 to nops(lsys) do
				v[u]:=v[u]+coeffs(op(dualsupport)[i],ord)*diff(lsys[u],L)/coefficient;
			od;
		od;
	else
		L:=[];
		coefficient:=1;
		dif:=Conv(dualsupport,ord);
		for j from 1 to nops(ord) do
			for k from 1 to dif[j] do
				L:=[op(L),ord[j]];
			od;
			coefficient:=coefficient*dif[j]!;
		od;
		for u from 1 to nops(lsys) do
			v[u]:=v[u]+coeffs(dualsupport,ord)*diff(lsys[u],L)/coefficient;
		od;
	fi;

	RETURN(v);

end proc:
	
		
#####   Evalue Vector   #################################################################################
##   
##	Calling sequence:  EvalVector(ord)
##      
##	Input:  An order
##
##	Output:  EvalVector at the zero
##
#########################################################################################################


EvalVector:=proc(ord)

	local L,i;

	L:=[];
	for i from 1 to nops(ord) do
		L:=[op(L),ord[i]=0];
	od;

	RETURN(L);

end proc:
	
		
#####   Differential Operator   #########################################################################
##   
##	Calling sequence:  DifOperator(dualsupp,ord,i)
##      
##	Input:  A dualsupport, an order and ith
##
##	Output:  The differential functional
##
#########################################################################################################


DifOperator:=proc(dualsupp,ord,i)


	RETURN(expand(dualsupp*ord[i]));

end proc:
	
		
#####   Convert to Differential Functional   ############################################################
##   
##	Calling sequence:  Conv(mon,ord)
##      
##	Input:  A monomial and an order
##
##	Output:  The Differential Functional respect to the monomial
##
#########################################################################################################


Conv:=proc(mon,ord)

	local L,i;

	L:=[];
	for i from 1 to nops(ord) do
		L:=[op(L),degree(mon,ord[i])];
	od;

	RETURN(L);

end proc:
	
		
#####   Transform A Polynomial System Based an A Matrix   ###############################################
##   
##	Calling sequence:  TransPol(lsys,TransMatrix,ord)
##      
##	Input:  A polynomial system, a matrix and an order
##
##	Output:  The transformed system
##
#########################################################################################################


TransPol:=proc(lsys,TransMatrix,ord)

	local i,L,v;

	L:=[];
	v:=TransMatrix.convert(ord,Vector);
	for i from 1 to nops(ord) do
		L:=[op(L),ord[i]=v[i]];
	od;

	RETURN(expand(subs(L,lsys)));

end proc:
	
		
#####   Form A Full-Rank Matrix   #######################################################################
##   
##	Calling sequence:  FullRank(nullvector)
##      
##	Input:  A vector
##
##	Output:  A full-rank matrix based on the vector 
##
#########################################################################################################


FullRank:=proc(nullvector)

	local ns,i,A;

	ns:=NullSpace(Transpose(convert(nullvector,Matrix)));
	A:=convert(nullvector,Matrix);
	for i from 1 to nops(ns) do
		A:=Matrix([A,convert(ns[i],Matrix)]);
	od;

	RETURN(A);
	
end proc:
	
		
#####   Separate the Differential Functional into Part   ################################################
##   
##	Calling sequence:  Part(rest,ord,j)
##      
##	Input:  A differential functional, an order and jth
##
##	Output:  The part of the differential functional
##
#########################################################################################################


Part:=proc(rest,ord,j)

	local L,k,r;

	L:=0;
	r:=rest+1/ord[1];
	if nops(r)>2 then
		for k from 1 to nops(rest) do
			if degree(op(rest)[k],ord[j])=0 then
				L:=L+op(rest)[k];
			fi;
		od;
	else
		if degree(rest,ord[j])=0 then
			L:=L+rest;
		fi;
	fi;

	RETURN(expand(L));

end proc:
	
		
#####   Compute the Coefficient Before A Differential Operator   ########################################
##   
##	Calling sequence:  MultiFactorial(ind)
##      
##	Input:  An index
##
##	Output:  The coefficient
##
#########################################################################################################


MultiFactorial:=proc(ind)

	local mf,i;

	mf:=1;
	for i from 1 to nops(ind) do
		mf:=mf*(ind[i]!);
	od;

	RETURN(mf);
	
end proc:
	
		
#####   Transform the Dual Basis Back to Original System   ##############################################
##   
##	Calling sequence:  TransBack(dualbasis,TransMatrix,ord)
##      
##	Input:  A dualbasis, a matrix and an order
##
##	Output:  The dualbasis respect to the original system
##
#########################################################################################################


TransBack:=proc(dualbasis,TransMatrix,ord)

	local i,j,k,DualBasis,poly,tmid,mon,mid,ttmid;

	DualBasis:=[];
	for i from 1 to nops(dualbasis) do
		poly:=0;
		tmid:=dualbasis[i]+1/ord[1];
		if nops(tmid)>2 then
			for j from 1 to nops(dualbasis[i]) do
				mon:=0;
				mid:=TransPol(op(dualbasis[i])[j],Transpose(TransMatrix),ord);
				ttmid:=mid+1/ord[1];
				if nops(ttmid)>2 then
					for k from 1 to nops(mid) do
						mon:=mon+MultiFactorial(Conv(op(mid)[k],ord))*op(mid)[k];
					od;
				else
					mon:=mon+MultiFactorial(Conv(mid,ord))*mid;
				fi;
				poly:=poly+mon/MultiFactorial(Conv(op(dualbasis[i])[j],ord));
			od;			
		else
			mon:=0;
			mid:=TransPol(dualbasis[i],Transpose(TransMatrix),ord);
			ttmid:=mid+1/ord[1];
			if nops(ttmid)>2 then
				for k from 1 to nops(mid) do
					mon:=mon+MultiFactorial(Conv(op(mid)[k],ord))*op(mid)[k];
				od;
			else
				mon:=mon+MultiFactorial(Conv(mid,ord))*mid;
			fi;
			poly:=poly+mon/MultiFactorial(Conv(dualbasis[i],ord));
		fi;
		DualBasis:=[op(DualBasis),poly];
	od;

	RETURN(DualBasis);
	
end proc:
	
		
#####   Refine the Mutiple Root   #######################################################################
##   
##	Calling sequence:  RootRefine(lsys,ord,appx,tol)
##      
##	Input:  One polynomial system, an order, an approxiate root and a tolerance
##
##	Output:  The multiplity, the index and the refined root 
##
#########################################################################################################


MultiStructure:=proc(pol,ord,mroot)

	local lsys,E,F,v,s,t,Jac,nullvector,TransMatrix,dualbasis,DualBasis,
	adddual,HilbertFunction,depth,multiplicity,involutive,dualsupp,DualSupp,
	DualSupport,P,L,U,i,j,CoeffMatrix,rest,a,b,p;

	lsys:=pol;
	F:=[];
	for i from 1 to nops(ord) do
		F:=[op(F),ord[i]=mroot[i]];
	od;
	Jac:=eval(VectorCalculus[Jacobian](lsys,ord),F);
	nullvector:=NullSpace(Jac)[1];
	TransMatrix:=FullRank(nullvector);
	userinfo(2,MultiStructure,print(`The Transpose Matrix is:`,TransMatrix));
	lsys:=TransPol(lsys,TransMatrix,ord);
	userinfo(2,MultiStructure,print(`The New Polynomial is:`,lsys));
	v:=MatrixInverse(TransMatrix).mroot;
	E:=[];
	for i from 1 to nops(ord) do
		E:=[op(E),ord[i]=v[i]];
	od;	
	s:=nops(ord);
	t:=nops(lsys);
	Jac:=eval(VectorCalculus[Jacobian](lsys,ord),E);
	P,L,U:=LinearAlgebra[LUDecomposition](SubMatrix(Jac,1..-1,2..-1));
	userinfo(2,MultiStructure,print(`The LUDecomposition is:`,P,L,U));
	DualSupp:=[ord[2]];
	for i from 3 to s do
		DualSupp:=[op(DualSupp),ord[i]];
	od;
	dualbasis:=[1,ord[1]];
	adddual:=ord[1];
	HilbertFunction:=[1,1];
	depth:=1;
	multiplicity:=2;
	involutive:=0;
	while involutive=0 do
		dualsupp:=DifOperator(adddual,ord,1);
		if multiplicity>3 then
			for i from multiplicity-1 by -1 to 2 do
				rest:=Part(dualbasis[i],ord,1);
				for j from 2 to s do
					dualsupp:=dualsupp+eval(diff(dualbasis[multiplicity-i+2],ord[j]),EvalVector(ord))*DifOperator(rest,ord,j);
					rest:=Part(rest,ord,j);
				od;
			od;
		fi;
		userinfo(3,MultiStructure,print(`The Dual Support is:`,dualsupp));
		p:=-eval(DifOnPol(lsys,ord,dualsupp),E);
		b:=LinearAlgebra[LinearSolve] (L,Transpose(P).p,method='subs');
		if b[-1]=0 then
			a:=LinearAlgebra[LinearSolve] (U,b,method='subs');
			adddual:=dualsupp+(Transpose(convert(a,Matrix)).convert(DualSupp,Vector))[1];
			dualbasis:=[op(dualbasis),adddual];
			depth:=depth+1;
			multiplicity:=multiplicity+1;
		else
			involutive:=1;
		fi;
	od;
	DualBasis:=TransBack(dualbasis,TransMatrix,ord);
	userinfo(2,MultiStructure,print(`The dual basis of new polynomial at isolated zero is:`,dualbasis));
	userinfo(1,MultiStructure,print(`The multiplicity is:`,multiplicity));
	userinfo(2,MultiStructure,print(`The depth is:`,depth));
	userinfo(1,MultiStructure,print(`The dual basis of original polynomial at isolated zero is:`,DualBasis));
	for i from 1 to multiplicity do
		userinfo(1,MultiStructure,print(`The dual basis `,i,` acts on original polynomial at isolated zero is:`,map(expand,eval(DifOnPol(pol,ord,DualBasis[i]),F))));
	od;

end proc: