#-------------------------------------------------------------------------------------------------------------------------
# author :  Zhengfeng Yang
# date   :  Jan. 16,  2006.
#--------------------------------------------------------------------------------------------------------------------------


########################################################################################################################### 
## Description:
##     Given the polynomial f  and given the integer k, we want to compute the $tilde(f)$ 
##     that has a root of multiplicity at least k. 
##     This procedure is to compute the $k$-th multiple polynomials Sylvester matrix $tilde(M)$ 
##     generated by \tilde{f},diff(\tilde{f},x),....diff(\tilde{f),x$(k-1)$, and rank deficiency of $tilde(M)$ is 1.  
##     Here, $tilde{f} \in C$.
##     In C_multiple_root,  M=[B,A], we want to compute the matrix res=[f+B,A+E] such that $f+B \in Range(A+E)$.
##     M is a Sylvester matrix by f,diff(f,x),....diff(f,x$(k-1))
##     B is a column of M
##     A is a remaining columns of M 
## Input:
##         f            ----     the polynomial
##         k            ----     a root of multiplicity at least k
##         e            ----     the tolerance  
## Output 
##         res_f        ----     the polynomial has a root of multiplicity at least k. 
##         X            ----     the vector by the cofactors of ff, and diff(ff,x)..diff(ff,x$(k-1))
###################################################################################################################################

C_multiple_root:=proc(f,k,e)
   local 
     
     r                    # the number of the polynomials which we input
     ,La                  # the list of the polynomials
     ,deg_L               # the list of degree of the polynomials 
     ,Row_dim             # the row dimension of the matrix
     ,Col_dim             # the column dimension of the matrix
     ,elta,pp             # the vectors
     ,v_sub               # the subvector
     ,fv                  # the vector by f
     ,var_dim             # the dimension of the vector which represents the Sylvester matrix
     ,M,Ma                # the  Sylvester matrix 
     ,elta_R,elta1_R      # the vector
     ,elta_I,elta1_I      # the vector
     ,A,Aa,A_R,A_I        # the submatrices of the matrix M 
     ,b,ba                # the column of M
     ,X_R,X_I             # the least squares solution of A.X=b 
     ,ww                  # the weight value
     ,P0_Matrix           # the matrix such that $ P0_Matrix.elta$ is the column of M  
     ,XY                  # the matrix
     ,XY_R,XY_I           # the matrix
     ,delta_x,f_A         # the vector 
     ,f_A_R,f_A_I         # the vector
     ,num                 # the number of the iterations
     ,EE,E                # the matrix
     ,EE_R,EE_I           # the matrix
     ,E_R,E_I             # the matrix
     ,i,kk,uu,u,v,vv,j,jj # the parameter in for loop statement 
     ,XYY                 # the matrix
     ,M_M                 # the whole matrix
     ,M_sub,M_sub2        # the sub matrix of M_M
     ,f0,g0,fa,ga         # arrange the polynomials f ang g
     ,m,n                 # the degree of f and g   
     ,error_orgin         # the error of $b-A.X$.
     ,I_matrix            # the identity matrix
     ,zero_matrix         # the zero matrix
     ,zero_matrix_sub     # the sumatrix of the zero matrix
     ,D_elta_sub          # the sub vector of the vector D_elta 
     ,R                   # the residual i.e. b-A.X
     ,R_R,R_I             # the real part and the image part  of the residual
     ,A_A                 # the matrix which denote the change of the matrix A
     ,X,b_R,b_I           # the vector
     ,elta_A_R,elta_A_I   # the vector
     ,A_A_R,A_A_I         # the matrix  
     ,A_A_sub1,A_A_sub2   # the submatrix
     ,X_int               # the matrix which compute the initial point
     ,M_tem                   # the matrix
     ,UU,_SS,_VV          # the singular values of the matrix
     ,num_S               # the number of the singular values
     ,M_new_inverse       # the inverse matrix
     ,elta_s              # the vector
     ,XY_sub,A_A_sub      # the submatrix
     ,XY_sub1,XY_sub2     # the submatrix
     ,D_elta              # the vector 
     ,X_A_Y               # the solution of LS problem
     ,X_A                 # the vector which denote the change of the vector X
     ,X_A_R,X_A_I         # the vector 
     ,elta_A              # the vector which denote the change of the vector elta
     ,L                   # the list of the polynomials [f,diff(f,x),...diff(f,x$(k-1))  
     ,vec_org             # the vector by the L
     ,vec_res             # the vector by the list of the polynomials which we compute
     ,res_f               # the polynomial which we compute
     ,res_M;              # the matrix of the  Sylvester matrix by f,diff(f,x)....diff(f,x$(k-1))
     
     m:=degree(f);
     
     ## we get the list of the polynomial f, diff(f,x),...diff(f,x$(k-1))
     L:=[];
     for i from 1 to k-1 do
       L:=[op(L),diff(f,x$i)];
     end do;
       L:=[op(L),f];              
       deg_L:=[];
       for i from 1 to k do
         deg_L:=[op(deg_L),degree(L[i],x)];
       end do;
     
     ## construct the  Sylvster matrix, and choose the column $b$ such that $A.X appxo = b$.    
       
     M:=generate_syl(L,x,1);
     pp:= LinearAlgebra[SingularValues](M, output=['Vt']);
     vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
     v_sub:=LinearAlgebra:-SubVector(vv,1..degree(L[1])):
     v_sub:=convert(v_sub,list):
     member(max(op(map(x ->abs(x),v_sub))), map(x->abs(x),v_sub),'choo_k'):
     A:=LinearAlgebra:-DeleteColumn(M,choo_k);
     b:=LinearAlgebra:-Column(M,choo_k);
     if LinearAlgebra:-SingularValues(evalf(A))[LinearAlgebra:-ColumnDimension(A)]<10^(-9) then
        userinfo(1,C_multiple_root,print(`the polynomial f  has a root of multiplicity at least k, no iteration is needed`));
         return(M);
     fi;  
     
     #  we compute the initial vector $X$ via the Lagrangian multiplier method, then do the iteration
     
     X:=LinearAlgebra:-LeastSquares(A,b,optimize);
     error_orgin:=LinearAlgebra:-Norm(b-A.X,infinity);
     fv:=con_vec(evalf(f));
     
     # We add complex randomized preconditioning of the inputs
     if LinearAlgebra:-Norm(map(x->Im(x),fv),infinity)<10^(-4) then
        fa:=expand(min(0.5*error_orgin/k,0.02/k)*I*randp(f,x));
        fa:=f+fa;
        La:=[];
        for i from 1 to k-1 do
          La:=[op(La),diff(fa,x$i)];
        end do;
          La:=[op(La),fa];
          Ma:=generate_syl(La,x,1);
          pp:= LinearAlgebra[SingularValues](Ma, output=['Vt']);
          vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
          if nargs>=4 and args[4]=monic then
            X_int:=Int_point_root(vv,f,k).monic_M(m);
          else
            X_int:=Int_point_root(vv,f,k);
          fi;
            M_tem:=X_int.LinearAlgebra:-HermitianTranspose(X_int);
            UU,_SS,_VV:=LinearAlgebra:-SingularValues(M_tem,output=['U','S', 'Vt']);
            num_S:=0;
            for i from 1 to LinearAlgebra:-ColumnDimension(M_tem) do
              if _SS[LinearAlgebra:-ColumnDimension(M_tem)-i+1]<min(10^(-4),error_orgin) then
                num_S:=num_S+1;
              else
                break;
              fi;
            end do;
              _SS:=Vector([map(x->1/x,LinearAlgebra:-SubVector(_SS,1..LinearAlgebra:-ColumnDimension(M_tem)-num_S)),Vector(num_S)]);
              _SS:=LinearAlgebra:-DiagonalMatrix(_SS,LinearAlgebra:-RowDimension(M_tem),LinearAlgebra:-RowDimension(M_tem));
              M_new_inverse:=LinearAlgebra:-HermitianTranspose(_VV).LinearAlgebra:-SubMatrix(_SS.LinearAlgebra:-HermitianTranspose(UU),1..LinearAlgebra:-ColumnDimension(M_tem),1..-1);
              elta_s:=-LinearAlgebra:-HermitianTranspose(X_int).M_new_inverse.Ma.vv;
     else
       if  error_orgin > 0.1 then
         fa:=expand(min(error_orgin,0.3)*I*randp(f,x)); 
         fa:=f+fa;
         La:=[];
         for i from 1 to k-1 do
           La:=[op(La),diff(fa,x$i)];
         end do;
           La:=[op(La),fa];
           Ma:=generate_syl(La,x,1);
           pp:= LinearAlgebra[SingularValues](Ma, output=['Vt']);
           vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
           if nargs>=4 and args[4]=monic then
             X_int:=Int_point_root(vv,f,k).monic_M(m);
           else
             X_int:=Int_point_root(vv,f,k);
           fi;
             M_tem:=X_int.LinearAlgebra:-HermitianTranspose(X_int);
             UU,_SS,_VV:=LinearAlgebra:-SingularValues(M_tem,output=['U','S', 'Vt']);
             num_S:=0;
             for i from 1 to LinearAlgebra:-ColumnDimension(M_tem) do
               if _SS[LinearAlgebra:-ColumnDimension(M_tem)-i+1]<min(10^(-4),error_orgin) then
                 num_S:=num_S+1;
               else
                 break;
               fi;
             end do;
                _SS:=Vector([map(x->1/x,LinearAlgebra:-SubVector(_SS,1..LinearAlgebra:-ColumnDimension(M_tem)-num_S)),Vector(num_S)]);
                _SS:=LinearAlgebra:-DiagonalMatrix(_SS,LinearAlgebra:-RowDimension(M_tem),LinearAlgebra:-RowDimension(M_tem));
                M_new_inverse:=LinearAlgebra:-HermitianTranspose(_VV).LinearAlgebra:-SubMatrix(_SS.LinearAlgebra:-HermitianTranspose(UU),1..LinearAlgebra:-ColumnDimension(M_tem),1..-1);
                elta_s:=-LinearAlgebra:-HermitianTranspose(X_int).M_new_inverse.Ma.vv;
       else
         Ma:=M;
         pp:= LinearAlgebra[SingularValues](Ma, output=['Vt']);
         vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
         if nargs>=4 and args[4]=monic then
             X_int:=Int_point_root(vv,f,k).monic_M(m);
         else
           X_int:=Int_point_root(vv,f,k);
         fi;
           M_tem:=X_int.LinearAlgebra:-HermitianTranspose(X_int);
           UU,_SS,_VV:=LinearAlgebra:-SingularValues(M_tem,output=['U','S', 'Vt']);
           num_S:=0;
           for i from 1 to LinearAlgebra:-ColumnDimension(M_tem) do
             if _SS[LinearAlgebra:-ColumnDimension(M_tem)-i+1]<min(10^(-4),error_orgin) then
               num_S:=num_S+1;
             else
               break;
             fi;
           end do;
             _SS:=Vector([map(x->1/x,LinearAlgebra:-SubVector(_SS,1..LinearAlgebra:-ColumnDimension(M_tem)-num_S)),Vector(num_S)]);
             _SS:=LinearAlgebra:-DiagonalMatrix(_SS,LinearAlgebra:-RowDimension(M_tem),LinearAlgebra:-RowDimension(M_tem));
             M_new_inverse:=LinearAlgebra:-HermitianTranspose(_VV).LinearAlgebra:-SubMatrix(_SS.LinearAlgebra:-HermitianTranspose(UU),1..LinearAlgebra:-ColumnDimension(M_tem),1..-1);
             elta_s:=-LinearAlgebra:-HermitianTranspose(X_int).M_new_inverse.Ma.vv;
       fi;
     fi;            
       X:=-vv/vv[choo_k];
       X:=convert(LinearAlgebra:-DeleteRow(convert(X,Matrix),choo_k),Vector);
       error_orgin:=LinearAlgebra:-Norm(b-A.X,infinity);  
       userinfo(1,C_multiple_root,print('error_orgin'=error_orgin));
       X_R:=map(x->Re(x),X);
       X_I:=map(x->Im(x),X);
       A_R:=map(x->Re(x),A);
       A_I:=map(x->Im(x),A);
       b_R:=map(x->Re(x),b);
       b_I:=map(x->Im(x),b);  
       if nargs>=4 and args[4]=monic then
         elta1_R:=map(x->Re(x),elta_s);
         elta1_I:=map(x->Im(x),elta_s);
         elta_R:=monic_V(elta1_R,m);
         elta_I:=monic_V(elta1_I,m);
       else
           elta:=elta_s;
           elta_R:=map(x->Re(x),elta);
           elta_I:=map(x->Im(x),elta); 
       fi;   
         ww:=evalf(10^(min(4-ceil(log10(min(error_orgin,0.1)*(k-1)*norm(f,2))),7)));
         Row_dim:=LinearAlgebra:-RowDimension(A);
         Col_dim:=LinearAlgebra:-ColumnDimension(A);
         if nargs>=4 and args[4]=monic then
           P0_Matrix:=P_mul_root_monic(f,choo_k,k);
           XY:=Matrix(Row_dim,m);
           I_matrix:=evalf(LinearAlgebra:-IdentityMatrix(m));
           I_matrix:=LinearAlgebra:-DiagonalMatrix([(I_matrix)$2]);
           zero_matrix:=Matrix(2*(m),2*Col_dim); 
         else
           P0_Matrix:=P_mul_root(f,choo_k,k);
           XY:=Matrix(Row_dim,m+1);
           I_matrix:=evalf(LinearAlgebra:-IdentityMatrix(m+1));
           I_matrix:=LinearAlgebra:-DiagonalMatrix([(I_matrix)$2]);
           zero_matrix:=Matrix(2*(m+1),2*Col_dim);
         fi; 
         if nargs>=4 and args[4]=monic then
           f_A_R:=P0_Matrix.elta1_R;  
           f_A_I:=P0_Matrix.elta1_I;
         else   
            f_A_R:=P0_Matrix.elta_R;  
            f_A_I:=P0_Matrix.elta_I;
         fi;  
           delta_x:=1;
           num:=0;
     
     
     # the following, we compute the $\delta f$ and $\delta E$ such that
     # $b+\delta f \in Range(A+\delta E)$
  
     while delta_x>10^(-e)  do
           num:=num+1;
           if num>1000 then
             break;
           fi;  
         
           ## we construct the matrix $E$  and $Y_{k}$ in our paper, according to X and elta. 
         
           EE_R:=Con_mul_root(elta_R,f,k);
           EE_I:=Con_mul_root(elta_I,f,k);
           EE_I:=LinearAlgebra:-DeleteColumn(EE_I,choo_k);
           EE_R:=LinearAlgebra:-DeleteColumn(EE_R,choo_k);
           XY_R:=Y_mul_root(X_R,f,choo_k,k);
           XY_I:=Y_mul_root(X_I,f,choo_k,k);
           if nargs>=4 and args[4]=monic then
             XY_R:=XY_R.monic_M(m);
             XY_I:=XY_I.monic_M(m);
           fi; 
             A_A_R:=A_R+EE_R;
             A_A_I:=A_I+EE_I;
             R_R:=b_R+f_A_R-A_A_R.X_R+A_A_I.X_I;
             R_I:=b_I+f_A_I-A_A_R.X_I-A_A_I.X_R;
 
           ## we construct the big matrix M
           A_A_sub1:=Matrix([A_A_R,-A_A_I]);
           A_A_sub2:=Matrix([A_A_I,A_A_R]);
           A_A_sub:=Matrix([[A_A_sub1],[A_A_sub2]]);
           XY_sub1:=Matrix([XY_R-P0_Matrix,-XY_I]);
           XY_sub2:=Matrix([XY_I,XY_R-P0_Matrix]);
           XY_sub:=Matrix([[XY_sub1],[XY_sub2]]);
           M_sub:=ww*Matrix([A_A_sub,XY_sub]);
           M_sub2:=Matrix([zero_matrix,I_matrix]);
           M_M:=Matrix([[M_sub],[M_sub2]]);
           if nargs>=4 and args[4]=monic then
             D_elta:=Vector([ww*R_R,ww*R_I,-elta1_R,-elta1_I]);
           else   
             D_elta:=Vector([ww*R_R,ww*R_I,-elta_R,-elta_I]);
           fi;
           userinfo(3,C_multiple_root,print('num'=num));
           
           ## we compute LS solution  of the matrix M_M.
           X_A_Y:=LinearAlgebra:-LeastSquares(M_M,D_elta,optimize);
           delta_x:=LinearAlgebra:-Norm(X_A_Y,infinity);
           userinfo(2,C_multiple_root,print('delta_x'=delta_x));  
           if nargs>=4 and args[4]=monic then
             elta_A_R:=LinearAlgebra:-SubVector(X_A_Y,2*Col_dim+1..2*Col_dim+m);
             elta_A_I:=LinearAlgebra:-SubVector(X_A_Y,2*Col_dim+m+1..2*Col_dim+2*m);
           else
             elta_A_R:=LinearAlgebra:-SubVector(X_A_Y,2*Col_dim+1..2*Col_dim+m+1);
             elta_A_I:=LinearAlgebra:-SubVector(X_A_Y,2*Col_dim+m+2..2*Col_dim+2*m+2);      
           fi;
             X_A_R:=LinearAlgebra:-SubVector(X_A_Y,1..Col_dim);
             X_A_I:=LinearAlgebra:-SubVector(X_A_Y,Col_dim+1..2*Col_dim);
             X_R:=X_R+X_A_R;
             X_I:=X_I+X_A_I;
             
             ## compute the column of E.
             if nargs>=4 and args[4]=monic then
               elta1_R:=elta1_R+elta_A_R;
               elta_R:=monic_V(elta1_R,m);
               elta1_I:=elta1_I+elta_A_I;
               elta_I:=monic_V(elta1_I,m,n);
               f_A_R:=P0_Matrix.elta1_R;
               f_A_I:=P0_Matrix.elta1_I;
             else
               elta_R:=elta_R+elta_A_R;
               elta_I:=elta_I+elta_A_I;
               f_A_R:=P0_Matrix.elta_R;
               f_A_I:=P0_Matrix.elta_I;
             fi; 
          # userinfo(3,C_multiple_root,print(`check the distance between f and ff`, LinearAlgebra:-Norm(elta_R+I*elta_I,2))); 
     end do;
           # end of the iteration
          
           userinfo(1,C_multiple_root,print(`number of the iteration` =num));
           vec_org:=con_vec(f);
           userinfo(1,C_multiple_root,print(`check the distance between f and ff`, (LinearAlgebra:-Norm(elta_R+I*elta_I,2))^2)); 
           vec_res:=vec_org+elta_R+I*elta_I;
           res_f:=op(Con_Poly(vec_res,[f]));
           
           # compute the cofactor vector.
           if args[-1]='cofa' then
             X:=X_R+I*X_I;
             if choo_k=1 then
                X:=Vector([-1,X]);
             else
               X:=Vector([LinearAlgebra:-SubVector(X,1..choo_k-1),-1,LinearAlgebra:-SubVector(X,choo_k..-1)]);   
             fi;
             return res_f,X;
           else
             return res_f;
           fi;    
end proc;

#############################################################################################################################
## This procedure is to used to compute the matrix $Y_{k}$.
##Input: 
##       L       -       a vector which is the  LS solution of Ax=B.
##       f      -        the polynomial ; 
##       choo_k  -       the number of the column $b$ such that $A.X approx= b$. 
##        k      -        the number of multiple roots.
##Output
##      MM-- the matrix satisfies the condition: MM.X=E.elta.    
#################################################################################################################################
Y_mul_root:=proc(L,f,choo_k,k)
  local r,deg_LL,bb_dim,bb,poly_list,ff,gg,m,u,v,num_dim,M_list1,MM1
    ,M,M_sub1,MM_sub1,M_list2,vec_add,i,j,M_sub2,MM,diff_vec,XY_diff,LL;
    deg_LL:=[];
    LL:=[];
    for i from 1 to k-1 do
      LL:=[op(LL),diff(f,x$(i))];
    end do;
      LL:=[op(LL),f];
      for i from 1 to k do
        deg_LL:=[op(deg_LL),degree(LL[i],x)];
      end do;
        bb_dim:=sum('deg_LL[jj]','jj'=1..k);
        bb:=vector(bb_dim);
        poly_list:=[];
        num_dim:=0;
        for i from 1 to k do
          ff:=0;
          m:=degree(LL[i],x)-1; 
          for v from 0 to m do
            ff:=ff+bb[num_dim+v+1]*x^(m-v);
          end do;
            num_dim:=num_dim+m+1;
            poly_list:=[op(poly_list),ff];
        end do;
        
        # we compute the cauchy matrix by $diff(f,x)$,....$diff(f,x$(k-1))$.
        M_list1:=[];     
        for j from 1 to k-1 do   
          ff:=poly_list[k];
          gg:=x^(deg_LL[k]+1);
          MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
          M:=LinearAlgebra:-Transpose(MM1);
          MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[k]+1);
          diff_vec:=Vector([seq((degree(f,x)-i)!/(degree(f,x)-i-j)!,i=0..degree(f,x)-j)]);
          diff_vec:=Vector([diff_vec,Vector(j)]);
          XY_diff:=LinearAlgebra:-DiagonalMatrix(diff_vec);
          MM_sub1:=LinearAlgebra:-SubMatrix(MM_sub1.XY_diff,1..2*deg_LL[k]-j,1..-1);
          M_list1:=[op(M_list1),MM_sub1];
        end do;
          
          ## we compute the Cauchy matrix by $f$
          M_list2:=[]; 
          for v from 1 to k-1 do  
            ff:=poly_list[v];
            gg:=x^(deg_LL[k]+1);
            MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
            M:=LinearAlgebra:-Transpose(MM1);
            MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[k]+1);
            M_list2:=[op(M_list2),MM_sub1];  
          end do;
           
            if k=2 then
              MM:=M_list1[1]+M_list2[1];
            else
              MM:=M_list1[1]+M_list2[1];
              for j from 2 to  k-1 do
                MM:=Matrix([[MM],[M_list1[j]+M_list2[j]]]);
              end do;
            fi;        
              if choo_k=1 then
                vec_add:=Vector([0,L]);
              else
                vec_add:=Vector([LinearAlgebra:-SubVector(L,1..choo_k-1),0,LinearAlgebra:-SubVector(L,choo_k..-1)]);
              fi;           
                 MM:=subs({'bb[j]=vec_add[j]'$'j'=1..bb_dim},MM);
                 return MM;
  end proc;  
 
 

#######################################################################################################################################
##  This procedure compute the matrix P0  as defined in the paper
##  In this case,  suppose b is the  choo_k column of  Sylvester matrix, 
##  M is the matrix we compute by the following procedure.
##  then  b=M.elta    
##  Input --
##          f       --  the polynomial
##          choo_k  --  the number of the column which satisfies that $Ax=b$
##          k     --  the number of the multiple roots of f
##  Output--
##            M-- the matrix 
########################################################################################################################################
 
P_mul_root:=proc(f,choo_k,k)
  local r,deg_L,i,M,j,m;
    m:=degree(f);
    M:=Matrix(2*(k-1)*m-k*(k-1)/2,m+1);
    for j from 1 to m+1 do
      M[j+choo_k-1,j]:=1;
    end do;
      return(M)
end proc;  


#######################################################################################################################################
##  This procedure  compute the matrix P0 as defined in the paper.  
##   Unlike the above procedure, in this case, the leading coefficient can not be perturbed.
##   Suppose b is the first column of k-th Sylvester matrix, 
##   M is the matrix we compute by the following procedure.
##   then  b=M.elta    
##  Input --
##          f,g --  the polynomials
##            k --  the given integer    
##  Output--
##            M-- the matrix 
########################################################################################################################################
 
P_mul_root_monic:=proc(f,choo_k,k)
  local n,m,j,M;
   m:=degree(f);
   M:=Matrix(2*(k-1)*m-k*(k-1)/2,m);
    for j from 1 to m do
      M[j+choo_k,j]:=1;
    end do;
      return(M)
end proc;



####################################################################################################################################
## The following procedure  compute  the matrix $E$ as defined in the paper, according to the vector L and [m,n,k].
## Input--
##        elta  - the vector represents free elements of a Sylvester matrix. 
##        f     -   the polynomial
##        k     -   the number of the multiple root      
## Output--
##        the matrix M which is  Sylvester matrix.         
###########################################################################################################################################    
   Con_mul_root:=proc(elta,f,k)
     local r,bb,LL,num_dim,i,j,v,m,ff,M;
       m:=degree(f);
       bb:=vector(LinearAlgebra:-Dimension(elta));
       LL:=[];
       ff:=add(bb[i+1]*x^(m-i),i=0..m);
       LL:=[]; 
       for j from 1 to k-1 do
         LL:=[op(LL),diff(ff,x$(j))];
       end do;
         LL:=[op(LL),ff];     
         M:=generate_syl(LL,x,1);
         M:=subs({'bb[j]=elta[j]' $ 'j' = 1..LinearAlgebra:-Dimension(elta)},M);
           return(M);
     end proc; 
     
     
# the procedure is to compute a simple matrix corresponding to a monic polynomial
monic_M:=proc(m)
  local i,j,M;
     M:=Matrix(m+1,m);
     for i from 1 to m do
       M[i+1,i]:=1;
     end do;
     return M;
  end proc;     
              


###########################################################################################
monic_V:=proc(V,m)
  local v1,VV;
    v1:=LinearAlgebra:-SubVector(V,[1..m]);
    VV:=Vector(m+1,[0,v1]);
    return VV;
  end proc;
       
#############################################################################################################################
## This procedure is to used to compute the matrix $H(c)$.
##Input: 
##       L       -       a vector which is the  LS solution of Ax=B.
##       f      -        the polynomial ; 
##        k      -        the number of multiple roots.
##Output
##      MM-- the matrix satisfies the condition: MM.X=M.elta.    
#################################################################################################################################
Int_point_root:=proc(L,f,k)
  local r,deg_LL,bb_dim,bb,poly_list,ff,gg,m,u,v,num_dim,M_list1,MM1
    ,M,M_sub1,MM_sub1,M_list2,vec_add,i,j,M_sub2,MM,diff_vec,XY_diff,LL;
    deg_LL:=[];
    LL:=[];
    for i from 1 to k-1 do
      LL:=[op(LL),diff(f,x$(i))];
    end do;
      LL:=[op(LL),f];
      for i from 1 to k do
        deg_LL:=[op(deg_LL),degree(LL[i],x)];
      end do;
        bb_dim:=sum('deg_LL[jj]','jj'=1..k);
        bb:=vector(bb_dim);
        poly_list:=[];
        num_dim:=0;
        for i from 1 to k do
          ff:=0;
          m:=degree(LL[i],x)-1; 
          for v from 0 to m do
            ff:=ff+bb[num_dim+v+1]*x^(m-v);
          end do;
            num_dim:=num_dim+m+1;
            poly_list:=[op(poly_list),ff];
        end do;
        
        # we compute the cauchy matrix by $diff(f,x)$,....$diff(f,x$(k-1))$.
        M_list1:=[];     
        for j from 1 to k-1 do   
          ff:=poly_list[k];
          gg:=x^(deg_LL[k]+1);
          MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
          M:=LinearAlgebra:-Transpose(MM1);
          MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[k]+1);
          diff_vec:=Vector([seq((degree(f,x)-i)!/(degree(f,x)-i-j)!,i=0..degree(f,x)-j)]);
          diff_vec:=Vector([diff_vec,Vector(j)]);
          XY_diff:=LinearAlgebra:-DiagonalMatrix(diff_vec);
          MM_sub1:=LinearAlgebra:-SubMatrix(MM_sub1.XY_diff,1..2*deg_LL[k]-j,1..-1);
          M_list1:=[op(M_list1),MM_sub1];
        end do;
          
          ## we compute the Cauchy matrix by $f$
          M_list2:=[]; 
          for v from 1 to k-1 do  
            ff:=poly_list[v];
            gg:=x^(deg_LL[k]+1);
            MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
            M:=LinearAlgebra:-Transpose(MM1);
            MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[k]+1);
            M_list2:=[op(M_list2),MM_sub1];  
          end do;
           
            if k=2 then
              MM:=M_list1[1]+M_list2[1];
            else
              MM:=M_list1[1]+M_list2[1];
              for j from 2 to  k-1 do
                MM:=Matrix([[MM],[M_list1[j]+M_list2[j]]]);
              end do;
            fi;        
                 vec_add:=L;
                 MM:=subs({'bb[j]=vec_add[j]'$'j'=1..bb_dim},MM);
                 return MM;
  end proc;