Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
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. Nov 1 '12 at 17:45
Ugh. I suppose a more specialized algorithm for solving a Pell equation is in order... –  J. M. 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 Speyer Feb 15 '13 at 14:42
show 4 more comments

1 Answer 1

up vote 6 down vote accepted

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_] := 
      {cornGCD[m, #], Sqrt[(m - cornGCD[m, #]^2) /d]} &,
      roots[-d, m]],

(* 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_] :=
                  #*sfCorn[d, m/#^2]&,
                      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}}
share|improve this answer
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
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.