[#118] Add test to divModE

The tests for `divModE` fail, so the reason for the function's
incorrect behavior is under investigation.

`divModE'` is a similar division algorith used for
`GaussianIntegers` that works, if no other way to fix the
current process used by `divModE`, `divModE'` will be
used instead.

Math/NumberTheory/EisensteinIntegers.hs

@@ -19,13 +19,16 @@ module Math.NumberTheory.EisensteinIntegers
   , conjugate
   , divE
   , divModE
+  , divModE'
+  , eisensteinToComplex
   , modE
   , norm
   ) where
 
-import qualified Data.Complex as C
-import Data.Ratio             ((%))
-import GHC.Generics           (Generic)
+import qualified Data.Complex                     as C
+import Data.Ratio                                 ((%))
+import Debug.Trace                                (trace)
+import GHC.Generics                               (Generic)
 
 infix 6 :+
 
@@ -79,6 +82,18 @@ instance Num EisensteinInteger where
 realEisen :: EisensteinInteger -> Rational
 realEisen (a :+ b) = fromInteger a - (b % 2)
 
+divModE'
+    :: 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)
+
+
 -- | 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>
@@ -87,51 +102,64 @@ divModE
     :: EisensteinInteger
     -> EisensteinInteger
     -> (EisensteinInteger, EisensteinInteger)
-divModE alfa beta | norm rho1 < norm rho2 = (gamma1, rho1)
-                  | norm rho1 > norm rho2 = (gamma2, rho2)
+divModE alfa beta | norm rho1 < norm rho2 = res1
+                  | norm rho1 > norm rho2 = res2
                   | norm rho1 == norm rho2 && 
-                    realEisen gamma1 > realEisen gamma2 = (gamma1, rho1)
-                  | otherwise = (gamma2, 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 ω
+    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
+    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
+    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)
+    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 Alrgorithm from
+    -- 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')
+    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')
 
 
 -- | Eisenstein integer division, truncating toward negative infinity.

test-suite/Math/NumberTheory/EisensteinIntegersTests.hs

@@ -38,9 +38,21 @@ 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.
+divModProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
+divModProperty1 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
+-- 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)
+
 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
   ]