[#118] Fix Eisenstein integer's division algorithm

Math/NumberTheory/EisensteinIntegers.hs

@@ -17,19 +17,20 @@ module Math.NumberTheory.EisensteinIntegers
   ( EisensteinInteger(..)
   , ω
   , conjugate
+  , norm
+
+  -- * Division and remainder functions
   , divE
   , divModE
-  , divModE'
-  , eisensteinToComplex
   , modE
-  , norm
+  , quotRemE
+  , quotE
+  , remE
   ) where
 
-import qualified Data.Complex                     as C
-import Data.Ratio                                 ((%))
-import Debug.Trace                                (trace)
 import GHC.Generics                               (Generic)
 
+
 infix 6 :+
 
 -- |An Eisenstein integer is a + bω, where a and b are both integers.
@@ -43,12 +44,6 @@ data EisensteinInteger = (:+) { real :: !Integer, imag :: !Integer }
 ω :: EisensteinInteger
 ω = 0 :+ 1
 
-eisensteinToComplex :: forall a . RealFloat a => EisensteinInteger -> C.Complex a
-eisensteinToComplex (a :+ b) = (a' - b' / 2) C.:+ ((b' * sqrt 3) / 2)
-  where
-    a' = fromInteger a
-    b' = fromInteger b
-
 instance Show EisensteinInteger where
     show (a :+ b)
         | b == 0     = show a
@@ -76,91 +71,29 @@ instance Num EisensteinInteger where
         | a == 0 && b == 0 = z               -- hole at origin
         | otherwise        = z `divE` abs z
 
--- | Function to return the real part of an @EisensteinInteger@ @α + βω@ when
--- written in the form @a + bι@, where @α, β@ are @Integer@s and @a, b@ are
--- real numbers.
-realEisen :: EisensteinInteger -> Rational
-realEisen (a :+ b) = fromInteger a - (b % 2)
-
-divModE'
+-- |Simultaneous 'quot' and 'rem'.
+quotRemE
     :: EisensteinInteger
     -> EisensteinInteger
     -> (EisensteinInteger, EisensteinInteger)
-divModE' g h =
-    let nr :+ ni = g * conjugate h
-        denom = norm h
-        q = div nr denom :+ div ni denom
-        p = h * q
-    in (q, g - p)
+quotRemE = divHelper quot
 
+-- |Eisenstein integer division, truncating toward zero.
+quotE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
+n `quotE` d = q where (q,_) = quotRemE n d
+
+-- |Eisenstein integer remainder, satisfying
+--
+-- > (x `quotE` y)*y + (x `remE` y) == x
+remE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
+n `remE`  d = r where (_,r) = quotRemE n d
 
 -- | Simultaneous 'div' and 'mod' of Eisenstein integers.
--- The algorithm used here was derived from
--- <http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.892.1219&rep=rep1&type=pdf NTRU over the Eisenstein Integers>
--- by K. Jarvis, Theorem 4.1.1.
 divModE
     :: EisensteinInteger
     -> EisensteinInteger
     -> (EisensteinInteger, EisensteinInteger)
-divModE alfa beta | norm rho1 < norm rho2 = res1
-                  | norm rho1 > norm rho2 = res2
-                  | norm rho1 == norm rho2 && 
-                    realEisen gamma1 > realEisen gamma2 = res1
-                  | otherwise = res2
-  where
-    res1 = trace show' (gamma1, rho1)
-    res2 = trace show' (gamma2, rho2)
-    show' = "norm ρ1: " ++ show (norm rho1) ++ "\n" ++
-            "norm ρ2: " ++ show (norm rho2) ++ "\n" ++
-            "ρ1': " ++ show rho1' ++ "\n" ++
-            "ρ2': " ++ show rho2' ++ "\n" ++
-            "ρ1' * β': " ++ show (beta' * rho1') ++ "\n" ++
-            "ρ2' * β': " ++ show (beta' * rho2') ++ "\n" ++
-            "(γ1,ρ1): " ++ show (gamma1, rho1) ++ "\n" ++
-            "(γ2,ρ2): " ++ show (gamma2, rho2) ++ "\n" ++
-            "Re(γ1): " ++ show (realEisen gamma1) ++ "\n" ++
-            "Re(γ2): " ++ show (realEisen gamma2) ++ "\n"
-
-    -- @sqrt 3@ is used many times throughout the function, as such it's
-    -- calculated here once.
-    sqrt3                :: Double
-    sqrt3                = sqrt 3
-
-    -- First step of assignments performed in the Division Algorithm from
-    -- Theorem 4.1.1.
-    alfa'              = eisensteinToComplex alfa
-    beta'              = eisensteinToComplex beta
-    a1 C.:+ b1         = alfa' / beta'
-    a2 C.:+ b2         = (alfa' / beta') - eisensteinToComplex ω
-
-    -- @round@s a @Double@ and returns that as another @Double@.
-    dToIToD            :: Double -> Double
-    dToIToD            = (fromIntegral :: Integer -> Double) . round
-
-    -- Second step of assignments performed in the Division Algorithm from
-    -- Theorem 4.1.1.
-    properFraction'    :: Double -> (Integer, Double)
-    properFraction'    = properFraction
-    rho'               :: Double -> Double -> C.Complex Double
-    rho' a' b'         = (snd $ properFraction' a') C.:+ (b' - sqrt3 * dToIToD (b' / sqrt3))
-    rho1'              = rho' a1 b1
-    rho2'              = rho' a2 b2
-
-    -- Converts a complex number in the form @a + bι@, where @a, b@ are
-    -- @Double@s, into an @EisensteinInteger@ in the form @α + βω@, where
-    -- @α, β@ are @Integer@s.
-    toEisen            :: C.Complex Double -> EisensteinInteger
-    toEisen (x C.:+ y) = round (x + y / sqrt3) :+ round (y * 2 / sqrt3)
-
-    -- Third step of assignments performed in the Division Algorithm from
-    -- Theorem 4.1.1.
-    rho1               = toEisen $ beta' * rho1'
-    b1sqrt3'           = round $ b1 / sqrt3
-    gamma1             = ((round a1) + b1sqrt3') :+ (2 * b1sqrt3')
-    rho2               = toEisen $ beta' * rho2'
-    b2sqrt3'           = round $ b2 / sqrt3
-    gamma2             = ((round a2) + b2sqrt3') :+ (1 + 2 * b2sqrt3')
-
+divModE = divHelper div
 
 -- | Eisenstein integer division, truncating toward negative infinity.
 divE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
@@ -172,6 +105,22 @@ n `divE` d = q where (q,_) = divModE n d
 modE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
 n `modE` d = r where (_,r) = divModE n d
 
+-- | Function that does most of the underlying work for @divMod@ and
+-- @quotRem@, apart from choosing the specific integer division algorithm.
+-- This is instead done by the calling function (either @divMod@ which uses
+-- @div@, or @quotRem@, which uses @quot@.)
+divHelper
+    :: (Integer -> Integer -> Integer)
+    -> EisensteinInteger
+    -> EisensteinInteger
+    -> (EisensteinInteger, EisensteinInteger)
+divHelper divide g h =
+    let nr :+ ni = g * conjugate h
+        denom = norm h
+        q = divide nr denom :+ divide ni denom
+        p = h * q
+    in (q, g - p)
+
 -- | Conjugate a Eisenstein integer.
 conjugate :: EisensteinInteger -> EisensteinInteger
 conjugate (a :+ b) = (a - b) :+ (-b)

test-suite/Math/NumberTheory/EisensteinIntegersTests.hs

@@ -38,21 +38,47 @@ absProperty z = isOrigin || (inFirstSextant && isAssociate)
     inFirstSextant = x' > y' && y' >= 0
     isAssociate = z' `elem` map (\e -> z * (1 E.:+ 1) ^ e) [0 .. 5]
 
--- | Verify that `divModE` produces the right quotient and remainder.
+-- | Verify that @divE@ and @modE@ are what `divModE` produces.
 divModProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
-divModProperty1 x y = (y == 0) || (x `E.divE` y) * y + (x `E.modE` y) == x
+divModProperty1 x y = y == 0 || (q == q' && r == r')
+  where
+    (q, r) = E.divModE x y
+    q'     = E.divE x y
+    r'     = E.modE x y
+
+-- | Verify that @divModE` produces the right quotient and remainder.
+divModProperty2 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
+divModProperty2 x y = (y == 0) || (x `E.divE` y) * y + (x `E.modE` y) == x
 
--- | Verify that `divModE` produces a remainder smaller than the divisor with
+-- | Verify that @divModE@ produces a remainder smaller than the divisor with
 -- regards to the Euclidean domain's function.
 modProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
 modProperty1 x y = (y == 0) || (E.norm $ x `E.modE` y) < (E.norm y)
 
+-- | Verify that @quotE@ and @reE@ are what `quotRemE` produces.
+quotRemProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
+quotRemProperty1 x y = (y == 0) || q == q' && r == r'
+  where
+    (q, r) = E.quotRemE x y
+    q'     = E.quotE x y
+    r'     = E.remE x y
+
+-- | Verify that @quotRemE@ produces the right quotient and remainder.
+quotRemProperty2 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
+quotRemProperty2 x y = (y == 0) || (x `E.quotE` y) * y + (x `E.remE` y) == x
+
 testSuite :: TestTree
 testSuite = testGroup "EisensteinIntegers" $
   [ testSmallAndQuick "forall z . z == signum z * abs z" signumAbsProperty
   , testSmallAndQuick "abs z always returns an @EisensteinInteger@ in the\
                       \ first sextant of the complex plane" absProperty
-  , testSmallAndQuick "divModE works properly" divModProperty1
-  , testSmallAndQuick "The remainder's norm is smaller than the divisor's"
-                      modProperty1
+  , testGroup "Division"
+    [ testSmallAndQuick "divE and modE work properly" divModProperty1
+    , testSmallAndQuick "divModE works properly" divModProperty2
+    , testSmallAndQuick "The remainder's norm is smaller than the divisor's"
+                        modProperty1
+
+    , testSmallAndQuick "quotE and remE work properly" quotRemProperty1
+    , testSmallAndQuick "quotRemE works properly" quotRemProperty2
+    ]
   ]