Q:

Factor bivariate polynomial over the complex numbers

This is very much like Factoring polynomials to factors involving complex coefficients except that I'm concerned about bivariate polynomials, not univariate polynomials. Take for example the polynomial
$$3 - 6 x^2 + 3 x^4 - 10 y^2 + 6 x^2 y^2 + 3 y^4$$
I'd like to factor that over the complex numbers. If I know where this came from, I can do Factor[poly, Extension -> {Sqrt[3]}] to find that the above equals
$$\tfrac13 (3 - 3 x^2 + 2 \sqrt3 y - 3 y^2) (3 - 3 x^2 - 2 \sqrt3 y - 3 y^2)$$
but to find this I have to know, a priori, that $\sqrt3$ might be a useful extension.
Even more important, that I haven't forgotten any additional extensions which would have allowed me to factor this even more. (The latter is easy in the case above, but hard if the degrees get higher.)
Edit: I found that this question can be answered using IrreduciblePolynomialQ[poly, Extension → All].
Is there any way I can compute a factorization of a polynomial from $\mathbb Z[x,y]$ into factors from $\mathbb C[x,y]$ (or equivalently $\bar{\mathbb Q}[x,y]$ i.e. with algebraic coefficients), without the use of any knowledge except for the polynomial itself?

A:

Here is an approach based on finding an approximate root, bumping to an approximate factor using GroebnerBasis, and resolving as an exact factor using RootApproximant.
poly = 3 - 6*x^2 + 3*x^4 - 10*y^2 + 6*x^2*y^2 + 3*y^4;

x0 = 11/7;
roots = y /. NSolve[poly /. x -> x0, WorkingPrecision -> 400];
root1 = First[roots];
fac = First[
   GroebnerBasis[{poly, (y - root1)^10, (x - x0)^17}, {y, x}, 
    MonomialOrder -> DegreeReverseLexicographic, 
    CoefficientDomain -> InexactNumbers]];
dtl = Chop[GroebnerBasis`DistributedTermsList[fac, {y, x}]];
newdtl = MapAt[RootApproximant, dtl, 
   Thread[{1, Range[Length[dtl[[1]]]], 2}]];
algfactor = GroebnerBasis`FromDistributedTermsList[newdtl]

(* Out[28]= -1 + x^2 + (2 y)/Sqrt[3] + y^2 *)

One can divide out by this factor, rinse, repeat. I'll skip that last and anyway I suspect it will not reduce further.
Factor[poly/algfactor, Extension -> Automatic]

(* Out[31]= -3 + 3 x^2 - 2 Sqrt[3] y + 3 y^2 *)

An explanation of the method, and the code from which I thoughtlessly cribbed, may be found here.

