##################################################################################################################################
##  The following procedure  compute the k-th sylvster matrix.  
##    Input--
##              L -- the list of the polynomials
##              x -- the variable
##              k -- the integer
##    Output--
##            MM --  the k-th Sylvester matrix of f and g
###################################################################################################################################    
generate_syl:=proc(L,x,k)
  local r,n,M_list1,M_list2,i,j,MM1,M,m,t,MM,MM_sub1,MM_sub2,M_sub1,M_sub2;
    r:=nops(L);
    m:=degree(L[r],x);
    M_list1:=[];
    M_list2:=[];
    for i from 1 to r-1 do
      MM1:=LinearAlgebra:-SylvesterMatrix(L[r],L[i],x);
      M:=LinearAlgebra:-Transpose(MM1);
      n:=degree(L[i],x);
      t:=n+1-k;
      MM:=LinearAlgebra:-DeleteRow(LinearAlgebra:-DeleteColumn(LinearAlgebra:-DeleteColumn(M,[m+t+1..n+m]),[t+1..n]),[m+t+1..n+m]);
      MM_sub1:=LinearAlgebra:-SubMatrix(MM,1..-1,1..n-k+1);
      MM_sub2:=LinearAlgebra:-SubMatrix(MM,1..-1,n-k+2..-1);
      M_list1:=[op(M_list1),MM_sub1];
      M_list2:=[op(M_list2),MM_sub2];
    end do;
      M_sub1:=LinearAlgebra:-DiagonalMatrix(M_list1);
      if r=2 then
        M_sub2:=M_list2[1];
      else
        M_sub2:=M_list2[1];
        for j from 2 to  r-1 do
          M_sub2:=Matrix([[M_sub2],[M_list2[j]]]);
        end do;
      fi;        
      MM:=Matrix([M_sub1,M_sub2]);
      return(MM);
  end proc; 



# This routine computes the approximate GCD of univariate polynomials.

mul_ugcd:=proc(L,x,e)
local n2,m2,f0,g0,n,m,MM1,j,gcd_part,num_dim,prim_list,r,M,cc,t,i,MM,PP,QQ,RR,x0,xp,p0,p1,q1,q2,xq,pv,p2;
  r:=nops(L);
  M:=generate_syl(L,x,1);
  cc:=LinearAlgebra:-SingularValues(evalf(M),output='list');
  t:=0;
  for i from LinearAlgebra:-ColumnDimension(M) to 1 by -1  do
    if evalf(cc[i]/cc[1])<10^(-e)  then 
      t:=t+1;
    else
      break;
    fi;   
  end do;
    if t=0 then 
      return 1,L;
    fi; 
      MM:=generate_syl(L,x,t);
      (PP,QQ,RR):=leastsv(MM,e);
      x0:=PP[2];
      prim_list:=[];
      num_dim:=0; 
      for j from 1 to r do
        xp:=LinearAlgebra:-SubVector(x0,num_dim+1..num_dim+degree(L[j])-t+1);
        xp:=xp*lcoeff(L[j])/xp[1];
        p0:=convert([seq(x^(degree(L[j])-t-i),i=0..degree(L[j])-t)],Vector);
        p1:=LinearAlgebra:-DotProduct(p0,xp,conjugate=false);
        prim_list:=[op(prim_list),p1];
        num_dim:=num_dim+degree(L[j])-t+1;
      end do;   
      gcd_part:=mul_polydiv(L,prim_list);
      return gcd_part, prim_list;
 end proc;

# This routine computes the approximate GCD of univariate polynomials.

mul_ugcd_M:=proc(L,X,t)
local n2,m2,f0,g0,n,m,MM1,j,gcd_part,num_dim,prim_list,r,M,cc,i,MM,PP,QQ,RR,x0,xp,p0,p1,q1,q2,xq,pv,p2;
  r:=nops(L);
  MM:=generate_syl(L,x,t);
  x0:=X;
  prim_list:=[];
  num_dim:=0; 
  for j from 1 to r do
    xp:=LinearAlgebra:-SubVector(x0,num_dim+1..num_dim+degree(L[j])-t+1);
    xp:=xp*lcoeff(L[j])/xp[1];
    p0:=convert([seq(x^(degree(L[j])-t-i),i=0..degree(L[j])-t)],Vector);
    p1:=LinearAlgebra:-DotProduct(p0,xp,conjugate=false);
    prim_list:=[op(prim_list),p1];
    num_dim:=num_dim+degree(L[j])-t+1;
  end do;   
      gcd_part:=mul_polydiv(L,prim_list);
      return gcd_part, prim_list;
 end proc;



#########################################################################################################################################
## This procedure  obtain  the list of  polynomials from  a vector, here we know the list of degree of the polynomials
## Input--
##        L--  the vector which is constructed by the coefficients of the polynomials
##        LL-- the list of the polynomials
## Output--
##         res_L-- the polynomials
##################################################################################################################################    
Con_Poly:=proc(L,LL)
  local m,n,i,j,ff,num_dim,poly_list,r;
    poly_list:=[];
    r:=nops(LL);
    num_dim:=0;
    for j from 1 to r do
      m:=degree(LL[j]); 
      ff:=0;
      for i from 0 to m do
        ff:=ff+L[i+1+num_dim]*x^(m-i);
      end do;
         num_dim:=num_dim+m+1;
         poly_list:=[op(poly_list),ff];
    end do;    
       return poly_list;
  end proc;
        




##################################################################################################################################
## According to the primary part LL, and the list of the polynomials L  we compute the polynomial $p$.
## i.e. we compute the polynomial $p$ such that $f=pu$ and $f=pv$.
## Input--
##        L --   the list of the polynomials
##        LL  --  the list of the prime part of L
## Output--
##         gcd_res-- gcd of L       
#############################################################################################################################
 mul_polydiv:=proc(L,LL)
   local  f_vec,u_dim,r,p_dim,
      A_1,A_2,i,j,M,V,X,xp,gcd_res;
      r:=nops(L);
      f_vec:=con_vec(L[1]);
      for i from 1 to r-1 do
        f_vec:=Vector([f_vec,con_vec(L[i+1])]);
      end do;
        p_dim:=degree(L[1])-degree(LL[1])+1;
        u_dim:=LinearAlgebra:-Dimension(con_vec(LL[1]));
        A_1:=Matrix(p_dim+u_dim-1,p_dim);
        for j from 1 to p_dim do
          A_1[j..j+u_dim-1,j]:=con_vec(LL[1]);
        end do;
          for i from 2 to r do
            u_dim:=LinearAlgebra:-Dimension(con_vec(LL[i]));
            A_2:=Matrix(p_dim+u_dim-1,p_dim);
            for j from 1 to p_dim do
              A_2[j..j+u_dim-1,j]:=con_vec(LL[i]);
            end do;                               
               A_1:=Matrix([[A_1],[A_2]]);  
          end do;       
            X:=LinearAlgebra:-LeastSquares(A_1,f_vec);
            xp:=convert([seq(x^(p_dim-1-i),i=0..p_dim-1)],Vector);
            gcd_res:=LinearAlgebra:-DotProduct(X,xp,conjugate=false);
              return gcd_res;
 end proc;       
  
##################################################################################################################################
## According to the primary part $u$ and $v$,we compute the polynomial $p$.
## i.e. we compute the polynomial $p$ such that $f=pu$ and $f=pv$.
## Input--
##        f,g -- the polynomials
##        u,v -- the prime part of f and g
## Output--
##         gcd_res-- the GCD of f and g       
#############################################################################################################################
 polydiv1:=proc(f,g,u,v)
   local  f_vec,g_vec,u_dim,v_dim,p_dim,
     A_1,A_2,i,j,M,V,X,xp,gcd_res;
     f_vec:=con_vec(f);
     g_vec:=con_vec(g);
     if evalf((f_vec[1]*g_vec[1])/(u[1]*v[1]))<0 then
       f_vec:=-f_vec;
     fi; 
       u_dim:=LinearAlgebra:-Dimension(u);
       v_dim:=LinearAlgebra:-Dimension(v);
       p_dim:=degree(f)-u_dim+2;
       A_1:=Matrix(p_dim+u_dim-1,p_dim);
       A_2:=Matrix(p_dim+v_dim-1,p_dim);
       for i from 1 to p_dim do
         A_1[i..i+u_dim-1,i]:=u;
       end do;                               
       for j from 1 to p_dim do
         A_2[j..j+v_dim-1,j]:=v;
       end do;
         M:=Matrix([[A_1],[A_2]]);
         V:=Vector([f_vec,g_vec]);
         X:=LinearAlgebra:-LeastSquares(M,V);
         xp:=convert([seq(x^(p_dim-1-i),i=0..p_dim-1)],Vector);
         gcd_res:=LinearAlgebra:-DotProduct(X,xp,conjugate=false);
              return gcd_res;
 end proc;       
  
   
 
  
##############################################################################################################################
## This procedure is to convert a polynomial to  a coefficient  vector.
## Input--
##        f- the polynomial
## Output-- the vector
##################################################################################################################################    
con_vec:=proc(f)
  local L,V,n;
    n:=degree(f);
    L:=[seq(coeff(f,x,n-i),i=0..n)];
    V:=convert(L,Vector);
    return V;
 end proc;
 
 
#Polynomial division, check the exact divisibility.
polydiv:=proc(f::polynom, g::polynom,x)
  local n,m,A,i,j,sol,cf,q;
    n:=degree(f);
    m:=degree(g);
    A:=Matrix(n+1,n-m+1);
    if n<m then 
      RETURN(0);
    fi;
      for i from 1 to n+1 do
        for j from 1 to n-m+1 do 
          if j<=i and  m+j-i>=0 then
            A[i,j]:=coeff(g,x,m+j-i);
          fi;
        od;
      od;
      cf:=Vector(n+1,i->coeff(f,x,n-i+1));
      sol:=LinearAlgebra[LeastSquares](A,cf);
      q:=add(sol[i]*x^(n-m-i+1),i=1..n-m+1);
      RETURN(q);
end proc:
  
 


 
# This routine computes the least singular value of the matrix and the associated singular vector
 leastsv:=proc(M,e)
 local QQ,RR,kk,x0,r,DD,bb1,dimen,QQ1,RR1,epsilon,i,bb2,bb,zz,rk,pp,sigma,dimC,M1,M2,b1,b2,z0,RK; 
   (QQ,RR,rk):=LinearAlgebra[QRDecomposition](M, output=['Q', 'R', 'rank'],fullspan);
   userinfo(1,leastsv,print('rk'=rk));
   dimC:=LinearAlgebra[ColumnDimension](M);
   RK:=LinearAlgebra:-Rank(M);
   if RK<>rk and RK<dimC then
       z0:=op(1,LinearAlgebra:-NullSpace(RR));
       z0:=z0/LinearAlgebra[Norm](z0,2);
       epsilon:=evalf(10^(-21));
       pp:=[epsilon,z0];
       RETURN(pp,QQ,RR);
   fi;    
   if rk<dimC then
     if RR[rk,rk]<10^(-e) or RR[rk,rk]<10^(-5) then
       z0:=op(1,LinearAlgebra:-NullSpace(RR));
       z0:=z0/LinearAlgebra[Norm](z0,2);
       epsilon:=evalf(10^(-21));
       pp:=[epsilon,z0];
       RETURN(pp,QQ,RR);
     else    
       M1:=LinearAlgebra[SubMatrix](RR,[1..rk],[1..rk]);
       M2:=LinearAlgebra[SubMatrix](RR,[1..rk],[rk+1..dimC]);
       x0:=LinearAlgebra[RandomVector](dimC-rk, generator= (rand(51..99)/100) );
     
       b1:=-LinearAlgebra[MatrixVectorMultiply](M2,x0);
       b1:=convert(b1,Vector);
       b2:=LinearAlgebra[BackwardSubstitute](M1,b1);
       z0:=Vector(dimC,[b2,x0]);
       z0:=z0/LinearAlgebra[Norm](z0,2);
       sigma:=LinearAlgebra[VectorNorm](LinearAlgebra[MatrixVectorMultiply](M,z0),2);
       epsilon:=evalf(10^(-21));
       pp:=[epsilon,z0];
       RETURN(pp,QQ,RR);
     fi;  
   fi;    
     kk:=LinearAlgebra[RowDimension](RR);
     x0:=LinearAlgebra[RandomVector](LinearAlgebra[ColumnDimension](RR), generator= (rand(-3..3)/10) );
     r:=evalf(LinearAlgebra[MatrixNorm](RR,infinity));
     i:=1;
     while i<3 do    # using iterative method to compute the least singular value.
       DD:=Matrix(kk+1,LinearAlgebra[ColumnDimension](RR),[[LinearAlgebra[Transpose](2*r*x0)],[RR]]);
       bb1:=LinearAlgebra[MatrixVectorMultiply](RR,x0);
       dimen:=LinearAlgebra[Dimension](bb1);
       bb2:=r*LinearAlgebra[Multiply](LinearAlgebra[Transpose](x0),x0)-r;
       bb:=Vector(dimen+1,[[bb2],bb1]);
       (QQ1,RR1):=LinearAlgebra[QRDecomposition](DD);
       zz:=LinearAlgebra[BackwardSubstitute](RR1,LinearAlgebra[MatrixVectorMultiply](LinearAlgebra[HermitianTranspose](QQ1),bb));
       x0:=x0-zz;
       epsilon:=evalf(LinearAlgebra[Norm](LinearAlgebra[MatrixVectorMultiply](RR,x0),2)/LinearAlgebra[Norm](x0,2));     
       i:=i+1;
      end do; 
        x0:=x0/LinearAlgebra[Norm](x0,2);
        pp:=[epsilon,x0];
        RETURN(pp,QQ,RR);
end proc;
 
#The algorithm generate a random polynomial with the same support as the input polynomial.
randp:=proc(f,x)
  local g, d,L,i,norm_f;
     g:=0;
     d:=degree(f,x)+1;
     norm_f:=norm(f,2);
     L:=LinearAlgebra[RandomVector][row](d,generator=-0.5..0.5);
     L:=convert(L,Vector);
     L:=L*norm_f/LinearAlgebra:-Norm(L,2);
     for i from 1 to d  do 
      # if coeff(f,x,i-1)<>0 then
       g:=g+L[i]*x^(i-1);
     #  fi:
     od;
     RETURN(g);
 end proc:
 
 
 
 
 
################################################################################################## 
## The following procedure is added in Dec. 13.
# 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:=(n+m)!/(n!*m!);
      for i from 0 to N-1 do
        v:=combinat[inttovec](i,m);
        S:=[op(S),v];
      end do;
      RETURN(S);
 end proc:

 # 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:=(n+m)!/(n!*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;  

# 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!);
       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 call the Cauchy matrix routine to generate the multivariate version
# of the Sylvester matrix of f and g.  
# The parameter m decreases the size of the matrix. m=1 gives the full 
# Sylvester-style matrix.
#  Generalize the algorithm bvsyo written by John P May

mvsylo:= proc(f,g,ord,m)
        local M3 ;
    if nops(ord)=2 then
       M3:=Matrix([bvcauchyo(f,ord,degree(g,{op(ord)})-m),bvcauchyo(g,ord,degree(f,{op(ord)})-m)]);
    else    
      M3:=Matrix([mvcauchyo(f,ord,degree(g,{op(ord)})-m),mvcauchyo(g,ord,degree(f,{op(ord)})-m)]);
    fi;
    RETURN(M3);
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.

# Note: GenerateMatrix is probably not the most efficient way to do this.
#  Written by John P May

bvcauchyo := proc(f, ord,deg)
        local i, j, td, unknowns, g, eqns, x,y,
              __v; # Protect the unknowns
        x:=ord[1];
        y:=ord[2];
        td := degree(f, {x,y});

        unknowns:=[]: g:=0;

        for i from deg to 0 by -1 do
        for j from 0 to i do
                        g := g+__v[j,i-j]*x^j*y^(i-j);
                unknowns := [op(unknowns), __v[j,i-j]];
        end do;
        end do;

        eqns := expand(f*g);
        eqns := convert(bv2veco(eqns,ord),list);
        eqns := map(x->x=0, eqns);

        LinearAlgebra[GenerateMatrix](eqns, unknowns)[1];

end proc;


# This routine call the Cauchy matrix routine to generate the bivariate version
# of the Sylvester matrix of f and g.  This also called the Macauley Matrix?
# The parameter m decreases the size of the matrix. m=1 gives the full 
# Sylvester-style matrix.
#  Written by John P May

bvsylo := proc(f,g,ord,m)
        local M, M1, M2, r, c,x,y;
    x:=ord[1];
    y:=ord[2];
    Matrix([bvcauchyo(f,ord,degree(g,{x,y})-m), bvcauchyo(g,ord,degree(f,{x,y})-m)]);
end proc;

# This routine converts an input bivariate polynomial into a vector
 #  Written by John P May
 
 bv2veco := proc(f,ord)
        local i, j, tmpcoeff, td, fvec,x,y;
        x:=ord[1];
        y:=ord[2];
        tmpcoeff:=map(PolynomialTools[CoefficientList], PolynomialTools[CoefficientList](f,x) , y);
        td := degree(f, {x,y});
        fvec := [];
        for i from td to 0 by -1 do
           for j from 0 to i do
                # Insert coeff of x^j*y^(i-j)
                        try
                                fvec := [op(fvec), tmpcoeff[j+1][i-j+1]];
                        catch "invalid subscript selector":
                                fvec := [op(fvec), 0];
                        end try;
           end do;
       end do;
       Vector(fvec);
end proc;


# This routine converts the output of bv2vec back into a bivariate polynomial
#  Written by John P May

vec2bvo := proc(v,ord)
        local i, j, l, deg, f,x,y;
        x:=ord[1];
        y:=ord[2];
        l := LinearAlgebra[Dimension](v);
        f := 0;
        i:=-1;
        for deg from 0 to l do
           for j from deg to 0 by -1 do
                        try
                                f := f + v[i]*x^(j)*y^(deg-j);
                                i := i - 1;
                        catch:
                                break;
                        end try;
            end do;
        end do;
        f;
end proc;


##################################################################################################################################
## According to the primary part LL, and the list of the polynomials L  we compute the polynomial $p$.
## i.e. we compute the polynomial $p$ such that $f=pu$ and $f=pv$.
## Input--
##        L --   the list of the polynomials, where the polynomials in C[x1,x2....xn].
##        LL  --  the list of the prime part of L
## Output--
##         gcd_res-- gcd of L       
#############################################################################################################################
 multi_div:=proc(f,g,ord)
   local  f_vec,n,M_sub,M,N,r,i,X,gcd_res,V,
     k;
      n:=nops(ord);
      f_vec:=convert(con_vec_fg([f],ord),Vector);
      k:=degree(f)-degree(g);
      if n=2 then
        M:=bvcauchyo(g,ord,k);
      else
        M:=mvcauchyo(g,ord,k);
      fi;  
          X:=LinearAlgebra:-LeastSquares(M,f_vec);
          N:=(k+n)!/(k!*n!);
          gcd_res:=0;
          for i from 1 to N do
            V:=combinat[inttovec](N-i,n);
            gcd_res:=gcd_res+X[i]*mul(ord[jj]^V[jj],jj=1..n);
          end do; 
              return gcd_res;
end proc;