# FactorInteger over UFDs

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 :)

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.

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 = InternalHGCD[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 ---

• How does InternalHGCD work? Sep 15, 2020 at 12:06
• @niloderoock It is an implementation of the "half GCD" for integer pairs. Implementation details here or here Sep 15, 2020 at 13:14