# Author: Lihong Zhi
# Date: 2002--2008

#Order the monomials according to the inttovect(combinat)

selectcolumn_xi:=proc(ordi, tdeg, nvars)
  local  tterms, L, j;
  
   tterms:=combinat[binomial](nvars+tdeg+1, nvars);
   L:=[];
   for j from   ordi to tterms do 
     if  combinat[inttovec](tterms-j,nvars)[nvars-ordi+1]<> 0 then
       L:=[op(L), j];
     fi;
   od;
   RETURN(L);
  end proc:
    
   
      

#Set up the order of monomials  w.r.t Greg's monomial order

Gordmons:=proc()
   local n,m,q,N,Nq,dyd, indeps,deps,sys, i, psys,dOrder, DRing, vars,
         MaxNumDiffs;
   #   Make sure the matrix has its cols order according to class (not done yet)
   sys:=args[1]:  # the differential system [p1,p2,...]
   indeps:=args[2]:  # the variables [x,y,z..]
   deps  :=args[3]:   # the dependent variable U
  
     i:=args[4];   # Inv. Test passed at i=
   
    q:= PDEtools[difforder](sys);
   
       
       psys:=Prolong(sys,indeps,i-q);
       dOrder:=PDEtools[difforder](psys):
       DRing:= diffalg[differential_ring](derivations=indeps, ranking=deps, notation=diff);

       psys:= map(z -> z = 0, reorder(map(diffalg[denote], psys, jet, DRing))): 
       vars:= DEtools[checkrank]([seq(deps[j](op(indeps)), j = 1 .. nops(deps))],
                                      deps, indep=indeps,degree=dOrder);
       vars:= map(diffalg[denote], vars, jet, DRing):
   
    RETURN(vars);
    end proc:
    
    
    
    
 #Set up the multiplication table w.r.t the monomial order defined by Greg
   Gselectcolumn_xi:=proc()
         local n,m,q,N,Nq,dyd, indeps,deps,sys, i,j, psys,dOrder, DRing, vars,
         mons,vvars,xi, dlsys, L, LL, temp;
     xi:=args[1]:     # the multiplication table w.r.t xi
     dlsys:=args[2]:    # the system of differential equations
     indeps:=args[3]:  # the variables [x,y...]
     deps  :=args[4]:  # the dependent variable [U]
      i     :=args[5]:    # Inv. Test passed at i=
      j     :=args[6]:    #  Inv. Test passed at j=
     
      
      
     L:=Gordmons(dlsys, indeps, deps, i);
    
     mons:= binomial(nops(indeps)+j-1, nops(indeps));               # the total number of mons degree <= j-1
   
     L:=L[(nops(L)-mons+1)..nops(L)];
     
     LL:=[];
     for i from 1 to nops(L) do 
        temp:=[op(L[i])];
        if member(xi,temp)  then
          LL:=[op(LL), i];
        fi;
      od;
     RETURN(LL);
     end proc:
         
         
         
         
         
         
 #Solve the polynomial system by  SNED and  eigenvalue method 

polysolver:=proc()
   local sI, ps,vars, depvars, MaxNumDiffs,
       RankTol, dlsys, Sdlsys, Condv, ZZ,TT,L, Mx,LL,Lx,Tx, a, k,E1, E2,B,BB,randnumb, M,MM, e1, ev1, ev, sol, EB1, RM,RN,numsol,i, j, SE21, SE22, SE23, SE24,U1,S1,V1;
    
    
     ps:=args[1]:  # the polynomial system [p1,p2...]
     vars:=args[2]:  # the variables [x,y, z...]
     depvars:=args[3]: # the dependent variables [U...]
     if nargs=3 then 
       MaxNumDiffs:=5:	        	#DEFAULT
       RankTol:= 0.5e-3;		#DEFAULT
     elif nargs=4 then
       MaxNumDiffs:=args[4]:            #USER DEFINED
       RankTol:= 0.5e-3;		#DEFAULT
     elif nargs=5 then
       MaxNumDiffs:=args[4]:	        #USER DEFINED
       RankTol:= args[5];		#USER DEFINED
     end if;
   #PolynomialToPDE has problem to transfer polynomial with complex coefficients, here we introduce a new parameter sI
     ps:=subs(I=sI,ps);
     dlsys:=subs(sI=I, PolynomialToPDE(ps, vars, depvars));
     Sdlsys:=Sned(dlsys, vars, depvars, MaxNumDiffs, RankTol);
    # lprint(`DimMtx=`, Sdlsys[DimMtx]);
    # lprint(`RankSymbolMtx=`, Sdlsys[RankSymbolMtx]);     
   # check if the system is involutive and the number of solutions
     if Sdlsys[InvolutiveTest]=true  then
        i:=Sdlsys[Invrow];
        j:=Sdlsys[Invcol];
      
           
        numsol:=Sdlsys[DimMtx][i+1,j+1];
        print(` number of solutions`, numsol);
           # choose the smallest i such that DimMtx[i,j]=DimMtx[i, j-1]=numsol
          while Sdlsys[DimMtx][i,j]=numsol  and Sdlsys[DimMtx][i,j-1]=numsol do 
               i:=i-1;
           od:
         # form the multiplication matrix
         
        #  choose E1, E2 of  row dimension as numsol.
           E1:= HermitianTranspose(Sdlsys[ProjSpanSets][i, j-2]);
           E2:= HermitianTranspose(Sdlsys[ProjSpanSets][i, j-1]);
           
            
         # There is only one solution
           if numsol=1  then
                print(`there is only one solution`);
                 sol:=[[1,{ 'vars[nops(vars)-t+1]=E2[RowDimension(E2)-t,1]/E2[RowDimension(E2),1]' $ 't'=1..nops(vars)}]]; 
            
               return(sol); 
            fi;
                
            U1,S1,V1:=LinearAlgebra[SingularValues](E1, output=['U', 'S','Vt']):
            
            U1:=SubMatrix(U1,1..RowDimension(U1), 1..numsol):
            V1:=SubMatrix(V1, 1..numsol, 1..ColumnDimension(V1)):
            S1:=DiagonalMatrix(['S1[k]' $'k'=1..numsol]):
            
             randnumb:=randnumber(nops(vars));
           
            
            B[1]:=SubMatrix(E2, Gselectcolumn_xi(vars[1],dlsys, vars, depvars,i,j), 1..ColumnDimension(E2));
            BB:=B[1]:
            for k from 2 to nops(vars) do 
            B[k]:=SubMatrix(E2, Gselectcolumn_xi(vars[k],dlsys, vars, depvars,i,j), 1..ColumnDimension(E2));
            BB:=BB+randnumb[k].B[k]:
            od;
           
           
            M:= HermitianTranspose(U1).BB. HermitianTranspose(V1).MatrixInverse(S1):
           
            
            Condv:=EigenConditionNumbers(M):
           
            
            if min(op(convert(Condv[1],list)), op(convert(Condv[2],list)))>10^(-4) or numsol=1 then #good condition, without multiple roots
               #lprint(`good case`,M);
               e1,ev1:=Eigenvectors(M):
               ev:=U1.ev1;
               sol:=[];
               for k from 1 to numsol do 
               sol:=[op(sol),[1,{ 'vars[nops(vars)-t+1]=ev[RowDimension(ev)-t,k]/ev[RowDimension(ev),k]' $ 't'=1..nops(vars)}]]; od;
               return(sol);
            else
                print(`the system has multiple roots`):
     
                M:=Matrix( M, datatype=complex(sfloat)):
                
                ZZ:= ReorderSchur(M,10^(-3)): # slight bigger than the ranktol
               
                TT:= HermitianTranspose(ZZ).M.ZZ;
               
                L:=groupSchur(TT, 10^(-3)):    #group the eigenvalues
                
                LL:=[]:
                for j from 1 to nops(vars) do
                      Mx:= HermitianTranspose(U1).B[j]. HermitianTranspose(V1).MatrixInverse(S1): # the multiplication table w.r.t xj
                     
                      Tx:= HermitianTranspose(ZZ).Mx.ZZ:
                      userinfo(3,polysolver,print(`Check whether the transformed matrices are reorder upper block matrix`, Tx));
                      Lx:=averSchur(Tx, L):
                      LL:=[op(LL),Lx]:
                od: 
                sol:=[]:
                for k from 1 to nops(L) do 
                    sol:=[op(sol),[L[k],{ 'vars[t]=LL[t][k][2]' $ 't'=1..nops(vars)}]]; 
                od:
             return(sol);
        
                
                
             fi:  
               
               
            
            
       fi:
         return [];
end proc;
          
 
 
#Group the diagonal elements of a matrix, output the averages
groupSchur:=proc(T, tol)
local dim, L, i, a,k,j;

   dim:=RowDimension(T);  
  
   L:=[]:
   i:=1:
   while i <= dim do 
         a:=T[i,i]:
         k:=1:
         if i=dim then
           L:=[op(L),1]: break:
         else
           j:=i+1:
          
           while abs(a-T[j,j])<=tol and j <=dim  do
             if j=dim then
               k:=k+1:
               RETURN([op(L),k]):
             fi:
             k:=k+1:
           
             j:=j+1:
           od:
           
            L:=[op(L),k]:
           i:=j:
          fi:
   od:
   RETURN(L):
  end proc:     
         
     
 # Compute the average of the elements w.r.t group L
  averSchur:=proc(T, L)
     local dim, LL, i, k, j,suma,a;
     dim:=RowDimension(T); 
     LL:=[]:
     i:=1;
     k:=1;
     while i<=dim do
        a:=T[i,i]:
        if i=dim  then
           LL:=[op(LL),[L[nops(L)],a]]: break
         elif L[k]=1 then
            LL:=[op(LL),[1,a]]:
            k:=k+1:
            i:=i+1:
         else
             suma:=a:
             for j from 1 to L[k]-1 do
                suma:=suma+T[i+j,i+j]:
             od:
             LL:=[op(LL), [L[k],suma/L[k]]]:
             i:=i+L[k]:
             k:=k+1:
          fi:
      od:
     RETURN(LL):
  end proc:     
 
 
 
 #Choose a random list [a1, a2...] a1+a2+...=1
 randnumber:=proc(n)
  local i,a,die2,suma,die1;
    a:=[];
    die2:=rand(-100..100);
    for i from 1 to n do 
     
      a:=[op(a), die2()/100];
    od;
      suma:=sum('abs(a[i])','i'=1..nops(a));
      a:=['a[i]/suma' $ 'i'=1..nops(a)];
      RETURN(a);
end proc:
      
      
      
 
     
 
# Compute c,s, for exchanging  tii, tjj, j=i+1
Schurgiven:=proc(A, i,tol)
  local IR, R, ei,ej,c,s,a,e,b,r;
      
      IR:=IdentityMatrix(RowDimension(A)); # the identity matrix of size dim
      a:=A[i,i]:
      e:=A[i+1,i+1]:
      b:=A[i,i+1]:
      if abs(a-e)<tol then
        R:=IR:
      elif abs(b)< tol then
         R:=IR:
         R[i,i]:=0;
         R[i,i+1]:=1:
         R[i+1,i]:=1:
         R[i+1,i+1]:=0:
      else  
        r:=sqrt(abs(a-e)^2+abs(b)^2):
        c:=abs(b)/r:
        s:=(a-e)/r*(conjugate(b))/abs(b):
       
        ei:=Column(IR,i);
        ej:=Column(IR,i+1); 
        R:=IR+(c-1)*ei.Transpose(ei)-conjugate(s)*ei.Transpose(ej)+s*ej.Transpose(ei)+(conjugate(c)-1)*ej.Transpose(ej):
      fi:
     
  
     RETURN(R);
  end proc:
     

# Recursive rotation to move the j-th row to the i-th row 
  RotationSchur:=proc(A, i, j, tol)
      local R, T, k,R1;
       R:=IdentityMatrix(RowDimension(A)); # the identity matrix of size dim
       T:=A:
      
       if j <= i then 
          RETURN([T,R]):
       else
          for k from j-1 by -1  to i do 
             R1:=Schurgiven(T,k,tol):
             R:=R1.R:
             T:= R1.T.HermitianTranspose(R1);
          od:
         
          for k from i+1 to j-1 do
                R1:=Schurgiven(T,k,tol):
                R:=R1.R:
             T:= R1.T.HermitianTranspose(R1);
          od:
           
              
          
         RETURN([T,R]):
       fi
    end proc:
       


 
 # Reorder the Schur decomposition
 ReorderSchur:=proc(A, tol)
   local TT, ZZ, L,i,dim,j,a,R,Bn,k;
     
       
     (TT, ZZ) := SchurForm(A, output=['T', 'Z'] ):
     
     dim:=RowDimension(TT):
     i:=1:
     Bn:=[]: # the size of blocks
     while i <dim do 
         a:=TT[i,i]:
         k:=1:
         j:=i+1:
         while abs(a-TT[j,j])<=tol and i <dim-1  do
           i:=i+1:
           k:=k+1:
           j:=i+1:
         od:
        
         if   i <dim-1 then
            a:=TT[i,i]:
           for j from i+2 to dim do 
               if abs(a-TT[j,j])<= tol then
               L:= RotationSchur(TT,i+1, j, 10^(-2)*tol):
              
               TT:=L[1]:
               R:=L[2]:
               ZZ:=ZZ.HermitianTranspose(R):
               i:=i+1:
               k:=k+1:
              fi:
            od:
         fi:
         
        
         if i=dim-1 then 
         Bn:=[op(Bn),k,1]:
         else
         Bn:=[op(Bn),k]:
        
         fi:
          i:=i+1:
         
       od:
        
    RETURN(ZZ):
    
    end proc:
       
          
                            
        
#Solve the polynomial system  with only one multiple solution by  SNED and  eigenvalue method

polysolverone:=proc()
   local sI, ps,vars, depvars, MaxNumDiffs,  MList,#######################
       RankTol, dlsys, Sdlsys, Condv, ZZ,TT,L, Mx,LL,Lx,Tx, a, k,E1, E2,B,BB,randnumb, M,MM, e1, ev1, ev, sol, EB1, RM,RN,numsol,i, j, SE21, SE22, SE23, SE24,U1,S1,V1;
    
    MList:=[];
     ps:=args[1]:  # the polynomial system [p1,p2...]
     vars:=args[2]:  # the variables [x,y, z...]
     depvars:=args[3]: # the dependent variables [U...]
     if nargs=3 then 
       MaxNumDiffs:=3:	        	#DEFAULT
       RankTol:= 0.5e-3;		#DEFAULT
     elif nargs=4 then
       MaxNumDiffs:=args[4]:            #USER DEFINED
       RankTol:= 0.5e-3;		#DEFAULT
     elif nargs=5 then
       MaxNumDiffs:=args[4]:	        #USER DEFINED
       RankTol:= args[5];		#USER DEFINED
     end if;
   #PolynomialToPDE has problem to transfer polynomial with complex coefficients, here we introduce a new parameter sI
     ps:=subs(I=sI,ps);
     dlsys:=subs(sI=I, PolynomialToPDE(ps, vars, depvars));
    
     Sdlsys:=Sned(dlsys, vars, depvars, MaxNumDiffs, RankTol);
        
   # check if the system is involutive and the number of solutions
     if Sdlsys[InvolutiveTest]=true  then
        i:=Sdlsys[Invrow];
        j:=Sdlsys[Invcol];
       
        numsol:=Sdlsys[DimMtx][i+1,j+1];
      
           # choose the smallest i such that DimMtx[i,j]=DimMtx[i, j-1]=numsol
          while Sdlsys[DimMtx][i,j]=numsol  and Sdlsys[DimMtx][i,j-1]=numsol do 
               i:=i-1;
           od:
         # form the multiplication matrix
          
        #  choose E1, E2 of  row dimension as numsol.
           E1:= HermitianTranspose(Sdlsys[ProjSpanSets][i, j-2]);
           E2:= HermitianTranspose(Sdlsys[ProjSpanSets][i, j-1]);
          
            
         # There is only one solution
           if numsol=1  then
                print(`there is only one solution`);
                 sol:=[1,{ 'vars[nops(vars)-t+1]=E2[RowDimension(E2)-t,1]/E2[RowDimension(E2),1]' $ 't'=1..nops(vars)}]; 
          
               return(sol); 
            fi;
                
            U1,S1,V1:=LinearAlgebra[SingularValues](E1, output=['U', 'S','Vt']):
            
            U1:=SubMatrix(U1,1..RowDimension(U1), 1..numsol):
            V1:=SubMatrix(V1, 1..numsol, 1..ColumnDimension(V1)):
            S1:=DiagonalMatrix(['S1[k]' $'k'=1..numsol]):
            sol:={};
            for k from 1 to nops(vars) do 
            B:=SubMatrix(E2, Gselectcolumn_xi(vars[k],dlsys, vars, depvars,i,j), 1..ColumnDimension(E2));
            M:= HermitianTranspose(U1).B. HermitianTranspose(V1).MatrixInverse(S1):
            MList:=[op(MList),M];
            ev:=Eigenvalues(M);
            sol:={ vars[k]=sum('ev[i]', 'i'=1..numsol)/numsol, op(sol)};
            od;
           # Mlist is the list of multiplication matrices, can be returned  
           # U1 is the basis generator, can be returned
            return([numsol, sol]);
       fi; return [];
              
end proc;
          
 
 

         
 

 
 
################################################################################################## 
## The following procedures are written   by Zhengfeng Yang 
# This routine computes the matrix which represents the multiplication of 
# f by a polynomial of total degree k.  This is called the Cauchy matrix
# by some.
# This routine generalize John's definition for bivariate Cauchy matrix
mvcauchyo := proc(f, ord,deg)
        local  N,td,i, unknowns, g, eqns, pp,V
             ,a ,n,A; # Protect the unknowns
       td := degree(f,{op(ord)});
       n:=nops(ord);
      # N:=(n+deg)!/(n!*deg!);
        N:=binomial(n+deg,deg);
       g:=0;
       V:=Vector(N);
       for i from 1 to N do
         V[i]:=a[i];
       end do;
         pp:=convert(V,list);
         g:=veccoe(V,ord);
         eqns := expand(f*g);
         eqns := convert(convec(eqns,ord),list);
         eqns := map(x->x=0, eqns);
         A:=LinearAlgebra[GenerateMatrix](eqns, pp)[1];
         RETURN(A);
end proc;




# This routine computes the matrix which represents the multiplication of 
# f by a polynomial of total degree k.  This is called the Cauchy matrix
# by some.
# This routine generalize John's definition for bivariate Cauchy matrix
mvcauchyok := proc(f, ord,deg,k)
        local  N,td,i, unknowns, g, eqns, pp,V
             ,a ,n,A; # Protect the unknowns
       td := degree(f,{op(ord)});
       n:=nops(ord);
        N:=binomial(n+deg,deg);
       g:=0;
       V:=Vector(N);
       for i from 1 to N do
         V[i]:=a[i];
       end do;
         pp:=convert(V,list);
         g:=veccoe(V,ord);
         eqns := expand(f*g);
         eqns := convert(conveck(eqns,ord,k),list);
         eqns := map(x->x=0, eqns);
         A:=LinearAlgebra[GenerateMatrix](eqns, pp)[1];
         RETURN(A);
end proc;






# This routine converts a vector to a multivariate polynomial.
veccoe:=proc(v,ord)
  local i,n,N,f0;
    n:=nops(ord);
    N:=LinearAlgebra[Dimension](v);
    f0:=0;
    for i from 1 to N do
      f0:=f0+v[i]*vecmo(combinat[inttovec](N-i,n),ord);
    end do;
 end proc;   
 
# This routine converts a vector to a  monomial.
vecmo:=proc(v,ord)
 local i,n,s;
   n:=nops(ord);
   s:=1;
   for i from 1 to n do
     s:=s*ord[i]^(v[i]);
   end do;
   RETURN(s);
end proc:

# The routine  computes the coefficients w.r.t. given degrees of variables.
 convs:=proc(f,ord,V)
   local i,n,f0;
     n:=nops(ord);
     f0:=f;
     for i from 1 to n do;
       f0:=coeff(f0,ord[i],V[i]);
     end  do;
       RETURN(f0);
 end proc;
 
 # The routine  computes a basis of the multivariate polynomial. 
numbve:=proc(n,m)
  local i,S,v,N;
    S:=[];
     N:=binomial(n+m,m);
      for i from 0 to N-1 do
        v:=combinat[inttovec](i,m);
        S:=[op(S),v];
      end do;
      RETURN(S);
 end proc:

 # The routine converts the coefficients of a  multivariate polynomial to   a vector.
 convec:=proc(f,ord)
   local i,m,n,vba,ve,N;
   m:=degree(f,{op(ord)});
   n:=nops(ord);
    N:=binomial(n+m,m);
   vba:=numbve(m,n);
   ve:=Vector(N);
   for i from 1 to N   do
     ve[N-i+1]:=convs(f,ord,vba[i]);
   end do;
   RETURN(ve);
end proc;  


 # The routine converts the coefficients of a  multivariate polynomial to   a vector.
 conveck:=proc(f,ord,k)
   local i,m,n,vba,ve,N;
   m:=k;
   n:=nops(ord);
    N:=binomial(n+m,m);
   vba:=numbve(m,n);
   ve:=Vector(N);
   for i from 1 to N   do
     ve[N-i+1]:=convs(f,ord,vba[i]);
   end do;
   RETURN(ve);
end proc;  



# The routin  formulates the coefficient vector w.r.t polynomials up to given degree 

 coeffveck:=proc(f, ord, k)
    local i,m,n,vba,ve,N;
     
       m:=degree(f,{op(ord)});
       n:=nops(ord);
       if m=k then 
       RETURN(convec(f, ord));
       else
          m:=k;
           N:=binomial(n+m,m);
         vba:=numbve(m,n);
        ve:=Vector(N);
         for i from 1 to N   do
        ve[N-i+1]:=convs(f,ord,vba[i]);
       end do;
       RETURN(ve);
       fi;
      end proc:

 
 
 
 
#  The routin formulates the coefficient vector w.r.t polynomials up to given degree 

 coeffvector:=proc(f, ord, k)
    local i,m,n,vba,ve,N,Lc,Lt,s;
     
      
       n:=nops(ord);
         N:=binomial(n+k,k);
        
         ve:=Vector(N);    
         Lc:=[coeffs(expand(f),ord, 'Lt')];
         
         Lt:=[Lt];
         for i from 1 to nops(Lc) do 
          
            if degree(Lt[i],{op(ord)})<=k then
              s:=combinat[vectoint](dvector(Lt[i],ord));
           
              ve[N-s]:=Lc[i];
           fi;
         od;
      RETURN(ve);
      end proc:
           
 
 #  The routine computes the degree vector of a monomial
       
 dvector:=proc(t,ord)
 local i, L;
    L:=[];
    for i from 1 to nops(ord) do 
    L:=[op(L),degree(t,ord[i])];
   
    od;
   RETURN(L);
   end proc:
    
 

#  The routine computes all the monomials up to a given  degree m

mons:=proc(ord, m)
  local  n, v, i, mon;     
      n:=nops(ord);
      v:=numbve( m, n);
      mon:=[];
      for i from 1 to nops(v) do 
        mon:=[op(mon), vecmo(v[i],ord)];
      od;
      RETURN(mon);
  end proc:
         
    
    
#  The routine computes  the monomials of a given  degree m

monm:=proc(ord, m)
  local  n, v, i, mon, N1, N2;     
      n:=nops(ord);
      v:=numbve( m, n);
      mon:=[];
      N1:=binomial(m+n,m);
      N2:=binomial(m+n-1, m-1);
      for i from N2+1 to N1 do 
        
        mon:=[op(mon), vecmo(v[i],ord)];
      od;
      RETURN(mon);
  end proc:


 
 #  The routine prepares a polynomial system with given degree, i.e., prolonging the polynomials to the given degree
 #  The polynomials are seperated according to the degrees k
 
 prolongpolyk:=proc(lsys, ord, k)
    local ps, i, j,d,f, mon, ps1; 
      ps:={};
      ps1:={};
      for i from 1 to nops(lsys) do 
            f:=expand(lsys[i]);
            d:=degree(f,{op(ord)});
           
         if d<k then 
            mon:=mons(ord, k-d);
            for j from 1 to nops(mon) do 
              if degree(mon[j],{op(ord)})<k-d then
              ps1:={op(ps1), expand(mon[j]*f)};
              else
              ps:={op(ps), expand(mon[j]*f)};
              fi;
            od;
          else
             ps:={op(ps), f};
          fi;

       od;
           RETURN([[op(ps1)], [op(ps)]]);
       end proc:

 
 
 #  The routine prepares a polynomial system with given degree, i.e., prolonging the polynomials t times
 
 
 prolongpolyt:=proc(lsys, ord, t)
    local ps, i, j,d,f, mon, ps1; 
      ps:=[];
      ps1:=[];
      for i from 1 to nops(lsys) do 
            f:=expand(lsys[i]);
            mon:=mons(ord, t);
            for j from 1 to nops(mon) do
              ps:=[op(ps), expand(mon[j]*f)];
            od;
       od;
           RETURN(ps);
       end proc:
 



 
 #  The routine  prolongs the polynomials by multiplying all variables in ord.
 

 prolongpoly1:=proc(lsys, ord)
     local ps, i, j,d,f, mon, ps1; 
      ps:=[];
      ps1:=[];
      for i from 1 to nops(lsys) do 
            f:=expand(lsys[i]);
            for j from 1 to nops(ord) do
              ps:=[op(ps), expand(ord[j]*f)];
            od;
       od;
           RETURN(ps);
       end proc:
   






#   The routine computes the matrix Mk0, if the polynomial in lsys has degree less than k then, we prolong it to have degree k
#  if the polynomial has degree larger than k, we truncate it to be of degree k. 

 startmatrix0:=proc(lsys, ord, k)
    local A,ps,i; 
      ps:=prolongpolyk(lsys, ord, k);
      ps:=[op(ps[1]), op(ps[2])];
      A:=Matrix([coeffvector(ps[1], ord, k)]);
      for i from 2 to nops(ps) do 
      A:= Matrix([A,coeffvector(ps[i], ord, k)]);
      od;
      RETURN(LinearAlgebra[Transpose](A));
     end proc:




 
 #  Compute the maximal degree of a polynomial system
 
 maxdeg:=proc(lsys, ord)
        local d,d1,i;
          d:=degree(lsys[1],{op(ord)});
          for i from 2 to nops(lsys) do 
            d1:=degree(lsys[i],{op(ord)});
            if d1 >d then
            d:=d1;
            fi;
          od;
          RETURN(d);
          end proc:



 
 #  Compute the minimal degree of a polynomial system
 
mindeg:=proc(lsys, ord)
        local d,d1,i;
          d:=degree(lsys[1],{op(ord)});
          for i from 2 to nops(lsys) do 
            d1:=degree(lsys[i],{op(ord)});
            if d1 <d then
            d:=d1;
            fi;
          od;
          RETURN(d);
          end proc:





# Prolong the matrix from d to d+1, suppose d >k

prolongM:=proc(Ma, ord,mon, k)
   local t, n, d0,d1,i,d, ps, ts, qs,A,j,A1, A2, N,M,md, mind, maxd, mon1, mon2,cd,AA,s;
   
                  
                  mon1:=mon;
                    A:=Ma;
                  n:=nops(ord);
                         
                     M:=[];
                     for i from 1 to n do 
                       M:=[op(M), mvcauchyok(ord[i],ord, k,k)];
                     od;             
                    
                 
                        mon2:=[];
                                           
                        A1:=Matrix([M[1].A]);

                        if nops(mon1)=0 then
                           mon2:=ord;
                           for j from 2 to n do                       
                             A1:=Matrix([A1, M[j].Column(A,1)]);
                             
                          od;
                        RETURN(A1,mon2);
                       fi;



                        for i from 1 to nops(mon1) do 
                          mon2:=[op(mon2), ord[1]*mon1[i]];
                       od;
                       for j from 2 to n do 
                        for i from 1 to nops(mon1) do 
                          if not {ord[j]*mon1[i]} subset {op(mon2)} then 
                             
                             A1:=Matrix([A1, M[j].Column(A,i)]);
                             mon2:=[op(mon2), ord[j]*mon1[i]];
                             
                          fi;
                        od
                       od;
                     
                  RETURN(A1,mon2);
                  
                
                end proc:

   
   
   
   
   
#Prolong the polynomial to be of degree k, if k>deg(f), then output the truncated coefficient matrix.

prolongp:=proc(f,ord, k)
   local t, n, d0,d1,i,d, ps, ts, qs,A,j,A1, A2, N,M,md, mind, maxd,mon, mon1, mon2,cd,AA,s;
                  mon1:=ord;
                  mon2:=[];
                  mind:=degree(f, {op(ord)});     
                  A:=Matrix([coeffvector(f, ord, k)]);
                  A2:=A;
                  n:=nops(ord);
                  if mind <k  then                   
                     M:=[];
                     for i from 1 to n do 
                       M:=[op(M), mvcauchyok(ord[i],ord, k,k)];
                     od;             
                 
                    A1:=Matrix([M[1].A]);                 
                    for j from 2 to n do 
                      A1:=Matrix([A1, M[j].A]);
                    od;
                    d:=mind+1;
                    mon2:=mon1;
                    A:=A1;
                    A2:=Matrix([A2,A1]);
                    
                     while d <k do 
                        mon2:=[];
                       
                       d:=d+1;
                       A1:=Matrix([M[1].A]);
                      
                       for i from 1 to nops(mon1) do 
                         mon2:=[op(mon2), ord[1]*mon1[i]];
                       od;
                       for j from 2 to n do 
                        for i from 1 to nops(mon1) do 
                          if not {ord[j]*mon1[i]} subset {op(mon2)} then 
                             
                             A1:=Matrix([A1, M[j].Column(A,i)]);
                             mon2:=[op(mon2), ord[j]*mon1[i]];
                          
                          fi;
                        od
                       od;
                       mon1:=mon2;
                       A:=A1;
                     
                      
                      A2:=Matrix([A2, A1]);
                     
                  od;
                  RETURN([A2,A],mon2);
                  
                else
                  RETURN([A,A],mon2);
                fi;
                end proc:


   
  




########################################################################
#   Compute the dimension of the prolonged system and check involutivity w.r.t P^k
#
involutivetest:=proc(lsys, ord, k, tol)
    local t, n, d0,d1,i,d, ps, ts, qs,A,j,A1, A2, N,M,md, mind, maxd,mon, mon1, mon2,cd,AA,s,Ma,nt,L2,L;
       t:=1;
     
        lprint(`k=`,k);
        if nargs=5 then
          md:=args[5];
        else
          md:=1000;
        fi;
          n:=nops(ord);
          ts:=lsys;     
          nt:=nops(ts);
        
          A:=prolongp(ts[1],ord, k-1);
          A1:=A[1][1]; # the prolong matrix
          A2:=A[1][2]; #the top matrix
          L2:=[A2];
          mon:=[A[2]];
          
           for s from 2 to nt do 
            A:=prolongp(ts[s],ord, k-1);
            A1:=Matrix([A1,A[1][1]]);
            L2:=[op(L2),A[1][2]];
            mon:=[op(mon), A[2]];
           od;
           
            d0:=dimtol(A1,tol);
            lprint(`check the dimension d0`, d0);
            L:=[];
            mon1:=[];
           
            for i from 1 to nt do  #prolong once           
              Ma:=prolongM(L2[i], ord, mon[i], k-1);
              A1:=Matrix([A1,Ma[1]]);
              L:=[op(L), Ma[1]]; #the top matrix
              mon1:=[op(mon1),Ma[2]];
            od;
            L2:=L;
            mon:=mon1; #the new monomial set
           d1:=dimtol(A1,tol);
                  lprint(`check the dimension d1`, d1);
             while d1<d0  and d1>md  do 
                  lprint(` the system is not involutive yet, prolong the system`, t);
                  t:=t+1; 
                   d0:=d1;
                    L:=[];
                    mon1:=[];
                   for i from 1 to nt do  #prolong once           
                     Ma:=prolongM(L2[i], ord, mon[i], k-1);
                     A1:=Matrix([A1,Ma[1]]);
                     L:=[op(L), Ma[1]]; #the top matrix
                     mon1:=[op(mon1), Ma[2]];
                   od;
                   L2:=L;
                   mon:=mon1; #the new monomial set
                   d1:=dimtol(A1,tol);
                  lprint(`check the dimension d1`, d1);
                od;
               
             

    
         lprint(`the system is involutive after t prolongation`, d1, t);
           RETURN(d1, t, A1);
       
          
       end proc:

           










 ##############################################################
 #   Compute the rank deficiency by SVD
 #
dimtol:=proc(M, tol)
   local L, n, s0, t, time1, time2;
     lprint(`the dimension of the matrix`, LinearAlgebra[Dimension](M));
    
     L:=LinearAlgebra[SingularValues](evalf(M));
     
     n:=LinearAlgebra[Dimension](L); 
     if L[1]<tol then 
      RETURN(n);
     fi;
     s0:=L[n];
     t:=1;
     while s0<tol and t<n do 
       s0:=L[n-t];
       t:=t+1;
       
     od;
     if t<=n then 
     RETURN(t-1);
     else 
     RETURN(n);
     fi;
     end proc:

  
 #############################################################
 #   Compute the rank deficiency by eigenvalues
 #
dimtoleig:=proc(M, tol)

   local L, M1,   t, i;
     lprint(`the dimension of the matrix`, LinearAlgebra[Dimension](M));
   
      M1:=LinearAlgebra[HermitianTranspose](evalf(M)).evalf(M);
      L:=LinearAlgebra[Eigenvalues](M1);
    
      t:=0;
      for i from 1 to LinearAlgebra[Dimension](L) do 
        if evalf(sqrt(abs(L[i]))) < tol then 
         t:=t+1;
        fi;
      od;


     RETURN(t);
    
     end proc:

  ##################################################
  #  Compute the index and multiplicity of the primary component crosponding to the zero at origin. 
  #
  indexmulti:=proc(lsys, ord, tol)
     local k, L2, L1, d0,d1;
        k:=2;
        L1:=involutivetest(lsys,ord, k,tol);
        d0:=L1[1];
        L2:=involutivetest(lsys, ord, k+1, tol,d0);
        d1:=L2[1];
        if d1<d0 then 
          print(`the tolerance is too small, try big tolerance`);
          RETURN([]);
        fi;
        if d1=d0 then 
          
           if k<=d1 then 
            print(`the index and multiplicity are`, k, d1); 
           RETURN(k,d1,L1[3],L2[3]);
           else 
            lprint(`the index should be less or equal to multiplicity,`, d1, d1);
           RETURN(d1,d1,L1[3],L2[3]);
           fi;
        fi;
         k:=k+1;
        while d1>d0 do 
           k:=k+1;
           d0:=d1;
           L1:=L2;
           L2:=involutivetest(lsys, ord, k, tol,d0);
           d1:=L2[1];
           lprint(`the order of the primary ideal and dimension`, k, d1);
         od;
         if k<d1+1 then 
         lprint(`the index and multiplicity are `, k-1, d1);
         RETURN(k-1,d1,L1[3],L2[3]);
         else 
         lprint(`the index should be less or equal to multiplicity,`, d1, d1);
         RETURN(d1,d1,L1[3],L2[3]);
         fi;
        end proc:
          

 ###############################################################
 #  Suppose the mroot=(0), we compute the differential conditions
 #  A is returned as the last prolonged matrix, 
 #  d is the index, 
 #  m is the multiplicity

Diffcon:=proc(A,ord, d,m)
  local U1, S1, V1,V,L, L1, MonVec, j, i, N;


  U1, S1, V1:=SingularValues(Transpose(A), output=['U', 'S', 'Vt']);
  V:=HermitianTranspose(SubMatrix(V1, -m..-1,1..-1));
  
 
  L:=[];
  MonVec:=numbve(d-1,nops(ord));
  N:=nops(MonVec);
 
  for j from 1 to m do 
     L1:=0;
    for i from 1 to N do 
      if abs(Column(V,j)[i])>10^(-Digits+6) then 
      L1:=L1+c[op(MonVec[N-i+1])]*Column(V,j)[i];    
      fi;
    end do;
     L:=[op(L),L1];
  od;
  lprint(`we get the differential condition`, nops(L), L);
  RETURN(Vector(L));

end proc;



##############################################################################################################################
#
#  Select the rows of the null vectors  corresponding to the multiplication of xi*[x1^(d-1), x1^(d-1)x2,.......x1,x2,..., 1]
#
Choosexi:=proc(V, ord, d,xi)
   local n, N, Vt, MonVec, MonVec1, L, v, N1, i,vi;
     n:=nops(ord);
     N:=binomial(d+n,d);
     if RowDimension(V)<>N then
        lprint(`We have wrong null space`);
        RETURN([]);
     fi;
     Vt:=Transpose(V);
     MonVec:=numbve(d, n);
     MonVec1:=numbve(d-1, n);
     if  member(xi, ord, 'k') then     
       
        L:=Matrix(Column(Vt, N-(n-k+1)));   
        v:=convert(Column(IdentityMatrix(n),n-k+1),list);
        N1:=binomial(d-1+n,d-1);
        for i from 2 to N1 do 
            vi:=v+MonVec1[i];
           
            if member(vi, MonVec, 'kk') then 
             
              L:=Matrix([Column(Vt, N-kk+1), L]);
            fi;
        od;
        RETURN(Transpose(L));
      
     else
        lprint(`We have wrong variable`);
        RETURN([]);
     fi;
   

  end proc; 



##########################################################
#
#  Generate the multiplication matrix w.r.t xi and output the roots
#
Multimatrix:=proc(V, ord, d)
    local n, N, E1, E2, numsol, U1,S1, V1,B, M, MList, ev, sol,k;
       n:=nops(ord);
       N:=binomial(n+d-1,d-1);
       E1:=SubMatrix(V,-N..-1, 1..-1);
      E2:=SubMatrix(V, -binomial(n+d,d)..-1, 1..-1);
       MList:=[];
       numsol:=ColumnDimension(V);
        U1,S1,V1:=LinearAlgebra[SingularValues](E1, output=['U', 'S','Vt']):
            
            U1:=SubMatrix(U1,1..RowDimension(U1), 1..numsol):
            V1:=SubMatrix(V1, 1..numsol, 1..ColumnDimension(V1)):
            S1:=DiagonalMatrix(['S1[k]' $'k'=1..numsol]):
            sol:={};
            for k from 1 to n do 
            B:=Choosexi(E2, ord, d,ord[k]);        
            M:= HermitianTranspose(U1).B. HermitianTranspose(V1).MatrixInverse(S1):
            MList:=[op(MList),M];
             ev:=Trace(M)/numsol;
            sol:=[ev, op(sol)];
            od;

            return([numsol, sol,MList]);


  end proc; 

 
 
##############################################################################################################################
#
#  Select the rows of the null vectors  corresponding to the multiplication of xi*[x1^(d-1), x1^(d-1)x2,.......x1,x2,..., 1]
#
Choosexiexact:=proc(V, ord, d, Mon, xi)
   local n, N, Vt, MonVec, MonVec1, L, v, N1, i,vi;
     n:=nops(ord);
     N:=binomial(d+n,d);
     if RowDimension(V)<>N then
        lprint(`We have wrong null space`, RowDimension(V), N);
        RETURN([]);
     fi;
     Vt:=Transpose(V);
     MonVec:=numbve(d, n);
     MonVec1:=Mon;
      if  member(xi, ord, 'k') then     
        
        L:=Matrix(Column(Vt, N-(n-k+1)));   
        v:=convert(Column(IdentityMatrix(n),n-k+1),list);
        N1:=binomial(d-1+n,d-1);
        for i from 2 to nops(MonVec1) do 
            vi:=v+MonVec1[i];
           
            if member(vi, MonVec, 'kk') then 
             
              L:=Matrix([Column(Vt, N-kk+1), L]);
            fi;
        od;
        RETURN(Transpose(L));
      
     else
        lprint(`We have wrong variable`);
        RETURN([]);
     fi;
    
   

  end proc; 

 

##########################################################
#
#  Generate the multiplication matrix w.r.t xi and output the roots, for the exact case
#
Multimatrixexact:=proc(V, ord, d)
    local n, N, E1, E2, numsol, U1,S1, V1,B, M, MList, ev, sol,k,mon,B1,r,MonVec1,i,r1;
       n:=nops(ord);
       N:=binomial(n+d-1,d-1);
       E1:=SubMatrix(V,-N..-1, 1..-1) ;
       E2:=V;
       MList:=[];
       numsol:=ColumnDimension(V);
       B1:=Matrix(Row(E2,-1));    
       r:=Rank(B1);
       MonVec1:=numbve(d-1, n);
       mon:=[MonVec1[1]];
        for i from 2 to N do 
            r1:=Rank(Matrix([ [Row(E1,-i)],[B1]]));
           if r1>r  then
               B1:=Matrix([ [Row(E1,-i)], [B1]]);
               mon:=[op(mon), MonVec1[i]];
           fi;  
	   r:=r1;
               if r=numsol then      
		          sol:={};
                  for k from 1 to n do
                    B:=Choosexiexact(V, ord, d, mon, ord[k]);
                    M:=B.MatrixInverse(B1);
                    MList:=[M, op(MList)];
                     ev:=Trace(M)/numsol;
                    sol:=[ev, op(sol)];
                 od;
                RETURN([numsol, sol,MList]);             
               fi;              
        od;

end proc;
              
   


#####################################################
#
# Refine the multiple solution  near the origin.  
#
Polysolvenew:=proc(lsys, ord, tol)
   local L, U1,V1, S1,V;
       
      L:=indexmulti(lsys,ord,tol);
      U1, S1, V1:=SingularValues(Transpose(L[4]), output=['U', 'S', 'Vt']);
      V:=HermitianTranspose(SubMatrix(V1, -L[2]..-1,1..-1));
      
      RETURN( Multimatrix(V, ord, L[1]));
     end proc:





#########################################
#
# Refine the root by the orthogonal matrix, the second method in Wu-Zhi ISSAC2008.
#
#
Refineroot:=proc(lsys,ord, tol, mroot,it)
  local root, var, lsys1, newroot,i, ttol,k;
    root:=mroot;
    ttol:=tol;
    k:=0;
    while k<it do 
    var:={};
    for i from 1 to nops(ord) do 
     var:={op(var), ord[i]=ord[i]+root[i]};
    od;
    lsys1:=expand(subs(var,lsys));
    newroot:=Polysolvenew(lsys1,ord,ttol)[2];
    lprint(newroot+root);
    if Norm(Vector(newroot),2)<10^(-Digits+3) then
       lprint(`the number of iteration is`, k);
       RETURN(root+newroot);
    else
      root:=root+newroot;
      ttol:=10^(-2)*ttol;
      k:=k+1;
    fi;
    od;
       lprint(`the number of iteration is`, k);
       RETURN(root);
    end proc;




###########################################
#
#  Generate the matrix using Greg's subroutine Prolong
#
ProlongMatrix:=proc(lsys, ord, t)
 local dlsys, psys, dOrder, DRing, vars,A,b, deps, indeps, idep;
    
    dlsys:= PolynomialToPDE(lsys, ord, [U]);
    deps:=[U];
    indeps:=ord;
    psys:=Prolong(dlsys, ord, t);
    dOrder:=PDEtools[difforder](psys):
       DRing:= diffalg[differential_ring](derivations=ord, ranking=[U], notation=diff);
       psys:= map(z -> z = 0, reorder(map(diffalg[denote], psys, jet, DRing))): 
       vars:= DEtools[checkrank]([seq(deps[j](op(indeps)), j = 1 .. nops(deps))],
                                      deps, indep=indeps,degree=dOrder);
       vars:= map(diffalg[denote], vars, jet, DRing):
       # The jet variables at order i are vars[i]
       (A,b):= LinearAlgebra[GenerateMatrix](psys,vars):
      RETURN(A);
     end proc;
  


     
#################################################
# 
#  Compute the differential conditions at mroot
#  if tol<10^(-8), we use  exact routines.
#

DiffCond2:=proc(lsys,ord,  mroot, tol)
  local U1, S1, V1,V,L, L1, MonVec, j, i, N,m, DL,var,lsys1;

 #perform the linear transformation to move the mroot to the origin.
   var:=[];
   for i from 1 to nops(ord) do 
     var:={op(var), ord[i]=ord[i]+mroot[i]};
    od;
    lsys1:=expand(subs(var,lsys));
  #compute the multiplicity and index 
  L:=indexmulti(lsys1,ord,tol);
  
  if tol<10^(-8) then 
     # try nullspace for exact case
    V1:=Matrix([op(NullSpace(convert(Transpose(L[3]),rational)))]);
  
    if ColumnDimension(V1)=L[2] then 
      V:=V1;
    fi;
 else
  U1, S1, V1:=SingularValues(Transpose(L[3]), output=['U', 'S', 'Vt']);
  
  V:=HermitianTranspose(SubMatrix(V1, -L[2]..-1,1..-1));
 fi;

 
  MonVec:=numbve(L[1]-1,nops(ord));
  m:=L[2];
  N:=nops(MonVec);
  DL:=[];
 
  for j from 1 to m do 
     L1:=0;
    for i from 1 to N do 
      if abs(Column(V,j)[i])>10^(-Digits+6) then 
      L1:=L1+c[op(MonVec[N-i+1])]*Column(V,j)[i];    
      fi;
    end do;
     DL:=[op(DL),L1];
  od;
  lprint(`we get the differential condition`, nops(DL), DL);
  RETURN(Vector(DL));

end proc;

    
###############################################################################
# 
# Compute the local ring at origin, i.e., the multiplication matrices w.r.t an orthogonal basis
#
LocalRing2:=proc(lsys,ord, tol)
  local U1, S1, V1,V,L;



  #compute the multiplicity and index 
  L:=indexmulti(lsys,ord,tol);
   if tol<10^(-8) then 
     # try nullspace for exact case
    V1:=Matrix([op(NullSpace(convert(Transpose(L[4]),rational)))]);
   
    if ColumnDimension(V1)=L[2] then     
      RETURN( Multimatrixexact(V1, ord, L[1]));
    fi;
 else
  

  U1, S1, V1:=SingularValues(Transpose(L[4]), output=['U', 'S', 'Vt']);
  V:=HermitianTranspose(SubMatrix(V1, -L[2]..-1, 1..-1));
  RETURN( Multimatrix(V, ord, L[1]));
 fi;
 
  

end proc;