# Why Mathematica can not factorize polynomials over algebraic fields?

I noticed that a factorization over algebraic fields is useless in Mathematica. Here is the example over the field containing I*Sqrt:

Pol=4 (3 I Sqrt (-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 (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)))


This is not a large polynomial and the answer is simple:

Pol=fac1*fac2;

fac1=2 (-3 I - Sqrt + Sqrt x) y -
3 (6 I - 2 Sqrt -
2 (3 I + Sqrt) x + (I + Sqrt) x^2) y^2 z + (8 I -
4 (I + Sqrt) x + (-I + Sqrt) x^2) z^3;

fac2 = (2 (-3 I + Sqrt - Sqrt x) y +
3 (2 (3 I + Sqrt) +
2 (-3 I + Sqrt) x - (-I + Sqrt) x^2) y^2 z + (-8 I -
4 (-I + Sqrt) x + (I + Sqrt) x^2) z^3);


Factor[Pol,Extension->Sqrt[-3]] will never do it. I left it overnight, it consumed 15Gb of memory and couldn't do it.

Maple does such factorizations in seconds. Is there any way do it in Mathematica or its algorithm is simply not powerful enough ?

Pol2 = Resultant[Pol /. √3 -> a, a^2 - 3, a];
facs = Factor[Pol2, Extension -> √3 I]; // AbsoluteTiming
fac = Select[Subsets[facs, {2}], Expand[# == Pol] &][]
Pol == fac // Simplify


{1.3518, Null}

((6 I-2 √3) y+(-3 I-√3) x y-12 √3 y^2 z+(-18 I+6 √3) x y^2 z+6 I x^2 y^2 z+(-4 I+4 √3) z^3+8 I x z^3+(-I-√3) x^2 z^3) ((6 I+2 √3) y+(-3 I+√3) x y-12 √3 y^2 z+(18 I+6 √3) x y^2 z-6 I x^2 y^2 z+(4 I+4 √3) z^3-8 I x z^3+(I-√3) x^2 z^3)

True

• Nicely done. :) – J. M. is in limbo Aug 3 '17 at 3:53
• Unfortunately, this trick works only with very simple polynomials. By taking resultant you basically "square" back your polynomial. I tried to use your trick for something more complicated and got stuck immediately. The only solution which works is copy/paste to Maple, factorize and bring back to Mathematica. – Aus Man Aug 3 '17 at 3:58

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 (-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 (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 + Sqrt x) y -
3 (6 I - 2 Sqrt -
2 (3 I + Sqrt) x + (I + Sqrt) x^2) y^2 z + (8 I -
4 (I + Sqrt) x + (-I + Sqrt) x^2) z^3;

fac2 = (2 (-3 I + Sqrt - Sqrt x) y +
3 (2 (3 I + Sqrt) +
2 (-3 I + Sqrt) x - (-I + Sqrt) x^2) y^2 z + (-8 I -
4 (-I + Sqrt) x + (I + Sqrt) 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}], ! NumberQ[#[]] &];
newfac = fax[[1, 1]]

(* Out=
18060 I - 25284 Sqrt + (-64043 I + 82075 Sqrt) 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 = 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= 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 = 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[] // N

(* Out= (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[], fac1, {x, y, z},
CoefficientDomain -> InexactNumbers]

(* Out= \
{{-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[], {x, y, z},
CoefficientDomain -> InexactNumbers]]

(* Out= {{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[] /. coeff : (_Real | _Complex) :> RootApproximant[coeff]

(* Out=
1/3 (3 + I Sqrt) y + 1/6 (-3 + I Sqrt) x y +
2 I Sqrt y^2 z + (-3 - I Sqrt) x y^2 z + x^2 y^2 z +
1/3 (-2 - 2 I Sqrt) z^3 + (4 x z^3)/3 +
1/6 (-1 + I Sqrt) 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).

It is an input processing issue. Factor uses the extension generated by the algebraic numbers specified through the Extension option and algebraic numbers it finds in the coefficients of the input polynomial. Here it finds Sqrt and I separately, so it uses the degree 4 extension generated by Sqrt[-3], Sqrt, and I. As a workaround, you can replace algebraic number coefficients manually with a new variable that represents Sqrt[-3], use Factor with an undocumented option Modulus->irreducible univariate polynomial, and replace the new variable with Sqrt[-3].

In:= Pol1=Pol/.Complex[0, c_] Sqrt -> c w;

In:= Timing[Factor[Pol1, Modulus->w^2+3]/.w->I Sqrt]//InputForm

Out//InputForm=
{0.204969, (-6*y + (2*I)*Sqrt*y + 3*x*y + I*Sqrt*x*y -
(12*I)*Sqrt*y^2*z - 18*x*y^2*z + (6*I)*Sqrt*x*y^2*z + 6*x^2*y^2*z -
4*z^3 + (4*I)*Sqrt*z^3 + 8*x*z^3 - x^2*z^3 - I*Sqrt*x^2*z^3)*
(6*y + (2*I)*Sqrt*y - 3*x*y + I*Sqrt*x*y + (12*I)*Sqrt*y^2*z -
18*x*y^2*z - (6*I)*Sqrt*x*y^2*z + 6*x^2*y^2*z - 4*z^3 -
(4*I)*Sqrt*z^3 + 8*x*z^3 - x^2*z^3 + I*Sqrt*x^2*z^3)}

• Interesting. Are there any restrictions on the degree or coefficients of the polynomial that can be given to Modulus? – J. M. is in limbo Aug 5 '17 at 13:24
• The coefficients need to be integers and the polynomial needs to be irreducible. – Adam Strzebonski Aug 6 '17 at 21:03

Here is a much cleaner method than what I first posted. It's basically (Barry) Trager's method, and it is a cousin to the code posted by @chyaong.

algpoly =
4 (3 I Sqrt (-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 (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)));


Step 1: We can find the norm by multiplying it with its algebraic conjugate. I'm taking a short cut here in that I just replace sqrt(3) by its negative for that second factor. Really we are working in an extension that contains both sqrt(3) and sqrt(-1) and is of degree 4. It is possible that this more complicated realm is what is giving rise to the slowness (I have put in a question about that).

algpolyNorm = Expand[algpoly*(algpoly /. Sqrt -> -Sqrt)];

facs = Select[FactorList[algpolyNorm], ! NumberQ[#[]] &]

(* Out= {{12 y^2 - 6 x y^2 + 3 x^2 y^2 - 36 y^3 z + 54 x y^3 z -
36 x^2 y^3 z + 9 x^3 y^3 z + 108 y^4 z^2 - 108 x y^4 z^2 +
108 x^2 y^4 z^2 - 54 x^3 y^4 z^2 + 9 x^4 y^4 z^2 + 24 y z^3 -
24 x y z^3 + 12 x^2 y z^3 - 3 x^3 y z^3 - 72 y^2 z^4 +
72 x y^2 z^4 - 66 x^2 y^2 z^4 + 24 x^3 y^2 z^4 - 3 x^4 y^2 z^4 +
16 z^6 - 16 x z^6 + 12 x^2 z^6 - 4 x^3 z^6 + x^4 z^6,
1}, {12 y^2 - 6 x y^2 + 3 x^2 y^2 + 36 y^3 z - 54 x y^3 z +
36 x^2 y^3 z - 9 x^3 y^3 z + 108 y^4 z^2 - 108 x y^4 z^2 +
108 x^2 y^4 z^2 - 54 x^3 y^4 z^2 + 9 x^4 y^4 z^2 - 24 y z^3 +
24 x y z^3 - 12 x^2 y z^3 + 3 x^3 y z^3 - 72 y^2 z^4 +
72 x y^2 z^4 - 66 x^2 y^2 z^4 + 24 x^3 y^2 z^4 - 3 x^4 y^2 z^4 +
16 z^6 - 16 x z^6 + 12 x^2 z^6 - 4 x^3 z^6 + x^4 z^6, 1}} *)
algpolyNorm = Expand[algpoly*(algpoly /. Sqrt -> -Sqrt)];


Step 2: Factor the norm polynomial over the rationals.

facs = Select[FactorList[algpolyNorm], ! NumberQ[#[]] &]

(* Out= {{12 y^2 - 6 x y^2 + 3 x^2 y^2 - 36 y^3 z + 54 x y^3 z -
36 x^2 y^3 z + 9 x^3 y^3 z + 108 y^4 z^2 - 108 x y^4 z^2 +
108 x^2 y^4 z^2 - 54 x^3 y^4 z^2 + 9 x^4 y^4 z^2 + 24 y z^3 -
24 x y z^3 + 12 x^2 y z^3 - 3 x^3 y z^3 - 72 y^2 z^4 +
72 x y^2 z^4 - 66 x^2 y^2 z^4 + 24 x^3 y^2 z^4 - 3 x^4 y^2 z^4 +
16 z^6 - 16 x z^6 + 12 x^2 z^6 - 4 x^3 z^6 + x^4 z^6,
1}, {12 y^2 - 6 x y^2 + 3 x^2 y^2 + 36 y^3 z - 54 x y^3 z +
36 x^2 y^3 z - 9 x^3 y^3 z + 108 y^4 z^2 - 108 x y^4 z^2 +
108 x^2 y^4 z^2 - 54 x^3 y^4 z^2 + 9 x^4 y^4 z^2 - 24 y z^3 +
24 x y z^3 - 12 x^2 y z^3 + 3 x^3 y z^3 - 72 y^2 z^4 +
72 x y^2 z^4 - 66 x^2 y^2 z^4 + 24 x^3 y^2 z^4 - 3 x^4 y^2 z^4 +
16 z^6 - 16 x z^6 + 12 x^2 z^6 - 4 x^3 z^6 + x^4 z^6, 1}} *)


Step 3: Each of these is a product of a factor of the original and the conjugate of that factor. So extract the wanted parts by taking GCDs with the original polynomial.

Map[PolynomialGCD[#, algpoly, Extension -> Sqrt] &,
facs[[All, 1]]]

(* Out= {6 I y + 2 Sqrt y - 3 I x y + Sqrt x y -
12 Sqrt y^2 z + 18 I x y^2 z + 6 Sqrt x y^2 z -
6 I x^2 y^2 z + 4 I z^3 + 4 Sqrt z^3 - 8 I x z^3 + I x^2 z^3 -
Sqrt x^2 z^3,
6 I y - 2 Sqrt y - 3 I x y - Sqrt x y - 12 Sqrt y^2 z -
18 I x y^2 z + 6 Sqrt x y^2 z + 6 I x^2 y^2 z - 4 I z^3 +
4 Sqrt z^3 + 8 I x z^3 - I x^2 z^3 - Sqrt x^2 z^3} *)