[Bodigrim#118] Add findPrime function

In addition, the `isPrime` function has been fixed as it
contained two errors.

- Firstly, it checked if the norm of
an Eisenstein integer was congruent to 0 modulo 3, which
is not necessary because it is supposed to be a prime
integer.

- Secondly, it did not account for the possibility of an
Eisenstein integer being the product of an integer prime
and an Eisenstein unit. For example, `isPrime 11` would be
`True`, but `isPrime $ 11 * (1 + ω)` would be `False`, which
is incorrect.

Math/NumberTheory/EisensteinIntegers.hs

@@ -30,12 +30,17 @@ module Math.NumberTheory.EisensteinIntegers
   , gcdE
 
   -- * Primality functions
+  , findPrime
   , isPrime
   ) where
 
 import GHC.Generics                               (Generic)
 
+import qualified Math.NumberTheory.Moduli         as Moduli
+import Math.NumberTheory.Moduli.Sqrt              (FieldCharacteristic(..))
+import Math.NumberTheory.Primes.Types             (PrimeNat(..))
 import qualified Math.NumberTheory.Primes.Testing as Testing
+import Math.NumberTheory.Utils.FromIntegral       (integerToNatural)
 
 infix 6 :+
 
@@ -141,12 +146,13 @@ norm (a :+ b) = a*a - a * b + b*b
 -- page 12.
 isPrime :: EisensteinInteger -> Bool
 isPrime e@(a :+ b)
-    | e == 0                             = False
-    | b == 0 && a `mod` 3 == 2           = Testing.isPrime a
-    | a == 0 && b `mod` 3 == 2           = Testing.isPrime b
-    | nE `mod` 3 == 0 || nE `mod` 3 == 1 = Testing.isPrime nE
+    | e == 0                     = False
+    | b' == 0 && a' `mod` 3 == 2 = Testing.isPrime a
+    | a' == 0 && b' `mod` 3 == 2 = Testing.isPrime b
+    | nE `mod` 3 == 1            = Testing.isPrime nE
     | otherwise = False
-  where nE = norm e
+  where nE       = norm e
+        a' :+ b' = abs e
 
 -- | Compute the GCD of two Eisenstein integers. The result is always
 -- in the first sextant.
@@ -156,4 +162,15 @@ gcdE g h = gcdE' (abs g) (abs h)
 gcdE' :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
 gcdE' g h
     | h == 0    = g -- done recursing
-    | otherwise = gcdE' h (abs (g `modE` h))
+    | otherwise = gcdE' h (abs (g `modE` h))
+
+-- | Find an Eisenstein integer whose norm is the given prime number
+-- in the form @3k + 1@ using a modification of the
+-- <http://www.ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/S0025-5718-1972-0314745-6.pdf Hermite-Serret algorithm>.
+findPrime :: Integer -> EisensteinInteger
+findPrime p = case Moduli.sqrtModMaybe (9*k*k - 1) (FieldCharacteristic (PrimeNat . integerToNatural $ p) 1) of
+    Nothing   -> error "findPrime: argument must be prime p = 6k + 1"
+    Just zed  -> gcdE (p :+ 0) ((zed - 3 * k) :+ 1)
+    where
+        k :: Integer
+        k = p `div` 6

test-suite/Math/NumberTheory/EisensteinIntegersTests.hs

@@ -15,12 +15,13 @@ module Math.NumberTheory.EisensteinIntegersTests
   ) where
 
 import qualified Math.NumberTheory.EisensteinIntegers as E
-
+import Math.NumberTheory.Primes                       (primes)
 import Test.Tasty                                     (TestTree, testGroup)
 import Test.Tasty.HUnit                               (Assertion, assertEqual,
                                                       testCase)
 
-import Math.NumberTheory.TestUtils                    (testSmallAndQuick)
+import Math.NumberTheory.TestUtils                    (Positive (..),
+                                                       testSmallAndQuick)
 
 -- | Check that @signum@ and @abs@ satisfy @z == signum z * abs z@, where @z@ is
 -- an @EisensteinInteger@.
@@ -90,6 +91,13 @@ gcdEProperty2 z z1 z2
 gcdESpecialCase1 :: Assertion
 gcdESpecialCase1 = assertEqual "gcdE" 1 $ E.gcdE (12 E.:+ 23) (23 E.:+ 34)
 
+findPrimesProperty1 :: Positive Int -> Bool
+findPrimesProperty1 (Positive index) =
+    let -- Only retain primes that are of the form @6k + 1@, for some nonzero natural @k@.
+        prop prime = prime `mod` 6 == 1
+        p = (!! index) $ filter prop $ drop 3 primes
+    in E.isPrime $ E.findPrime p
+
 testSuite :: TestTree
 testSuite = testGroup "EisensteinIntegers" $
   [ testSmallAndQuick "forall z . z == signum z * abs z" signumAbsProperty
@@ -113,4 +121,11 @@ testSuite = testGroup "EisensteinIntegers" $
                         gcdEProperty2
     , testCase          "g.c.d. (12 :+ 23) (23 :+ 34)" gcdESpecialCase1
     ]
+
+  , testGroup "Primality"
+    [ testSmallAndQuick "Eisenstein primes found by the norm search used in\
+                        \ findPrime are really prime"
+                        findPrimesProperty1
+
+    ]
   ]