# Testing for primality in quadratic rings?

Testing for primality in $\mathbb{Z}[\sqrt{-1}]$ in Mathematica is easy:

PrimeQ[n, GaussianIntegers -> True]


But how can I test for primality in, say, $\mathbb{Z}[\sqrt{-7}]$? I'm interested in a large number of quadratic rings, not just that one.

I'll try to be more specific in case it helps. Of course much of the work is exploratory but some things remain roughly constant. $n$, the number to be tested for primality, is roughly 10 digits long. $d$ is almost always from 13 to 50, and in any case not more than a thousand.

I tried to express the problem $$p=(a+b\sqrt{-d})(c+e\sqrt{-d})$$ as a Diophantine equation to use Mathematica's general solving tools, but they seem to not be up to the task. Using

f[d_, p_] := FindInstance[a c - b d e == p && a e + b c == 0 && a^2 + b^2 > 1, {a, b, c, e}, Integers]


I get

FindInstance::nsmet: The methods available to FindInstance are insufficient to find the requested instances or prove they do not exist

If I use

f[d_, p_] := Resolve[Exists[{a, b, c, e}, a c - b d e == p && a e + b c == 0 && a^2 + b^2 > 1], Integers]


the program runs forever (or at least more than half an hour; the analogous calculation on PARI/GP takes about 2 milliseconds). Finally

f[d_, p_] := FullSimplify[Reduce[{a c - b d e == p, a e + b c == 0, a^2 + b^2 > 1}, {a, b, c, e}, Integers]]


returns the result unsimplified.

• Or, more generally: how can I factor a number over $\mathbb{Z}[\sqrt{-d}]$? – Charles Nov 1 '12 at 14:49
• You could use Reduce[] to solve a Pell equation, for starters... – J. M. is away Nov 1 '12 at 17:45
• Ugh. I suppose a more specialized algorithm for solving a Pell equation is in order... – J. M. is away Nov 1 '12 at 18:15
• I second that Pell equation approach. For d=7 (so you work over Z[sqrt(-7)]) and n=4234567 you could do Reduce[a^2 + 7*b^2 == n, {a, b}, Integers]. Or code the Smith-Cornacchia algorithm – Daniel Lichtblau Nov 1 '12 at 21:10
• Are you testing for primality or irreducibility? That is to say, if you give Mathematica the number $2$ in $\mathbb{Z}[\sqrt{-5}]$, do you want it to output Prime (because $2=a^2+5b^2$ has no solution) or NotPrime (because the ideal $(2)$ in the ring $\mathbb{Z}[\sqrt{-5}]$ is not prime). – David E Speyer Feb 15 '13 at 14:42

I had a clever idea for how to do this using LatticeReduce[], but I decided to code up the Smith-Cornacchia algorithm first to bench mark against, and it was effectively instant for inputs in your range. Here is a sloppy implementation. In particular, I am embarrassed by applying Divisors[] to something which is computed as a product. However, the result is fast enough that I am not going to bother to improve it.

(* The square roots of d modulo m. Returns {} if d is not square *)
roots[d_, m_] :=
With[{rr = Solve[{x^2 == d, Modulus == m}, x, Mode -> Modular]},
If[rr == {}, {}, x /. rr]]

cornGCD[m_, b_] :=
Module[{c = m, d = b, e},
(While[d^2 > m, (e = Mod[c, d]; c = d; d = e)]; d)]

(* Solves a^2+d b^2 == m with GCD[a,b] == 1 . *)
sfCorn[d_, m_] :=
Select[
Map[
{cornGCD[m, #], Sqrt[(m - cornGCD[m, #]^2) /d]} &,
roots[-d, m]],
IntegerQ[Last[#]]&]

(* Solves a^2+d b^2 == m. Works by calling sfCorn on the equation
aa^2+d bb^2 = m/s^2 for every s such that s^2 divides m. *)
Corn[d_, m_] :=
Union[Flatten[Map[
#*sfCorn[d, m/#^2]&,
Divisors[
Apply[Times,
Map[First[#]^Floor[Last[#]/2]&, FactorInteger[m]]]]], 1]]


Sample usage:

In[74]:= foo = 5336^2 + 151*24^2
Out[74]= 28559872

In[75]:=Corn[151, foo]
Out[75]={{4026, 286}, {5034, 146}, {5336, 24}}

• Thank you for this code. From V6 onwards you could use PowerModList[d,1/2,m]instead of your function roots[d,m]. The only difference occurs when the only solution is 0. Then PowerModList returns $\{0\}$, and roots returns $\{0,0\}$. – KennyColnago Feb 18 '13 at 20:09