#-------------------------------------------------------------------------------------------------------------------------
# author :  Zhengfeng Yang
# date   :  Jan. 16,  2006.
#--------------------------------------------------------------------------------------------------------------------------
##################################################################################################
## Description:
##     Given the polynomials L:=[f1,f2,....fn]  and given the integer k, we want to compute the $tilde{L}$=[\tilde{f1},\tilde{f2},....\tilde{fn}]$ 
##     such that degree(GCD(($\tilde{f1},\tilde{f2},....\tilde{fn})$)=k, and $\tilde{f1},\tilde{f2},....\tilde{fn} satisfy the linear constraints. If $f1,f2,....fn \in C$, then we want to find  \tilde{f1},\tilde{f2},....\tilde{fn} \in C field using the procedure "C_con_mulpoly".
##     If $f1,f2,....fn \in R$, we try to find the polynomials $\tilde{f1},\tilde{f2},....\tilde{fn} $\in R$ using the procedure "R_con_mulpoly". If it can not succeed, then we find the 
##     the polynomials $\tilde{f1},\tilde{f2},....\tilde{fn} $ \in C$ field using the procedure "C_con_mulpoly".
##     This procedure is to compute the $k$-th multiple polynomials Sylvester matrix $tilde(M)$ 
##     generated by $\tilde{f}$, $\tide{g}$, and rank deficiency of $tilde(M)$ is 1.      
## Input:
##         L=[f1,f2,....fn] ----    the list of  the polynomials
##         k                ----    degree of gcd of the polynomials;
##         e                ----    the tolerance 
##         Con_L            -----   the matrix of  linear constraints
##         DD               ----    the weight matrix 
## Output 
##     pp               ----     gcd of res_L
##     prim_list        ----      the list of the cofactors
##     res_L            ----     the list of the  polynomials which have non-trivial GCD and degree of GCD is k. 
###################################################################################################################################
con_mulpoly_gcd:=proc(L,k,e,Con_L,DD)
  local 
    res_M               # the result of  $k-$th Sylvester matrix
    ,res_L              # the result of  the list of the polynomials
    ,r                  # the number of the L
    ,vec_org            # the vector by L
    ,i                  # the parameter in for loop statement
    ,dd                 # the vector of the linear constraints 
    ,prim_list          # the list of the cofactors
    ,gcd_part           # the gcd of  res_L
    ,Col_dim            # the column dimension of the matrix
    ,_S;                # the singular values 
     
     r:=nops(L);
     vec_org:=con_vec(L[1]);
     for i from 1 to r-1 do
       vec_org:=Vector([vec_org,con_vec(L[i+1])]); 
     end do; 
     if args[5]='default' then
        dd:=Vector(LinearAlgebra:-Dimension(vec_org));
     else  
       dd:=LinearAlgebra:-Column(Con_L,-1);
     fi;
       if LinearAlgebra:-Norm(map(x->Im(x),V),infinity)>10^(-4) or LinearAlgebra:-Norm(map(x->Im(x),dd),infinity)>10^(-4)  then
         if nargs=6 then
           return C_con_mulpoly(L,k,e,Con_L,DD,args[6]);
         else
           return C_con_mulpoly(L,k,e,Con_L,DD);
         fi;    
       fi;    
    
      res_L,gcd_part,prim_list:=R_con_mulpoly(L,k,e,Con_L,DD,['ff','cofa']);
      res_M:=generate_syl(res_L,x,k);
      Col_dim:=LinearAlgebra:-ColumnDimension(res_M);
      _S:=LinearAlgebra:-SingularValues(res_M)[1..Col_dim];
      if _S[-1]<10^(-7) and _S[-2]>10^(-3) then 
        if nargs=6 then
          if args[6]='cofa' then
              return gcd_part,prim_list
          elif args[6]=['ff','cofa'] then
             return res_L,gcd_part,prim_list;  
          else
             return gcd_part;
          fi;
        else
           return gcd_part;
        fi;        
      else
        if nargs=6 then
           return C_con_mulpoly(L,k,e,Con_L,DD,args[6]);
        else
           return C_con_mulpoly(L,k,e,Con_L,DD);
        fi;    
      fi;    
 end proc;


########################################################################################################################### 
## Description:
##     Given the polynomial list L:=[f1,f2,....fn]  and given two integer k, we want to compute the $tilde{L}$=[\tilde{f1},\tilde{f2},....\tilde{fn}]$ 
##     such that degree($\tilde{f1},\tilde{f2},....\tilde{fn})$=k and $\tilde{f1},\tilde{f2},....\tilde{fn} satisfy the linear constraints.   
##     This procedure is to compute the $k$-th multiple polynomials Sylvester matrix $tilde(M)$ 
##     generated by $\tilde{f1},\tilde{f2},....\tilde{fn}]$, and rank deficiency of $tilde(M)$ is 1.  
##     Here, $\tilde{f1},\tilde{f2},....\tilde{fn}]$ in C field.
##     In C_con_mulpoly,  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 $k$th- multiple polynomials Sylvester matrix by f1,f2,....fn
##     B is a column of M
##     A is a remaining columns of M 
## Input:
##         L=[f1,f2,....fn] ----    the list of  the polynomials
##         k                ----    degree of gcd of the polynomials;
##         e                ----    the tolerance 
##         Con_L            -----   the matrix of  linear constraints
##         DD               ----    the weight matrix 
## Output 
##     pp               ----     gcd of res_L
##     prim_list        ----      the list of the cofactors
##     res_L            ----     the list of the  polynomials which have non-trivial GCD and degree of GCD is k. 
###################################################################################################################################

C_con_mulpoly:=proc(L,k,e,Con_L,DD)
   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
     ,pp                  # the vectors
     ,v_sub               # the subvector
     ,var_dim             # the dimension of the vector which represents the Sylvester matrix
     ,vec_org             # the vector by the L
     ,vec_res             # the vector by the list of the polynomials which we compute
     ,M,Ma                # the  Sylvester matrix 
     ,Con_M               # the matrix of the  linear constraints
     ,dd                  # the vector of the linear constraints 
     ,zeta                # the vectors
     ,zeta_R,zeta_I       # the real part and the image part of zeta
     ,Con_M_R,Con_M_I     # the real part and the image part of Con_M 
     ,dd_R,dd_I           # the real part and the image part of dd
     ,A,Aa,A_R,A_I        # the submatrices of the matrix M 
     ,int_vec             # the vector
     ,D_R,D_I             # the real part and the image part of the matrix DD
     ,xp,xq               # the vector generated by the term of  the cofactors  
     ,p0,pv               #  the cofators of ff and gg 
     ,zz                  # the vector
     ,zz_R,zz_I           # the real part and the image part of zz
     ,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  
     ,X_int               # the matrix
     ,M_tem               # the matrix
     ,UU,_SS,_VV          # the singular Values of the matrix
     ,num_S               # the integer
     ,M_new_inverse       # the inverse matrix
     ,elta_s              # the vector
     ,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
     ,A_A_R,A_A_I         # the matrix  
     ,A_A_sub1,A_A_sub2   # the submatrix
     ,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 
     ,res_M               # the matrix of the k-th sylvester matrix by the list of the polynomials
     ,LL                  # the list of the polynomials arranged by the degree
     ,LL_org              # the list of the polynomials
     ,f_var               # the free variables
     ,zeta_A_R,zeta_A_I   # the real part and the image part of zeta_A
     ,prim_list           # the list of the cofactors
     ,num_dim             # the number of the parameter
     ,p1                  # the vector constructed by the term
     ,gcd_part            # the gcd of  res_L                  
     ,res_L;              # the list of the polynomials which we compute
      
     r:=nops(L);
     deg_L:=[];
     for i from 1 to r do
       deg_L:=[op(deg_L),degree(L[i],x)];
     end do;
     member(min(op(deg_L)),deg_L, 'nd');
     if nd=r then
       LL:=L;
     elif nd=1 then
       LL:=[seq(L[i],i=2..r),L[1]];
     else
       LL:=[seq(L[i],i=1..nd-1),seq(L[i],i=nd+1..r),L[nd]];
     fi;     
       deg_L:=[];
       for i from 1 to r do
         deg_L:=[op(deg_L),degree(LL[i],x)];
       end do;    
       vec_org:=con_vec(LL[1]);
       for i from 1 to r-1 do
         vec_org:=Vector([vec_org,con_vec(LL[i+1])]); 
       end do; 
         
         # we get the matrix of linear constraints
         if args[4]='default' then
           Con_M:=LinearAlgebra:-IdentityMatrix(LinearAlgebra:-Dimension(vec_org));
           Con_M_R:=convert(map(x->Re(x),convert(Con_M,matrix)),Matrix);
           Con_M_I:=convert(map(x->Im(x),convert(Con_M,matrix)),Matrix);
           dd:=Vector(LinearAlgebra:-Dimension(vec_org));
         else  
           Con_M:=LinearAlgebra:-SubMatrix(Con_L,1..-1,1..-2);
           Con_M_R:=map(x->Re(x),Con_M);
           Con_M_I:=map(x->Im(x),Con_M);
           dd:=LinearAlgebra:-Column(Con_L,-1);
         fi;  
         
         var_dim:=sum('deg_L[j]','j'=1..r)+r;
         
          # we get the weight matrix.
          if  args[5]='default' then
              D_R:=LinearAlgebra:-IdentityMatrix(var_dim);
              D_I:=Matrix(var_dim,var_dim);
          else  
              D_R:=map(x->Re(x),DD);
              D_I:=convert(map(x->Im(x),convert(DD,matrix)),Matrix);
          fi;
       
          ## construct the $k$-th Sylvster matrix, and choose the column $b$ such that $A.X appxo = b$.    
          
          M:=generate_syl(LL,x,k);
          pp:= LinearAlgebra[SingularValues](M, output=['Vt']);
          vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
          v_sub:=LinearAlgebra:-SubVector(vv,1..degree(LL[1])-k+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);
          X:=LinearAlgebra:-LeastSquares(A,b,optimize);
          error_orgin:=LinearAlgebra:-Norm(M.vv,2);
          
          # We add complex randomized preconditioning of the inputs 
          
          if LinearAlgebra:-Norm(map(x->Im(x),dd+vec_org),infinity)<=10^(-4) then   
            int_vec:=[seq(LinearAlgebra[RandomVector][row](LinearAlgebra:-ColumnDimension(Con_M),generator=-0.5..0.5)[j],j=1..LinearAlgebra:-ColumnDimension(Con_M))];
            int_vec:=convert(int_vec,Vector);
            int_vec:=LinearAlgebra:-Norm(vec_org,2)*r*int_vec/(LinearAlgebra:-Norm(int_vec,2)*max(1,LinearAlgebra:-Norm(b-A.X,2)));
            LL_org:=Con_Poly(dd+vec_org+I*min(0.5*error_orgin,0.1)*Con_M.int_vec,LL); 
            La:=LL_org; 
            Ma:=generate_syl(LL_org+La,x,k);
          else 
            if error_orgin>0.1 then 
              int_vec:=[seq(LinearAlgebra[RandomVector][row](LinearAlgebra:-ColumnDimension(Con_M),generator=-0.5..0.5)[j],j=1..LinearAlgebra:-ColumnDimension(Con_M))];
              int_vec:=convert(int_vec,Vector);
              int_vec:=LinearAlgebra:-Norm(vec_org,2)*r*int_vec/(LinearAlgebra:-Norm(int_vec,2)*max(1,LinearAlgebra:-Norm(b-A.X,2)));
              LL_org:=Con_Poly(dd+vec_org+I*min(0.5*error_orgin,0.1)*Con_M.int_vec,LL); 
              La:=LL_org; 
              Ma:=generate_syl(LL_org+La,x,k);
            else   
              LL_org:=Con_Poly(dd+vec_org,LL);
              Ma:=generate_syl(LL_org,x,k);
            fi;  
          fi;
          
          #  we compute the initial vector $X$ via the Lagrangian multiplier method, then do the iteration
          pp:= LinearAlgebra[SingularValues](Ma, output=['Vt']);
          vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
          X_int:=Int_point_matrix(vv,LL,k).Con_M;
          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;
            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_con_mulpoly,print('error_orgin'=error_orgin));
            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); 
            X_R:=map(x->Re(x),X);
            X_I:=map(x->Im(x),X);
            f_var:=LinearAlgebra:-ColumnDimension(Con_M);
            zeta:=elta_s;
            zeta_R:=map(x->Re(x),zeta);
            zeta_I:=map(x->Im(x),zeta);
            ww:=evalf(10^(min(5-ceil(log10(0.01*error_orgin)),8)));
            Row_dim:=LinearAlgebra:-RowDimension(A);
            Col_dim:=LinearAlgebra:-ColumnDimension(A);
            P0_Matrix:=P_mul_gcd(LL,choo_k,k);
            XY:=Matrix(Row_dim,var_dim);
            I_matrix:=Matrix([D_R.Con_M_R-D_I.Con_M_I,-D_R.Con_M_I-D_I.Con_M_R]);
            I_matrix:=Matrix([[I_matrix],[Matrix([D_R.Con_M_I+D_I.Con_M_R,D_R.Con_M_R-D_I.Con_M_I])]]);
            zero_matrix:=Matrix(2*var_dim,2*Col_dim);
            delta_x:=1;
            zz:=Con_M.(zeta_R+I*zeta_I)+dd;
            zz_R:=map(x->Re(x),zz);
            zz_I:=map(x->Im(x),zz);
            dd_R:=map(x->Re(x),dd);
            dd_I:=map(x->Im(x),dd);
            f_A:=P0_Matrix.zz;
            f_A_R:=map(x->Re(x),f_A);
            f_A_I:=map(x->Im(x),f_A);
            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>300 then
             break;
           fi;  
         
           ## we construct the matrix $E$  and $Y_{k}$ in our paper, according to X and elta. 
         
           EE_R:=Con_mul_gcd(zz_R,LL,k);
           EE_R:=LinearAlgebra:-DeleteColumn(EE_R,choo_k);
           EE_I:=Con_mul_gcd(zz_I,LL,k);
           EE_I:=LinearAlgebra:-DeleteColumn(EE_I,choo_k);
           XY_R:=Y_mul_gcd(X_R,LL,choo_k,k);
           XY_I:=Y_mul_gcd(X_I,LL,choo_k,k);
           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;
           userinfo(3,C_con_mulpoly,print('num'=num));
           userinfo(3,C_con_mulpoly,print('R'=LinearAlgebra:-Norm(R_R+I*R_I,infinity)));
           
           ## 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.Con_M_R-XY_I.Con_M_I-P0_Matrix.Con_M_R,-(XY_R.Con_M_I+XY_I.Con_M_R-P0_Matrix.Con_M_I)]);
           XY_sub2:=Matrix([XY_I.Con_M_R+XY_R.Con_M_I-P0_Matrix.Con_M_I,XY_R.Con_M_R-XY_I.Con_M_I-P0_Matrix.Con_M_R]);
           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]]);
           D_elta:=Vector([ww*R_R,ww*R_I,-(D_R.zz_R-D_I.zz_I),-(D_R.zz_I+D_I.zz_R)]);
           userinfo(3,C_con_mulpoly,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_con_mulpoly,print('delta_x'=delta_x));  
           
           zeta_A_R:=LinearAlgebra:-SubVector(X_A_Y,2*Col_dim+1..2*Col_dim+f_var);
           zeta_A_I:=LinearAlgebra:-SubVector(X_A_Y,2*Col_dim+f_var+1..2*Col_dim+2*f_var);      
           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;
           zeta_R:=zeta_R+zeta_A_R;
           zeta_I:=zeta_I+zeta_A_I;
           zz_R:=Con_M_R.zeta_R-Con_M_I.zeta_I+dd_R;
           zz_I:=Con_M_R.zeta_I+Con_M_I.zeta_R+dd_I;
          
           ## compute the first column of E.
           f_A_R:=P0_Matrix.zz_R;
           f_A_I:=P0_Matrix.zz_I;
        end do;
          # end of the iteration
          
          userinfo(1,C_con_mulpoly,print(`number of the iteration` =num));
          userinfo(2,C_con_mulpoly,print(`check the distance between f and ff`, (LinearAlgebra:-Norm((D_R+I*D_I).(zz_R+I*zz_I),2))^2)); 
          X:=X_R+I*X_I;
          
          # we compute the cofactor vector X.
          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;
          vec_res:=vec_org+zz_R+I*zz_I;
          res_L:=Con_Poly(vec_res,LL);
          
          prim_list:=[];
          num_dim:=0; 
          for j from 1 to r do
            xp:=LinearAlgebra:-SubVector(X,num_dim+1..num_dim+degree(res_L[j])-k+1);
            xp:=xp*lcoeff(res_L[j])/xp[1];
            p0:=convert([seq(x^(degree(res_L[j])-k-i),i=0..degree(res_L[j])-k)],Vector);
            p1:=LinearAlgebra:-DotProduct(p0,xp,conjugate=false);
            prim_list:=[op(prim_list),p1];
            num_dim:=num_dim+degree(res_L[j])-k+1;
          end do;   
            gcd_part:=mul_polydiv(res_L,prim_list);
            if nd=r then
              res_L:=res_L;
              prim_list:=prim_list;
            elif nd=1 then
              res_L:=[res_L[r],seq(res_L[i],i=1..r-1)];
              prim_list:=[prim_list[r],seq(prim_list[i],i=1..r-1)];
           else
            res_L:=[seq(res_L[i],i=1..nd-1),res_L[r],seq(res_L[i],i=nd..r-1)];
            prim_list:=[seq(prim_list[i],i=1..nd-1),prim_list[r],seq(prim_list[i],i=nd..r-1)];
           fi;       
           if nargs=6 then
             if args[6]='cofa' then
               return gcd_part,prim_list;
             else
              return gcd_part;
             fi;
           else
             return pp;
           fi;         
end proc; 


########################################################################################################################### 
## Description:
##     Given the polynomial list L:=[f1,f2,....fn]  and given the integer k, we want to compute the $tilde{L}$=[\tilde{f1},\tilde{f2},....\tilde{fn}]$ 
##     such that degree(GCD($\tilde{f1},\tilde{f2},....\tilde{fn})$)=k, and and $\tilde{f1},\tilde{f2},....\tilde{fn} satisfy the linear constraints.    
##     This procedure is to compute the $k$-th multiple polynomials Sylvester matrix $tilde(M)$ 
##     generated by \tilde{f1},\tilde{f2},....\tilde{fn}]$, and rank deficiency of $tilde(M)$ is 1.  
##     
##     In R_con_mulpoly,  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 $k$th- multiple polynomials Sylvester matrix by f1,f2,....fn
##     B is a column of M
##     A is a remaining columns of M 
## Input:
##         L=[f1,f2,....fn] ----    the list of  the polynomials
##         k                ----    degree of gcd of the polynomials;
##         e                ----    the tolerance 
##         Con_L            -----   the matrix of  linear constraints
##         DD               ----    the weight matrix 
## Output 
##     pp               ----     gcd of res_L
##     prim_list        ----     the list of the cofactors
##     res_L            ----     the list of the  polynomials which have non-trivial GCD and degree of GCD is k. 
###################################################################################################################################

R_con_mulpoly:=proc(L,k,e,Con_L,DD)
   local 
      r                    # the number of the polynomials which we input
     ,M                    # the  Sylvester matrix 
     ,elta,pp              # the vectors
     ,v_sub                # the subvector
     ,var_dim              # the dimension of the vector which represents the Sylvester matrix
     ,vec_org              # the vector by the L
     ,vec_res              # the vector by the list of the polynomials which we compute
     ,LL                   # the list of the polynomials arranged by the degree
     ,res_L                # the list of the polynomials which we compute
     ,A                    # the submatrix of the matrix M 
     ,b                    #  the column of M
     ,X                    # the least squares solution of A.X=b 
     ,ww                   # the weight value
     ,DD_test              #  the weight matrix
     ,Con_M                # the matrix of the  linear constraints
     ,dd                   # the vector of the linear constraints 
     ,zeta                 # the vectors
     ,xp,xq                # the vector generated by the term of  the cofactors  
     ,p0,pv                #  the cofators of ff and gg 
     ,zz                   # the vector
     ,P0_Matrix            # the matrix such that $ P0_Matrix.elta$ is the column of M  
     ,P0_old               # the matrix 
     ,XY                   # the matrix
     ,delta_x,f_A          # the vectors 
     ,num                  # the number of the iterations
     ,EE,E                 # the matrices which has Sylvester structure
     ,i,kk,uu,u,v,vv,j,jj  # the parameter in for loop statement 
     ,M_M                  # the whole matrix
     ,M_sub,M_sub2         # the sub matrix of M_M
     ,X_int                # the matrix
     ,M_tem                # the matrix
     ,UU,_SS,_VV           # the singular values decomposition
     ,num_S                # the integer
     ,M_new_inverse        # the inverse matrix
     ,elta_s               # the vector             
     ,Row_dim              # the row  dimension of the matix
     ,Col_dim              # the column dimension of the matrix
     ,deg_L                # the list of degree of the polynomials    
     ,error_orgin          # the error of $b-A.X$.
     ,I_matrix             # the identity matrix
     ,zero_matrix          # the zero matrix
     ,R                    # the residual i.e. b-A.X
     ,A_A                  # the matrix which denote the change of the matrix A
     ,D_elta               # the vector 
     ,X_A_Y                # the solution of LS problem
     ,X_A                  # the vector which denote the change of the vector X
     ,elta_A               # the vector which denote the change of the vector elta
     ,LL_org               # the list of the polynomials
     ,f_var                # the free variables
     ,zeta_A               # the changed part of  zeta
     ,prim_list           # the list of the cofactors
     ,num_dim             # the number of the parameter
     ,p1                  # the vector constructed by the term
     ,gcd_part            # the gcd of  res_L                  
     ,res_M;               # the matrix of the k-th sylvester matrix by ff and gg

     r:=nops(L);
     deg_L:=[];
     for i from 1 to r do
       deg_L:=[op(deg_L),degree(L[i],x)];
     end do;
      member(min(op(deg_L)),deg_L, 'nd');
     if nd=r then
       LL:=L;
     elif nd=1 then
       LL:=[seq(L[i],i=2..r),L[1]];
     else
       LL:=[seq(L[i],i=1..nd-1),seq(L[i],i=nd+1..r),L[nd]];
     fi;     
       deg_L:=[];
       for i from 1 to r do
         deg_L:=[op(deg_L),degree(LL[i],x)];
       end do;
       var_dim:=sum('deg_L[j]','j'=1..r)+r; 
       
       ## construct the $k$-th Sylvster matrix, and choose the column $b$ such that $A.X appxo = b$.     
       
       M:=generate_syl(LL,x,k);
       pp:= LinearAlgebra[SingularValues](M, output=['Vt']);
       vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
       v_sub:=LinearAlgebra:-SubVector(vv,1..degree(LL[1])-k+1):
       v_sub:=convert(v_sub,list):
       member(max(op(map(x ->abs(x),v_sub))), map(x->abs(x),v_sub),'choo_k'):
       vec_org:=con_vec(LL[1]);
       for i from 1 to r-1 do
         vec_org:=Vector([vec_org,con_vec(LL[i+1])]); 
       end do; 
       
       # we get the matrix of linear constraints
       if args[4]='default' then
          Con_M:=LinearAlgebra:-IdentityMatrix(LinearAlgebra:-Dimension(vec_org)):
          dd:=Vector(LinearAlgebra:-Dimension(vec_org));
       else    
         Con_M:=LinearAlgebra:-SubMatrix(Con_L,1..-1,1..-2);
         dd:=LinearAlgebra:-Column(Con_L,-1);
       fi;
      
       LL_org:=Con_Poly(dd+vec_org,LL);
       M:=generate_syl(LL_org,x,k);
       pp:= LinearAlgebra[SingularValues](M, output=['Vt']);
       vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1);
       error_orgin:=LinearAlgebra:-Norm(M.vv,infinity);
       
       #  we compute the initial vector $X$ via the Lagrangian multiplier method, then do the iteration
       X_int:=Int_point_matrix(vv,LL,k).Con_M;
       
       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.M.vv;
           X:=-vv/vv[choo_k];
           X:=convert(LinearAlgebra:-DeleteRow(convert(X,Matrix),choo_k),Vector);
          
          
           M:=generate_syl(LL,x,k);
           A:=LinearAlgebra:-DeleteColumn(M,choo_k);
           b:=LinearAlgebra:-Column(M,choo_k);
           error_orgin:=LinearAlgebra:-Norm(b-A.X,infinity);
           userinfo(1,R_con_mulpoly,print('error_orgin'=error_orgin));
           f_var:=LinearAlgebra:-ColumnDimension(Con_M);
           #zeta:=Vector(f_var);
           zeta:=elta_s;
           ww:=evalf(10^(min(5-ceil(log10(error_orgin)),6)));
           P0_old:=P_mul_gcd(LL,choo_k,k);
           P0_Matrix:=P0_old.Con_M;
           
           # we get the weight matrix
           if args[5]='default' then
             DD_test:=LinearAlgebra:-IdentityMatrix(var_dim);
           else
             DD_test:=DD;
           fi;
             I_matrix:=DD_test.Con_M;   
       
             Row_dim:=LinearAlgebra:-RowDimension(A);
             Col_dim:=LinearAlgebra:-ColumnDimension(A);
             XY:=Matrix(Row_dim,var_dim);
             I_matrix:=DD_test.Con_M;
             zero_matrix:=Matrix(var_dim,Col_dim);
             delta_x:=1;
             zz:=Con_M.zeta+dd;
             f_A:=P0_old.zz;
             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>50 then
             break;
           fi;  
         
           ## we construct the matrix $E$  and $Y_{k}$ in our paper, according to X and elta. 
         
           EE:=Con_mul_gcd(zz,LL,k);
           EE:=LinearAlgebra:-DeleteColumn(EE,choo_k);
           XY:=Y_mul_gcd(X,LL,choo_k,k);
           XY:=XY.Con_M;
           XY:=XY-P0_Matrix;
           A_A:=A+EE;
           R:=b+f_A-A_A.X;
           userinfo(3,R_con_mulpoly,print('num'=num));
           userinfo(3,R_con_mulpoly,print('R'=LinearAlgebra:-Norm(R,infinity)));
           
           ## we construct the big matrix M
          
           M_sub:=ww*Matrix([A_A,XY]);
           M_sub2:=Matrix([zero_matrix,I_matrix]);
           M_M:=Matrix([[M_sub],[M_sub2]]);
           D_elta:=Vector([ww*R,-DD_test.zz]);
          
           
           ## we compute LS solution  of the matrix M_M.
             
           X_A_Y:=LinearAlgebra:-LeastSquares(M_M,D_elta);
           delta_x:=LinearAlgebra:-Norm(X_A_Y,infinity);
           userinfo(2,R_con_mulpoly,print('delta_x'=delta_x));  
           zeta_A:=LinearAlgebra:-SubVector(X_A_Y,Col_dim+1..-1);
           X_A:=LinearAlgebra:-SubVector(X_A_Y,1..Col_dim);
           X:=X+X_A;
           zeta:=zeta+zeta_A;
           zz:=Con_M.zeta+dd;      
           ## compute the first column of E.
           f_A:=P0_old.zz;    
           
        end do;
          # end of the iteration
          
          userinfo(1,R_con_mulpoly,print(`number of the iterative` =num));
          vec_org:=con_vec(LL[1]);
          for i from 1 to r-1 do
            vec_org:=Vector([vec_org,con_vec(LL[i+1])]); 
          end do;
          userinfo(2,R_con_mulpoly_gcd,print(`check the distance between f and ff`, LinearAlgebra:-Norm(elta,2))); 
          vec_res:=vec_org+zz;
          res_L:=Con_Poly(vec_res,LL);
          
          # we compute the cofactor vector.
          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;
          prim_list:=[];
          num_dim:=0; 
          for j from 1 to r do
            xp:=LinearAlgebra:-SubVector(X,num_dim+1..num_dim+degree(res_L[j])-k+1);
            xp:=xp*lcoeff(res_L[j])/xp[1];
            p0:=convert([seq(x^(degree(res_L[j])-k-i),i=0..degree(res_L[j])-k)],Vector);
            p1:=LinearAlgebra:-DotProduct(p0,xp,conjugate=false);
            prim_list:=[op(prim_list),p1];
            num_dim:=num_dim+degree(res_L[j])-k+1;
          end do;   
            gcd_part:=mul_polydiv(res_L,prim_list);
            if nd=r then
              res_L:=res_L;
              prim_list:=prim_list;
            elif nd=1 then
              res_L:=[res_L[r],seq(res_L[i],i=1..r-1)];
              prim_list:=[prim_list[r],seq(prim_list[i],i=1..r-1)];
            else
              res_L:=[seq(res_L[i],i=1..nd-1),res_L[r],seq(res_L[i],i=nd..r-1)];
              prim_list:=[seq(prim_list[i],i=1..nd-1),prim_list[r],seq(prim_list[i],i=nd..r-1)];
           fi;       
           if nargs=6 then
             if args[6]='cofa' then
               return gcd_part,prim_list;
             elif args[6]=['ff','cofa'] then
               return res_L,gcd_part,prim_list
             else
              return gcd_part;
             fi;
           else
             return pp;
           fi;         
end proc;
  
  

#############################################################################################################################
## This procedure  compute the matrix $Y_{k}$.
##Input: 
##       L       -       a vector which is the  LS solution of Ax=B.
##       LL      -       a list of the polynomials [f1,f2,....fr]; 
##       choo_k  -       the number of the column $b$ such that $A.X approx= b$. 
##       k       -        degree of gcd
##        
##Output
##      MM-- the matrix satisfies the condition: MM.X=E.elta    
#################################################################################################################################
Y_mul_gcd:=proc(L,LL,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;
    r:=nops(LL);
    deg_LL:=[];
    for i from 1 to r do
       deg_LL:=[op(deg_LL),degree(LL[i],x)];
    end do;
      bb_dim:=sum('deg_LL[jj]','jj'=1..r)-r*(k-1);
      bb:=vector(bb_dim);
      poly_list:=[];
      num_dim:=0;
      for i from 1 to r do
        ff:=0;
        m:=degree(LL[i])-k; 
        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;
        M_list1:=[];     
        for j from 1 to r-1 do   
          ff:=poly_list[j];
          gg:=x^(deg_LL[r]+1);
          MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
          M:=LinearAlgebra:-Transpose(MM1);
          MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[r]+1);
          M_list1:=[op(M_list1),MM_sub1];
        end do;
          M_list2:=[];  
          for v from 1 to r-1 do  
            ff:=poly_list[r];
            gg:=x^(deg_LL[v]+1);
            MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
            M:=LinearAlgebra:-Transpose(MM1);
            MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[v]+1);
            M_list2:=[op(M_list2),MM_sub1];  
          end do;
            M_sub2:=LinearAlgebra:-DiagonalMatrix(M_list2);
            if r=2 then
              M_sub1:=M_list1[1];
            else
              M_sub1:=M_list1[1];
              for j from 2 to  r-1 do
                M_sub1:=Matrix([[M_sub1],[M_list1[j]]]);
              end do;
            fi;        
              MM:=Matrix([M_sub2,M_sub1]);
              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 $Y_{k}$.
##Input: 
##       L       -       a vector which is the  LS solution of Ax=B.
##       LL      -       a list of the polynomials [f1,f2,....fr]; 
##       k       -        degree of gcd
##        
##Output
##      MM-- the matrix satisfies the condition: MM.X=M.elta    
#################################################################################################################################
Int_point_matrix:=proc(L,LL,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;
    r:=nops(LL);
    deg_LL:=[];
    for i from 1 to r do
       deg_LL:=[op(deg_LL),degree(LL[i],x)];
    end do;
      bb_dim:=sum('deg_LL[jj]','jj'=1..r)-r*(k-1);
      bb:=vector(bb_dim);
      poly_list:=[];
      num_dim:=0;
      for i from 1 to r do
        ff:=0;
        m:=degree(LL[i])-k; 
        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;
        M_list1:=[];     
        for j from 1 to r-1 do   
          ff:=poly_list[j];
          gg:=x^(deg_LL[r]+1);
          MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
          M:=LinearAlgebra:-Transpose(MM1);
          MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[r]+1);
          M_list1:=[op(M_list1),MM_sub1];
        end do;
          M_list2:=[];  
          for v from 1 to r-1 do  
            ff:=poly_list[r];
            gg:=x^(deg_LL[v]+1);
            MM1:=LinearAlgebra:-SylvesterMatrix(ff,gg,x);
            M:=LinearAlgebra:-Transpose(MM1);
            MM_sub1:=LinearAlgebra:-SubMatrix(M,1..-1,1..deg_LL[v]+1);
            M_list2:=[op(M_list2),MM_sub1];  
          end do;
            M_sub2:=LinearAlgebra:-DiagonalMatrix(M_list2);
            if r=2 then
              M_sub1:=M_list1[1];
            else
              M_sub1:=M_list1[1];
              for j from 2 to  r-1 do
                M_sub1:=Matrix([[M_sub1],[M_list1[j]]]);
              end do;
            fi;        
              MM:=Matrix([M_sub2,M_sub1]);
              vec_add:=L;   
              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
##  Suppose b is the choo_k-th column of k-th Sylvester matrix, 
##  M is the matrix we compute by the following procedure.
##  then  b=M.elta    
##  Input --
##            L       --  the list of the polynomials
##            choo_k  --   the number of the column $b$
##            k       --  the given integer    
##  Output--
##            M-- the matrix 
########################################################################################################################################
 
P_mul_gcd:=proc(L,choo_k,k)
  local r,deg_L,i,M,j;
    r:=nops(L);
    deg_L:=[];
    for i from 1 to r do
       deg_L:=[op(deg_L),degree(L[i])];
    end do;
    
    M:=Matrix(sum('deg_L[jj]','jj'=1..r-1)+(r-1)*(deg_L[r]-k+1),sum('deg_L[j]','j'=1..r)+r);
    for j from 1 to deg_L[r]+1 do
      M[j+choo_k-1,j+sum('deg_L[jjj]+1','jjj'=1..r-1)]:=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. 
##        L     -  the list of the polynomials
##        k     -  degree of gcd      
## Output--
##         the matrix M which is $k$-th Sylvester matrix.         
###########################################################################################################################################    
   Con_mul_gcd:=proc(elta,L,k)
     local r,bb,LL,num_dim,i,j,v,m,ff,M;
       r:=nops(L);
       bb:=vector(LinearAlgebra:-Dimension(elta));
       LL:=[];
       num_dim:=0;
       for i from 1 to r do
         ff:=0;
         m:=degree(L[i]); 
         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;
           LL:=[op(LL),ff];
       end do;
            
       M:=generate_syl(LL,x,k);
       M:=subs({'bb[j]=elta[j]' $ 'j' = 1..LinearAlgebra:-Dimension(elta)},M);
         return(M);
     end proc;      

###########################################################################################
##  The subroutine transforms linear constraints  into a matrix about many polynomials(r>=2)
  mul_trans_mat:=proc(LL,L,var)
  local i,j,num_dim,b,v,r,deg_L,res_L,res_var,va,var_dim,L_sub,vec_org,A,fg_org,L_sol,bb,XX,res_M,res_b,vv;
    r:=nops(LL);
     deg_L:=[];
     for i from 1 to r do
       deg_L:=[op(deg_L),degree(LL[i],x)];
     end do;
     member(min(op(deg_L)),deg_L, 'nd');
     if nd=r then
       res_L:=LL;
       res_var:=var;
     elif nd=1 then
       res_L:=[seq(LL[i],i=2..r),LL[1]];
       res_var:=[seq(var[i],i=2..r),var[1]];
     else
       res_L:=[seq(LL[i],i=1..nd-1),seq(LL[i],i=nd+1..r),LL[nd]];
       res_var:=[seq(var[i],i=1..nd-1),seq(var[i],i=nd+1..r),var[nd]];
     fi;     
       deg_L:=[];
       for i from 1 to r do
         deg_L:=[op(deg_L),degree(res_L[i],x)];
       end do;    
       for i from 1 to r do
         va[i]:=res_var[i];
       end do;
        
      var_dim:=sum('deg_L[j]','j'=1..r)+r;
      v:=[seq(vv[i],i=1..var_dim)];
      L_sub:=subs({'va[1][deg_L[1]+1-j]=vv[j]'$'j'=1..deg_L[1]+1},L);
      
      num_dim:=deg_L[1]+1;
      for i from 2 to r do
        
        L_sub:=subs({'va[i][deg_L[i]+1-j]=vv[j+num_dim]'$'j'=1..deg_L[i]+1},L_sub);
        num_dim:=num_dim+deg_L[i]+1;
      end do;
       
      A,b:=LinearAlgebra:-GenerateMatrix(L_sub,v);
      vec_org:=con_vec(res_L[1]);
       for i from 1 to r-1 do
         vec_org:=Vector([vec_org,con_vec(res_L[i+1])]); 
       end do;
       b:=A.vec_org-b; 
       L_sol:=LinearAlgebra:-LinearSolve(A,b,free='s');
       bb:=[seq(L_sol[i]=0,i=1..var_dim)];
       res_M,res_b:=LinearAlgebra:-GenerateMatrix(bb,[op(indets(bb))]);
       XX:=LinearAlgebra:-LeastSquares(res_M,res_b);
       res_b:=res_b-res_M.XX;
       res_M:=Matrix([res_M,res_b]);
end proc;