FactorInteger over UFDs

Question

How can I factor 'integers' over other quadratic number fields (not just gaussian integers).

For instance, how could I factor $7 = (3 + ω)(2 − ω)$ over Eisenstein integers ($ω = \frac{-1+ I \sqrt{3}}{2}$)

What I want to achieve is something like Factor[n, Extension->Sqrt[n]]

Apologies for any misused terminology :)

Answer 1

I don't think there is anything built in, but here is a naive implementation -- trial division without any redundancy checks.

Note this only will factor over Euclidean domains of the form $\mathbb{Z}[\alpha]$.

(* FactorIntegerED[{a, b}, ω, Nrm] factors a + b ω over Z[ω] with norm Nrm  *)

FactorIntegerED[{m_, n_}, ω_, norm_] := Module[{a = m, b = n, div, res, cands, x, y, unit},
    div = Select[Divisors[norm[m, n]], PrimeQ];

    res = Reap[
        Do[

            cands = Reverse[{x, y} /. Solve[norm[x, y] == div[[i]], {x, y}, Integers]];

            Do[
                While[divides[{a, b}, cands[[j]], ω],
                    Sow[cands[[j]].{1, ω}];
                    {a, b} = divide[{a, b}, cands[[j]], ω]
                ],
                {j, Length[cands]}
            ],

            {i, Length[div]}
        ]
    ][[-1]];

    If[Length[res] == 0 || Length[First[res]] == 1, 
        {{m + n ω, 1}},

        res = Tally[First[res]];
        unit = Solve[(Times@@Power@@@res)(q + w ω) == m + n ω && {q, w} ∈ Integers, {q, w}];
        Prepend[res, {First[{q, w} /. unit].{1, ω}, 1}]
    ]
]

divide[{a_, b_}, {x_, y_}, ω_] := divide[{a, b}, {x, y}, ω] = 
  With[{sol = Solve[q + w ω == (a + b ω)/(x + y ω) && {q, w} ∈ Integers, {q, w}]},
    If[!MatchQ[sol, {{_, _}}], 
        $Failed,
        First[{q, w} /. sol]
    ]
]

divides[e__] := ListQ[divide[e]]

This is very slow but gets the job done (for you example at least, I haven't tested if out for anything else).

FactorIntegerED[{7, 0}, (-1)^(2/3), Function[{a, b}, a^2-a b+b^2]]//AbsoluteTiming

(* {1.840636, {{-(-1)^(2/3), 1}, {3+2 (-1)^(2/3), 1}, {3+(-1)^(2/3), 1}}} *)

Edit: We can use this slow code to see which prime integers are still prime over $\mathbb{Z}[\omega]$.

myCenterDot[e_] := e
myCenterDot[l__] := CenterDot[l]

mySuperScript[p_, e_] := HoldForm[p^e]

ω = (-1)^(2/3);
Nrm = Function[{a, b}, a^2 - a b + b^2];
P = Prime[Range[26]];

Grid[Transpose@Partition[
    # == myCenterDot @@ mySuperScript @@@ FactorIntegerUFD[{#, 0}, ω, Nrm] & /@ P, 
  13], Spacings -> {4, .75}, Alignment -> Left] /. ω -> "ω" // TraditionalForm

Luckily Wikipedia agrees with this output.

Answer 2

I don't exactly know what I'm doing, but...here is an outline using Smith's method to recast a number as a sum of squares.

the overall approach is this.

(1) Factor into ordinary primes.

(2) For each ordinary prime p, factor over Z[sqrt(-3)].

(3) Now replace sqrt(-3) with 2ω+1.

Clearly (1) and (3) are straightforward in Mathematica. For (2) one can do as follows.

(2a) Find a square root of -3 modulo p. Call it s

(2b) Find the first element in the Euclidean Remainder Sequence of p and s that is less than sqrt(p). Call it a. Then we'll have p = a^2+3*b^2 for some integer b which is now easy to compute.

(2c) We now factor p as (a - sqrt(-3)*b)*(a + sqrt(-3)*b)

Here is an example. I'll use the prime 127.

p = 127;
rt = PowerMod[-3, 1/2, p]

(* Out[150]= 39 *)

(* Now find the first element in the ERS less than sqrt(127). No need
to work hard, there's an internal function that can do this. *)

a = Internal`HGCD[p, rt][[-1, -1]]

(* Out[151]= 10 *)

(* Now compute *b* *)

b = Sqrt[(p - a^2)/3]

(* Out[155]= 3 *)

(* Do the replacement to factor in terms of ω *)

fax = Expand[{a - sqrtm3*b, a + sqrtm3*b} /. sqrtm3 -> 2 ω+ 1]

{7 - 6 ω, 13 + 6 ω}

(* Verify it is correct *)

Expand[Times @@ fax /. ω -> (-1 + I*Sqrt[3])/2]

(* Out[157]= 127 *)

--- edit ---

This is now built in as a resource function.

https://resources.wolframcloud.com/FunctionRepository/resources/EisensteinFactorInteger/

There is a more complete version, one that factors explicit Eisenstein integers, in a function called EisensteinIntegers. I'm still working on the internals of that function so I'll omit the link for now (tho' it's an obvious modification of the one above).

Also I placed HGCD there.

https://resources.wolframcloud.com/FunctionRepository/resources/HalfGCD/

It includes a reasonably efficient variation for GaussianIntegers. We already have internal code for that but I didn't want to mess with its interface to expose it, and I do not know for a certainty that it meets the "half-gcd" guarantee on norms.

--- end edit ---