#--------------------------------------------------------------------------------------------------------
# author : Nan Li
# date   : Dec. 18 2009
#--------------------------------------------------------------------------------------------------------

#####   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,appx)
##      
##	Input:  An order and an approximate root
##
##	Output:  EvalVector at the approximate root
##
#########################################################################################################


EvalVector:=proc(ord,appx)

	local L,i;
	
	L:=[];
	for i from 1 to nops(ord) do
		L:=[op(L),ord[i]=appx[i]];
	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:
	
		
#####   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
##
#########################################################################################################


TransPol:=proc(lsys,TransMatrix,ord)

	local i,L,v;

	v:=TransMatrix.convert(ord,Vector);
	L:=EvalVector(ord,v);

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

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)


	RETURN(expand(subs([ord[j]=0],rest)));

end proc:
	
		
#####      ################################################
##   
##	Calling sequence:  
##      
##	Input:  
##
##	Output:  
##
#########################################################################################################


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:
	

#####      ################################################
##   
##	Calling sequence:  
##      
##	Input:  
##
##	Output:  
##
#########################################################################################################		


MDifOnPol:=proc(pol,ord,mon)

	local g,dif,i;

	g:=pol;
	dif:=Conv(mon,ord);
	for i from 1 to nops(ord) do
		g:=coeff(g,ord[i],dif[i]);
	od;

	RETURN(coeffs(mon)*g);

end proc:
	
		
#####   Refine the Mutiple Root   #######################################################################
##   
##	Calling sequence:  MultipleRootRefinerBreadthOne(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 
##
#########################################################################################################


MultipleRootRefinerBreadthOne:=proc(lsys,ord,appx,tol)

	local n,Ev,appy,appz,nlsys,nv,Jac,Jaclhs,nullvector,TransMatrix,dualbasis,adddual,multiplicity,flag,dualsupp,dualsuppx,NewMatrix,i,j,rest;

	n:=nops(ord);
	Ev:=[seq(ord[i]=appx[i],i=1..n)]:
	Jac:=VectorCalculus[Jacobian](lsys,ord);
	Jaclhs:=eval(Jac,Ev):
	appy:=appx+LinearAlgebra[LeastSquares](HermitianTranspose(Jaclhs).Jaclhs+LinearAlgebra[SingularValues](Jaclhs)[n]*IdentityMatrix(n),-HermitianTranspose(Jaclhs).convert(eval(lsys,Ev),Vector),method='SVD'):
	userinfo(3,MultipleRootRefinerBreadthOne,print(`The new approximation:`,appy));
 	Ev:=[seq(ord[i]=appy[i],i=1..n)]:
 	Jaclhs:=eval(Jac,Ev):
	TransMatrix:=LinearAlgebra[SingularValues](Jaclhs,output='Vt');
	nv:=HermitianTranspose(SubMatrix(TransMatrix,-1,1..-1)):
	userinfo(3,MultipleRootRefinerBreadthOne,print(`The nullvector:`,nv));
	TransMatrix:=Matrix([nv,HermitianTranspose(SubMatrix(TransMatrix,1..-2,1..-1))]);
	appz:=HermitianTranspose(TransMatrix).appy;
	nlsys:=TransPol(lsys,TransMatrix,ord);
	Ev:=[seq(ord[i]=appz[i],i=1..n)]:
	Jaclhs:=SubMatrix(eval(VectorCalculus[Jacobian](nlsys,ord),Ev),1..-1,2..-1);
	userinfo(3,MultipleRootRefinerBreadthOne,print(`The new Jacobian:`,Jaclhs));	
	dualbasis:=[1,ord[1]];
	multiplicity:=2;
	flag:=0;
	while flag<tol do
		dualsuppx:=DifOperator(dualbasis[-1],ord,1);
		dualsupp:=DifOperator(dualbasis[-1],ord,1);
		rest:=dualbasis;
		for i from 2 to n do
			for j from 2 to nops(dualbasis)-2 do
				rest[j+1]:=Part(rest[j+1],ord,i-1);
				dualsupp:=dualsupp+expand(MDifOnPol(dualsuppx,ord,ord[1]^j*ord[i])/MDifOnPol(dualbasis[j+1],ord,ord[1]^j)*rest[j+1]*ord[i]);				
			od;
		od;
		NewMatrix:=Matrix([convert(eval(DifOnPol(nlsys,ord,dualsupp),Ev),Matrix),Jaclhs]);
		flag:=LinearAlgebra[SingularValues](NewMatrix)[n];
		userinfo(2,MultipleRootRefinerBreadthOne,print(`The smallest singular value:`,flag));
		if flag>tol then
			break;
		else
			nullvector:=HermitianTranspose(SubMatrix(LinearAlgebra[SingularValues](NewMatrix,output='Vt'),-1,1..-1));
			adddual:=Transpose(nullvector).convert([dualsupp,op(ord)[2..-1]],Vector);
			dualbasis:=[op(dualbasis),expand(adddual[1])];			
			multiplicity:=multiplicity+1;
		fi;
		userinfo(2,MultipleRootRefinerBreadthOne,print(`The multiplicity:`,multiplicity));
	od;
	userinfo(1,MultipleRootRefinerBreadthOne,print(`The multiplicity:`,multiplicity));
	
	RETURN(appy+LinearAlgebra[LeastSquares](NewMatrix,-convert(eval(DifOnPol(nlsys,ord,dualbasis[-1]),Ev),Vector))[1]/multiplicity*convert(nv,Vector));

end proc: