From 673db261e66f373566288dd2d1cb3de388748e53 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alexandre=20Bald=C3=A9?= <alexandrer_b@outlook.com>
Date: Thu, 16 Aug 2018 19:14:28 +0000
Subject: [PATCH] [#118] Add `gcd` for Eisenstein integers

---
 Math/NumberTheory/EisensteinIntegers.hs       | 14 +++++++-
 .../NumberTheory/EisensteinIntegersTests.hs   | 32 +++++++++++++++++++
 2 files changed, 45 insertions(+), 1 deletion(-)

diff --git a/Math/NumberTheory/EisensteinIntegers.hs b/Math/NumberTheory/EisensteinIntegers.hs
index 0a3497055..a38552e66 100644
--- a/Math/NumberTheory/EisensteinIntegers.hs
+++ b/Math/NumberTheory/EisensteinIntegers.hs
@@ -27,6 +27,8 @@ module Math.NumberTheory.EisensteinIntegers
   , quotE
   , remE
 
+  , gcdE
+
   -- * Primality functions
   , isPrime
   ) where
@@ -144,4 +146,14 @@ isPrime e@(a :+ b)
     | a == 0 && b `mod` 3 == 2           = Testing.isPrime b
     | nE `mod` 3 == 0 || nE `mod` 3 == 1 = Testing.isPrime nE
     | otherwise = False
-  where nE = norm e
\ No newline at end of file
+  where nE = norm e
+
+-- | Compute the GCD of two Eisenstein integers. The result is always
+-- in the first sextant.
+gcdE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
+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))
\ No newline at end of file
diff --git a/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs b/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs
index 359f31df3..f087e94f0 100644
--- a/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs
+++ b/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs
@@ -17,6 +17,8 @@ module Math.NumberTheory.EisensteinIntegersTests
 import qualified Math.NumberTheory.EisensteinIntegers as E
 
 import Test.Tasty                                     (TestTree, testGroup)
+import Test.Tasty.HUnit                               (Assertion, assertEqual,
+                                                      testCase)
 
 import Math.NumberTheory.TestUtils                    (testSmallAndQuick)
 
@@ -67,6 +69,27 @@ quotRemProperty1 x y = (y == 0) || q == q' && r == r'
 quotRemProperty2 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
 quotRemProperty2 x y = (y == 0) || (x `E.quotE` y) * y + (x `E.remE` y) == x
 
+-- | Verify that @gcd z2 z2@ always divides @z1@ and @z2@.
+gcdEProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
+gcdEProperty1 z1 z2
+  = z1 == 0 && z2 == 0
+  || z1 `E.remE` z == 0 && z2 `E.remE` z == 0 && z == abs z
+  where
+    z = E.gcdE z1 z2
+
+-- | Verify that a common divisor of @z1, z2@ is a always divisor of @gcdE z1 z2@.
+gcdEProperty2 :: E.EisensteinInteger -> E.EisensteinInteger -> E.EisensteinInteger -> Bool
+gcdEProperty2 z z1 z2
+  = z == 0
+  || (E.gcdE z1' z2') `E.remE` z == 0
+  where
+    z1' = z * z1
+    z2' = z * z2
+
+-- | A special case that tests rounding/truncating in GCD.
+gcdESpecialCase1 :: Assertion
+gcdESpecialCase1 = assertEqual "gcdE" 1 $ E.gcdE (12 E.:+ 23) (23 E.:+ 34)
+
 testSuite :: TestTree
 testSuite = testGroup "EisensteinIntegers" $
   [ testSmallAndQuick "forall z . z == signum z * abs z" signumAbsProperty
@@ -81,4 +104,13 @@ testSuite = testGroup "EisensteinIntegers" $
     , testSmallAndQuick "quotE and remE work properly" quotRemProperty1
     , testSmallAndQuick "quotRemE works properly" quotRemProperty2
     ]
+
+  , testGroup "g.c.d."
+    [ testSmallAndQuick "The g.c.d. of two Eisenstein integers divides them"
+                        gcdEProperty1
+    , testSmallAndQuick "A common divisor of two Eisenstein integers always\
+                        \ divides the g.c.d. of those two integers"
+                        gcdEProperty2
+    , testCase          "g.c.d. (12 :+ 23) (23 :+ 34)" gcdESpecialCase1
+    ]
   ]