#####################################################################################################################################
## The following code is to compute the approximate GCD of two multivariate polynomials
##   Input--
##          F,G--the polynomials
##          ord-- the list of variables
##          e  -- the tolerance
##   Output--
##            a list= [v,da*L1,fbara,gbara]
##              where v is an integer describes if the backward error is big.
##                    da*L1 is the GCD
##                    fbara is the prime part of the polynomial F
##                    gbara is the prime part of  the polynomial G    
###################################################################################################################################  
multigcdew:= proc(F::polynom, G::polynom,ord,e )
  local
     n                # the number of variables
    ,QQ,RR            # QR decomposition
    ,SS,PP,KK         # the least singular value and the vector
    ,f,g,f0,g0        # polynomials expanded
    ,k                # degree of GCD
    ,i,kk,j           # loop variables
    ,Mg               # bvcauchyo  matrix
    ,Mf               # bvcauchyo  matrix
    ,Mfg              # the matrix constructed by Mf and Mg
    ,x0               #  the null space corrsponding to the least singular values
    ,gbara            # subvector of x0
    ,fbara            # subvector of x0 
    ,tmp              # degree of GCD
    ,da               # GCD
    ,len1             # the number terms in a  polynomial 
     ,gv	            # the vector converted by polynomial g
     ,fv              # the vector converted by polynomial f
     ,u_vec           # the vector converted by the prime part of f
     ,v_vec           # the vector converted by the prime part of g
     ,fv_leading      # the leading coefficient of the polynomial f
     ,gv_leading      # the leading coefficient of the polynomial g
     ,u_leading       # the leading coefficient of the polynomial u
     ,v_leading       # the leading coefficient of the polynomial v
     ,dv              # the solution of Mg.dv=gv
     , er             # the backward error
     ,L1              # special factors x, y
     ,fgv;            # the vector by fv and gv
     
 
     f:=expand(F); 
     n:=nops(ord);
     userinfo(4,multigcdew,print('f'=f));
     g:=expand(G);
     userinfo(4,multigcdew,print('g'=g));
     L1:=1;
     #test  the total degree of the GCD
     KK:=testtdeg(f,g,ord,e);
     userinfo(2,multigcdew,print('KK'=KK));
     if KK[2]=1 then  #There is a big jump, we believe we have estimated correctly  the total degree. 
	      k:=KK[1];
	      if k=0 then
	         RETURN([1,L1,F,G]);
	      fi;  
        (SS,QQ,RR):=solve_x0(k,f,g,ord,e);
         x0:=SS[2]; # x0 is cofactor of the gcd.
     else # check rank deficiency of  the kth multivariate Sylvester matrix      
       PP:=testk(KK[1],f,g,ord,e);
       k:=PP[1];
       x0:=PP[2]; # x0 is cofactor of the GCD. 
       if k=0 then #the gcd is constant
           RETURN([1,L1,F,G]);
       fi; 
     fi;  
       tmp:=degree(g)-k:
       if tmp=0 then
         er:=testcofactor(f,g,ord);
         da:=g/gv[1];
         fv:=convec(f,ord);     
         gv := convec(g,ord); 
         len1:=(n+tmp)!/(n!*tmp!);
         tmp:=degree(f)-k:
         v_vec:=x0[1..len1];
         u_vec:=x0[len1+1..-1];
         gbara:=veccoe(v_vec,ord);
         fbara:=-veccoe(u_vec,ord);
       else
         fv:=convec(f,ord);     
         gv := convec(g,ord); 
         len1:=(n+tmp)!/(n!*tmp!);
         v_vec:=x0[1..len1];
         u_vec:=x0[len1+1..-1];
         tmp:=degree(f)-k:
         fbara:=-veccoe(u_vec,ord);
         gbara:=veccoe(v_vec,ord);
       fi; 
           if n=2 then
              Mg:=bvcauchyo(gbara,ord,degree(g)-degree(gbara));
           else
              Mg := mvcauchyo(gbara,ord,degree(g)-degree(gbara));
           fi;
         
           if n=2 then
             Mf:=bvcauchyo(fbara,ord,degree(f)-degree(fbara));
           else 
            Mf:=mvcauchyo(fbara,ord,degree(f)-degree(fbara));
          fi;
          Mfg:=Matrix([[Mf],[Mg]]);
          fgv:=Vector([fv,gv]);
          dv:=LinearAlgebra:-LeastSquares(Mfg,fgv);
          da := veccoe(dv,ord);
          er:=LinearAlgebra[Norm](Mfg.dv-fgv,2);
          userinfo(1,multigcdew,print('backwarderror'=er));
          userinfo(3,multigcdew,print('da'=da));
          userinfo(2,multigcdew,print('deg_gcd'=degree(da), 'deg_g'=degree(g), 'deg_h'=degree(f)) );
          if (degree(da) > degree(g) and degree(da) > degree(f)) then
            userinfo(1,multigcdew,`Something went wrong.  Degree of the GCD is too big.`); 
	         RETURN([0,L1,f,g]);
          else
            if er < 10^(-e) or er<10^(-6)  then
             userinfo(1,multigcdew,`The estimated degree of GCD is correct, we get a good gcd.`); 
	             RETURN([1,da*L1,fbara,gbara]);
	         elif er<10^(-e+1) then
	            userinfo(1,multigcdew,`The estimated degree of GCD  may be  correct`); 
	             RETURN([2,da*L1,fbara,gbara]); 
	         elif er<10^(-e+2) and er<10^(-2) then 
	              userinfo(1,multigcdew,`The estimated degree of GCD  may be  wrong`);
	              RETURN([3,da*L1,fbara,gbara,ord]);
	         else
	              userinfo(1,multigcdew,` The estimated degree of GCD maybe wrong, we get a bad gcd.`); 
	               RETURN([0, da*L1,fbara,gbara]);
	         fi; 
          fi;
end proc: