#------------------------------------------------------------------------------------------------------------------
# author : Xiaoli Wu
# date   : Nov. 09,  2007.
#------------------------------------------------------------------------------------------------------------------

#########  product of two sets ##########################################################################################
##   Calling sequence:
##             SMultiple(A,B)
##   Input:  Two arrays A,B with finite elements
##
##   Output:  The product set of A and B
##
#########################################################################################################################
 SMultiple:=proc(A,B)
 
 
  
  local C,i,j;

  C:={};
  for i from 1 to nops(A) do
      for j from 1 to nops(B) do C:= C union {expand(A[i]*B[j])};
      end do;
  end do;
  return(C);
  end proc;




#### computing the index rho  #############################################################################################
########################################################################################################################### 
##
##    Calling sequence:
##              
##                          PrimComp(lsys,z,t,var)
##    
## Input:
##         lsys=[f1,f2,...,fm] ----    the list of  the polynomials
##         z=[z1,z2,...,zn]    ----    an isolated multiple approximate zero of  polynomial system lsys 
##         t                   ----    the tolerance 
##         var                 ----    the list of variables
## Output:
##         rho                ----     the index of the primary component
##         U                   ----     a matrix that generates the basis 
##         mu                  ----     the multiplicity
##         M                   ----     the list of multiplicative matrices     
############################################################################################################################

   PrimComp:=proc(lsys,z,t,var)

   local 
         i,
         n,
         mu,
         rho,
         sol1,sol2,sol,
         M,
         U,
         P,PP,ps;
        
   ps:=evalf(lsys);
  
   P:=[];
          
   n:=nops(var);
   for i from 1 to n do
   P:=[op(P),evalf(var[i]-z[i])];
   end do;
   
   
   PP:=SMultiple(P,P);
   sol1:=polysolverone([op(ps),op(PP)], var, [U], 2, t);  
   PP:=SMultiple(PP,P);
   sol2:=polysolverone([op(ps),op(PP)], var, [U], 2, t);
   if sol1[1]=1 and sol2[1]=1 then 
     print(`The solution is simple`);
     return([1, op(sol1)]);
   fi;


     i:=2;
     while sol1[1]<>sol2[1] do
        PP:=SMultiple(PP,P);
        sol1:=sol2;
        sol2:=polysolverone([op(ps),op(PP)], var, [U], 2, t);
        i:=i+1;
     end do;
     
   
    rho:=i;
    mu:=sol2[1];
    return([rho, op(sol2)]);
 end proc;






##### Zhengfeng Yang##########################################################################
 # The routine  computes a basis of the multivariate polynomial. 
 numbve:=proc(n,m)
   local i,S,v,N;
     S:=[];
     N:=(n+m)!/(n!*m!);
       for i from 0 to N-1 do
         v:=combinat[inttovec](i,m);
         S:=[op(S),v];
       end do;
       RETURN(S);
  end proc:
##############################################################################################




ConM:=proc(M,DegList,z)

local 
i,
v,
m,
VarL,
NewM,
DegL;


i:=0;
v:=z;
VarL:=M;
m:=LinearAlgebra:-RowDimension(M[1]);
NewM:=Matrix(m,m,shape=identity);
DegL:=DegList;


 for i from 1 to nops(DegL) do
     
     if DegL(i)<>0 then
     NewM:=NewM.(VarL[i]-v[i])^(DegL[i]);
     fi;
 end do;

return(NewM);


end proc;



NewConM:=proc(M,r,DegList,mroot)

local 
i,j,
v,
m,
VarL,
NewM,
DegL;


i:=0;
v:=mroot;
VarL:=M;
m:=LinearAlgebra:-RowDimension(M[1]);
NewM:=r;
DegL:=DegList;



 for i from 1 to nops(DegL) do
     
      
     if DegL[i]>0 then
       
       for j from 1 to DegL[i] do 
         
     NewM:=NewM.(VarL[i]-v[i]);
        
       od;
   
     fi;
 end do;

 NewM:=map(simplify, map(fnormal, NewM, Digits-1));

return(NewM);


end proc;





Mgen:=proc(M,m,z)
local

MonVec,
i,
v,n,
VarL,
MM,
Tdegree,
Nop1;

v:=z;
Tdegree:=m;
n:=LinearAlgebra:-RowDimension(M[1]);
VarL:=M;
MonVec:=numbve(Tdegree,nops(VarL));
MM:=Matrix(n,n);
i:=0;


for i from 1 to nops(MonVec) do
   MM:=MM+c[op(MonVec[i])]*ConM(VarL,MonVec[i],v);    
end do;


 return(MM);

end proc;




     

NewMgen:=proc(M,r, m,z)
local

MonVec,
i,
v,n,
VarL,
MM,
Tdegree,
Nop1;

v:=z;
Tdegree:=m;
n:=LinearAlgebra:-RowDimension(M[1]);
VarL:=M;
MonVec:=numbve(Tdegree,nops(VarL));
MM:=Transpose(Vector(n));

for i from 1 to nops(MonVec) do
   MM:=MM+c[op(MonVec[i])]*NewConM(VarL,r,MonVec[i],v);    
end do;


 return(MM);

end proc;





#### computing differential conditions  ###################################################################################
########################################################################################################################### 
##
##    Calling sequence:
##              
##                          DiffCon(mu,M,U,z)
##    
## Input:
##         U                   ----     a matrix that generates the basis 
##         mu                  ----     multiplicity
##         M                   ----     the list of multiplicative matrices 
##         z                   ----     zero of the polynomial system        
## Output:
##         L                   ----     the set of differential conditions  
############################################################################################################################

 


 NewDiffCon:=proc(rho, mu,M,U,z)
   local  
      i,j,
      zero,
      S,
      m,n,mm,
      r,s,
      MM,
      L,LL,N;
  

 m:=LinearAlgebra:-RowDimension(U);
 LL:=[];
 n:=nops(M);
 S:=[];
 mm:=mu-1;
 zero:=z;
 r:=LinearAlgebra:-Row(U,m);
 MM:=NewMgen(M, r, rho-1,zero);
 return(MM);

end proc;






#### computing differential conditions  ###################################################################################
########################################################################################################################### 
##
##    Calling sequence:
##              
##                          DiffCon(mu,M,U,z)
##    
## Input:
##         U                   ----     a matrix that generates the basis 
##         mu                  ----     multiplicity
##         M                   ----     the list of multiplicative matrices 
##         z                   ----     zero of the polynomial system        
## Output:
##         L                   ----     the set of differential conditions  
############################################################################################################################

 DiffCon:=proc(rho, mu,M,U,z)
   local  
      i,j,
      zero,
      S,
      m,n,mm,
      r,s,
      MM,Mlist,
      L,LL,N;
  

 m:=LinearAlgebra:-RowDimension(U);
 LL:=[];
 n:=nops(M);
 S:=[];
 mm:=mu-1;
 zero:=z;
 Mlist:=M;


 r:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(U),m);
 
 MM:=Mgen(Mlist,rho-1,zero);






 L:=LinearAlgebra:-Transpose(Vector(mu));
 for i from 1 to mu do
 L:=L+expand(r[i]*LinearAlgebra:-Row(MM,i));


 
 end do;

 return(L);

end proc;




################################################## Using Diffrential Conditions to generate new system ##########################
  
 
################################################## Part I : Action ##############################################################
################################################################################################################################# 
##
##    Calling sequence:
##              
##                         Action(Poly,var,MonL)
##    
## Input:
##         Poly              ----     input polynomial 
##         var               ----     list of variables
##         MonL              ----     a differrential condition that has one differrential operator
##                                    
## Output:
##         PolyOut           ----     output polynomial after the function of the differential condition  
################################################################################################################################
 
  
  
  Action:=proc(Poly,var,MonL)
    local 
       i,j,
       SubI, 
       PolyOut,
       VarList,
       Test,
       Coef,
       Monl;
  Coef:=coeffs(expand(MonL));
  Monl:=MonL;
  VarList:=var;
  PolyOut:=Poly;
  Test:=type(Monl, `*`);
  SubI:=[op(Monl)];
    if Test then
        
           Monl:=indets(Monl)[1];
           SubI:=[op(Monl)];    
                        for j from 1 to nops(SubI) do
                             
                                 for i from 1 to SubI[j] do
                                 PolyOut:=diff(PolyOut, VarList[j]);
                           
                                 end do;
                                 PolyOut:=1/(SubI[j]!)*PolyOut;
                         end do;
                         
 
                  
    else  
         
                   for j from 1 to nops(SubI) do
                          for i from 1 to SubI[j] do
                                 PolyOut:=diff(PolyOut, VarList[j]);
                                 
                           end do;
                             
                            PolyOut:=1/(SubI[j]!)*PolyOut;

                     end do;
 
 
 
 
   end if;
 
 

        
 
 PolyOut:=Coef*PolyOut;

  return(PolyOut);
 end proc:



################################################## Part II : PolyAction###########################################################
################################################################################################################################## 
##
##    Calling sequence:
##              
##                         PolyAction(Poly,var,PolyL)
##    
## Input:
##         Poly               ----     input polynomial 
##         var                ----     list of variables
##         PolyL              ----     input differrential condition that has many differential operators  
##                                   
## Output:
##         PolyOut            ----     output polynomial after the function of the differrential condition  
##################################################################################################################################
 

 PolyAction:=proc(Poly,var,PolyL)
 
 
    local 
       i,j,
       List1,MonList, DPolyL,
       PolyOut,Polyin,
       VarList,
       A,Test;
 
 
 
 i:=0;
 j:=0;
 List1:=[];
 MonList:=[];
 PolyOut:=0;
 Polyin:=Poly;
 VarList:=var;

 DPolyL:=expand(PolyL);
  if type(DPolyL, `+`)  then 
        MonList:=[op(DPolyL)];              
  else 
       
       PolyOut:=Action(Polyin, VarList, DPolyL);
 end if;
 
 for j from 1 to nops(MonList) do
 PolyOut:=PolyOut+Action(Polyin,VarList,MonList[j]);
 
 end do;
 
 return(PolyOut);
 
 end proc:


############################  construct new equations  #####################################################################
 ########################################################################################################################### 
 ##
 ##    Calling sequence:
 ##              
 ##                          NewEqu(lsys,var,L)
 ##    
 ## Input:
 ##         lsys=[f1,f2,...,fn] ----    input list of  the polynomials
 ##         L                   ----    input a vector  of the differential conditions 
 ##         var                 ----    list of variables
 ## Output:
 ##         NewSys              ----    new list of polynomials
 ############################################################################################################################
 
   NewEqu:=proc(lsys,var,DiffL)
   
      local
         i,j,k,
         FNewSys,NewSys,Lsys, neweq,
         List,VarList;
  

 
 NewSys:=evalf(lsys);
 Lsys:=evalf(lsys);
 List:=DiffL;
 VarList:=var;
 FNewSys:=[];

 
 
 for i from 1 to nops(Lsys) do
   for j from 1 to Dimension(List) do
 
       NewSys:=[op( {op(NewSys)} union {PolyAction(Lsys[i],VarList,List[j])}minus{0.,op(Lsys)} )];
      
    end do;
 end do;
for k from 1 to nops(NewSys) do 

    neweq:=fnormal(NewSys[k],Digits-1);
    if type(neweq, `*`)  then
      FNewSys:={op(FNewSys), op(2,neweq)};
    elif type(neweq, `+`) then
       FNewSys:={op(FNewSys), neweq};
    fi;
    
 od;
 
 return([op(FNewSys)]);
 
 end proc:






##########  evaluate the matrix #############################################################################
#   Calling sequence:
#             evalM(M,v,var)
#   Input:   M     ---- matrix with indeterminants entries 
#            v     ---- input value vector        
#            var   ---- list of variables
#   Output:  N     ----  the constant matrix 
#
#############################################################################################################
 
 
 
 
   evalM:=proc(M,v,var)
 
 
 
   local i,n,N;
 
   N:=M;
   n:=nops(var);
 
   for i from 1 to n do
      N:=eval(N,var[i]=v[i]);
   end do;
 
   return(N);
   end proc:







##########  Gauss-Newton method to refine a root of a polynomial system ##########################################
##
##      Calling sequence:
##             GN(lsys,v,t,var)
##   Input:    lsys       ---- a polynomail system
##             v          ---- an initial vector value (close to the real zero)             
##             t          ---- a tolerance
##             var        ---- list of variables
##   Output:   w          ---- the refined zero of system
##             k          ---- the number of iterations
###################################################################################################################
   GN:=proc(lsys,v,t,var)
   
   
      local 
          k,u,w,r,rr,i,
          J,M,MM, ps; 
          
     ps:=lsys;
  
    k:=1;
    u:=Vector(v);
    J:=VectorCalculus:-Jacobian(ps,var);
    M:=evalM(J,u,var);
  
    if LinearAlgebra:-SingularValues(LinearAlgebra:-HermitianTranspose(M))[-1] >t  then
      r:=evalM(Vector(ps),u,var);
     w:=u-LinearAlgebra:-LeastSquares(M,r, optimize, conjugate=true);
   
     while LinearAlgebra:-Norm(w-u,2)>=10.^(-3*Digits+1) and k<20 do 
         
          u:=w; 
          MM:=evalM(J,u,var); 
          rr:=evalM(Vector(ps),u,var);   
          w:=u-LinearAlgebra:-LeastSquares(MM, rr, optimize, conjugate);
          k:=k+1;
    end do;
    return(w,k);
    else
    return("The Jocobian matrix is singular,  Gauss-Newton iteration does not converge.");
    end if;
  
   end proc: