Computing a gcd in Q[x,y]

f  = 
      1 + 2  x + 2  x^2 + x^3 + x y + 2 x^2  y + 2 x^3  y + 
        x^4  y

g  =  -y + x^3*y - x*y^2 + x^4*y^2

extract gcd of coefficients

fcoeff  =  CoefficientList [f,  y]

<<  poliesGCD.m

a  =  GCDPolies [fcoeff[[1]],  fcoeff[[2]],  x]

MatrixForm [a [[1]]]

gcdfcoeff  =  PoliesExtendedGCD [fcoeff[[1]],  
          fcoeff[[2]],  x] [[1]]

fprimitiv  =  (fcoeff / gcdfcoeff) . {y^0,  
          y^1}

fprimitiv  =  Simplify [fprimitiv]

gcoeff  =  CoefficientList [g,  y]

a  =  GCDPolies [gcoeff[[2]],  gcoeff[[3]],  x]

MatrixForm [a[[1]]]

gcdgcoeff  =  PoliesExtendedGCD [gcoeff[[2]],  
          gcoeff[[3]],  x] [[1]]

gprimitiv  =  (gcoeff / gcdgcoeff) . {y^0,  
          y^1,  y^2}

gprimitiv  =  Simplify [gprimitiv]

Compute the gcds of the extracted parts

MatrixForm [GCDPolies [gcdfcoeff,  
          gcdgcoeff,  x] [[1]]]

d1  =  PoliesExtendedGCD [gcdfcoeff,  gcdgcoeff,  
          x] [[1]]

Compute the gcds of the primitiv parts

MatrixForm [GCDPolies [fprimitiv,  gprimitiv,  
          y] [[1]]]

d2  =  PoliesExtendedGCD [fprimitiv,  gprimitiv,  
          y] [[1]]

CoefficientList [d2,  y]

so d2 is already primitive

DD  =  d1  * d2

Check

PolynomialGCD [f,  g]

Expand [DD]