This is far from a complete answer but it does contain parts of a general method. The idea is to do it numerically, using a univariate factorization with specific values plugged in for two variables, then "lifting" one variable at a time to a full factorization. The list step is not very efficient (it's done with Groebner basis computations) and is made more effective by trial-and-error choice of finite precision. Among other drawbacks this means that algebraic number coefficients will need to be recovered from the approximate numeric coefficients.
Most of the method is described here. One incomplete aspect is the lifting degree. To be correct one should (i) make a "random" linear change of variables and (ii) determine how far to lift for the case of more than two variables (the paper only really covers the bivariate case). So this can be regarded as unfinished business. My suspicion is that either the later bounds become unwieldy or else are as low as the bivariate lift bound, but I have not sat down to try to prove which is the case, or check just how the bivariate proof of lifting bount might extend to more variables.
So here we go with some code.
algpoly =
4 (3 I Sqrt[3] (-12 + 6 x - 4 x^2 + x^3) y^3 z +
9 (12 - 12 x + 12 x^2 - 6 x^3 + x^4) y^4 z^2 +
I Sqrt[3] (8 + x^3) y z^3 + (4 - 2 x + x^2)^2 z^6 -
3 y^2 (4 + 24 z^4 - 8 x^3 z^4 + x^4 z^4 - 2 x (1 + 12 z^4) +
x^2 (1 + 22 z^4)));
We don't really use these (that would be cheating...) but we will show that we recover one of them, up to a constant numeric factor.
fac1 = 2 (-3 I - Sqrt[3] + Sqrt[3] x) y -
3 (6 I - 2 Sqrt[3] -
2 (3 I + Sqrt[3]) x + (I + Sqrt[3]) x^2) y^2 z + (8 I -
4 (I + Sqrt[3]) x + (-I + Sqrt[3]) x^2) z^3;
fac2 = (2 (-3 I + Sqrt[3] - Sqrt[3] x) y +
3 (2 (3 I + Sqrt[3]) +
2 (-3 I + Sqrt[3]) x - (-I + Sqrt[3]) x^2) y^2 z + (-8 I -
4 (-I + Sqrt[3]) x + (I + Sqrt[3]) x^2) z^3);
Use "random" values for two variables and factor the resulting univariate over the algebraic extension.
vals = {x -> 13/4, y -> 7/3}; fax =
Select[FactorList[algpoly /. vals,
Extension -> {I, Sqrt[3]}], ! NumberQ[#[[1]]] &];
newfac = fax[[1, 1]]
(* Out[601]=
18060 I - 25284 Sqrt[3] + (-64043 I + 82075 Sqrt[3]) z + 33282 I z^3 *)
Set up an initial polynomial set consisting of the factor and the polynomials representing the evaluations at two variables. We do this to finite precision. Possibly we could replace algebraics by new variables and polynomial defining relations and just maybe that would be fast enough. I did not try it so I don't know offhand but I'm not so optimistic.
prec = 600;
gby[0] = Flatten[N[{newfac, {x, y} - ({x, y} /. vals)}, prec]];
Now "lift" the factorization in powers of x, leaving y at its given value. I use the bound from the paper which should be fine for this step.
n = Ceiling[Log[2, Exponent[newfac, z] + 1]]
(* Out[590]= 2 *)
Do[gby[j] =
GroebnerBasis[
Flatten[{algpoly, gby[j - 1]^2, y - 7/3}], {x, y, z},
MonomialOrder -> DegreeReverseLexicographic,
CoefficientDomain -> InexactNumbers];
, {j, n}];
Now initialize a new basis to the last one we computed and lift again to get corrected powers in y. Use of the bivariate lifting bound is suspect, at best. And indeed the degree-minimal polynomial in the result is not a correct factor (as would hold in the bivariate case). But the next one is.
gb[0] = gby[n];
Do[gb[j] =
GroebnerBasis[Flatten[{algpoly, gb[j - 1]^2}], {x, y, z},
MonomialOrder -> DegreeReverseLexicographic,
CoefficientDomain -> InexactNumbers];
, {j, n}];
Let's check that second element (the first is a quartic in y alone, specifically, the expanded (y-7/3)^4).
gb[2][[2]] // N
(* Out[606]= (1. + 0.57735026919 I) y - (0.5 -
0.288675134595 I) x y + (0. + 3.46410161514 I) y^2 z - (3. +
1.73205080757 I) x y^2 z + (1. +
0. I) x^2 y^2 z - (0.666666666667 +
1.15470053838 I) z^3 + (1.33333333333 +
0. I) x z^3 - (0.166666666667 - 0.288675134595 I) x^2 z^3 *)
To see that this is a multiple of one of the known factors, can use polynomial division.
PolynomialReduce[gb[2][[2]], fac1, {x, y, z},
CoefficientDomain -> InexactNumbers]
(* Out[563]= \
{{-0.14433756729740644112728719512548936391190043781753171900465058162\
09944180757333364234288488964373813063021784512838919191641620912499 \
+ 0.083333333333333333333333333333333333333333333333333333333333333333\
33333333333333333333333333333333333333333333333333333333333333333333 \
I}, 0} *)
The quotient, in the first part, is a numerical constant. The remainder is 0. In the real world one would instead check that it divides the original polynomial (to good numeric approximation), as below.
N[
PolynomialReduce[algpoly, gb[2][[2]], {x, y, z},
CoefficientDomain -> InexactNumbers]]
(* Out[611]= {{1. ((-36. + 20.7846096908 I) y + (18. +
10.3923048454 I) x y - (0. + 124.707658145 I) y^2 z - (108. -
62.3538290725 I) x y^2 z + (36. + 0. I) x^2 y^2 z - (24. -
41.5692193817 I) z^3 + (48. + 0. I) x z^3 - (6. +
10.3923048454 I) x^2 z^3)}, 0. + 0. I} *)
Last, one needs to recover the exact numerical coefficients. That's fairly easy in this case, using RootApproximant. In general this step could require more precision than that which was required to get an approximate factorization, so this is another soft spot in the overall method.
gb[2][[2]] /. coeff : (_Real | _Complex) :> RootApproximant[coeff]
(* Out[618]=
1/3 (3 + I Sqrt[3]) y + 1/6 (-3 + I Sqrt[3]) x y +
2 I Sqrt[3] y^2 z + (-3 - I Sqrt[3]) x y^2 z + x^2 y^2 z +
1/3 (-2 - 2 I Sqrt[3]) z^3 + (4 x z^3)/3 +
1/6 (-1 + I Sqrt[3]) x^2 z^3 *)
As noted above, making this effective would require figuring out the appropriate lift bounds for the third and later variables, getting the lifting code to be more efficient (perhaps by using something other than GroebnerBasis).
