## Testing for primality in quadratic rings?

8

1

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 – 2012-11-01T14:49:59.480

2You could use Reduce[] to solve a Pell equation, for starters... – J. M.'s ennui – 2012-11-01T17:45:51.927

Ugh. I suppose a more specialized algorithm for solving a Pell equation is in order... – J. M.'s ennui – 2012-11-01T18:15:11.380

1

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 – 2012-11-01T21:10:45.747

Now that @Daniel brought up Cornacchia, I seem to remember that Stan Wagon already implemented that method in Mathematica... – J. M.'s ennui – 2012-11-02T01:31:14.273

@J.M.: Heh, he just emailed me (on an unrelated topic) yesterday... – Charles – 2012-11-02T01:42:24.467

Yeah, I'd imagine Stan has done code for that. Actually I read up on it in a joint expository paper by Buhler and Wagon. – Daniel Lichtblau – 2012-11-02T16:30:22.487

1Are 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 – 2013-02-15T14:42:35.103

@DavidSpeyer: Actually I'd like to test for irreducibility. – Charles – 2013-02-15T23:26:49.723

7

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:= foo = 5336^2 + 151*24^2
Out= 28559872

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


3Thank 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 – 2013-02-18T20:09:18.160