Testing for primality in quadratic rings?

Question

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.

Answer 1

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}}