[Bodigrim#118] Add division for EisensteinIntegers

The `Num EisensteinInteger` instance now works properly, and
there is now division over the Euclidean domain of
`EisensteinIntegers`.

Math/NumberTheory/EisensteinIntegers.hs

@@ -1,6 +1,8 @@
 -- |
 -- Module:      Math.NumberTheory.EisensteinIntegers
+-- Copyright:   (c) 2018 Alexandre Rodrigues Baldé
 -- Licence:     MIT
+-- Maintainer:  Alexandre Rodrigues Baldé <[email protected]>
 -- Stability:   Provisional
 -- Portability: Non-portable (GHC extensions)
 --
@@ -9,14 +11,21 @@
 --
 
 {-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE RankNTypes    #-}
 
 module Math.NumberTheory.EisensteinIntegers
   ( EisensteinInteger(..)
-  ,  ω
+  , ω
+  , conjugate
   , divE
+  , divModE
+  , modE
+  , norm
   ) where
 
-import GHC.Generics (Generic)
+import qualified Data.Complex as C
+import Data.Ratio             ((%))
+import GHC.Generics           (Generic)
 
 infix 6 :+
 
@@ -31,6 +40,12 @@ 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
@@ -44,25 +59,95 @@ instance Show EisensteinInteger where
 instance Num EisensteinInteger where
     (+) (a :+ b) (c :+ d) = (a + c) :+ (b + d)
     (*) (a :+ b) (c :+ d) = (a * c - b * d) :+ (b * c + a * d - b * d)
-    -- An Eisenstein integer @a :+ b@, with @a, b@ integers, can we written as
-    -- @(2*a - b) / 2 + ((b * sqrt 3) * ι) / 2@, but this number is in the
-    -- same quadrant as @(2*a - b) / 2 + (b * ι) / 2@, and this one in the
-    -- same as @(2*a - b) + b * ι@. Divisions or floating points are not
-    -- necessary.
-    abs z@(x :+ y) = abs' (2*x - y) x
-      where
-        abs' a b
-            | a == 0 && b == 0 =   z             -- origin
-            | a >  0 && b >= 0 =   z             -- first quadrant: (0, inf) x [0, inf)ω
-            | a <= 0 && b >  0 =   b  :+ (-a)    -- second quadrant: (-inf, 0] x (0, inf)ω
-            | a <  0 && b <= 0 = (-a) :+ (-b)    -- third quadrant: (-inf, 0) x (-inf, 0]ω
-            | otherwise        = (-b) :+   a     -- fourth quadrant: [0, inf) x (-inf, 0)ω
+    abs z@(a :+ b)
+        | a == 0 && b == 0                                  = z            -- origin
+        | a > b && b >= 0                                   = z            -- first sextant: 0 ≤ Arg(η) < π/3
+        | b >= a && a > 0                                   = (-ω) * z     -- second sextant: π/3 ≤ Arg(η) < 2π/3
+        | b > 0 && 0 >= a                                   = (-1 - ω) * z -- third sextant: 2π/3 ≤ Arg(η) < π
+        | a < b && b <= 0                                   = - z          -- fourth sextant: -π < Arg(η) < -2π/3 or Arg(η) = π
+        | b <= a && a < 0                                   = ω * z        -- fifth sextant: -2π/3 ≤ Arg(η) < -π/3
+        | otherwise                                         = (1 + ω) * z  -- sixth sextant: -π/3 ≤ Arg(η) < 0
     negate (a :+ b) = (-a) :+ (-b)
     fromInteger n = n :+ 0
     signum z@(a :+ b)
         | 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)
+
+-- | 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 = (gamma1, rho1)
+                  | norm rho1 > norm rho2 = (gamma2, rho2)
+                  | norm rho1 == norm rho2 && 
+                    realEisen gamma1 > realEisen gamma2 = (gamma1, rho1)
+                  | otherwise = (gamma2, rho2)
+  where
+    -- @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 Alrgorithm 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')
+
+
 -- | Eisenstein integer division, truncating toward negative infinity.
 divE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
-divE = undefined
+n `divE` d = q where (q,_) = divModE n d
+
+-- | Eisenstein integer remainder, satisfying
+--
+-- > (x `divE` y)*y + (x `modE` y) == x
+modE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
+n `modE` d = r where (_,r) = divModE n d
+
+-- | Conjugate a Eisenstein integer.
+conjugate :: EisensteinInteger -> EisensteinInteger
+conjugate (a :+ b) = (a - b) :+ (-b)
+
+-- | The square of the magnitude of a Eisenstein integer.
+norm :: EisensteinInteger -> Integer
+norm (a :+ b) = a*a - a * b + b*b

arithmoi.cabal

@@ -134,6 +134,7 @@ test-suite spec
     Math.NumberTheory.ArithmeticFunctions.MertensTests
     Math.NumberTheory.ArithmeticFunctions.SieveBlockTests
     Math.NumberTheory.CurvesTests
+    Math.NumberTheory.EisensteinIntegersTests
     Math.NumberTheory.GaussianIntegersTests
     Math.NumberTheory.GCDTests
     Math.NumberTheory.Moduli.ChineseTests

test-suite/Math/NumberTheory/EisensteinIntegersTests.hs

@@ -0,0 +1,46 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+-- |
+-- Module:      Math.NumberTheory.EisensteinIntegersTests
+-- Copyright:   (c) 2018 Alexandre Rodrigues Baldé
+-- Licence:     MIT
+-- Maintainer:  Alexandre Rodrigues Baldé <alexandrer_b@outlook.
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.EisensteinIntegers
+--
+
+module Math.NumberTheory.EisensteinIntegersTests
+  ( testSuite
+  ) where
+
+import qualified Math.NumberTheory.EisensteinIntegers as E
+
+import Test.Tasty                                     (TestTree, testGroup)
+
+import Math.NumberTheory.TestUtils                    (testSmallAndQuick)
+
+-- | Check that @signum@ and @abs@ satisfy @z == signum z * abs z@, where @z@ is
+-- an @EisensteinInteger@.
+signumAbsProperty :: E.EisensteinInteger -> Bool
+signumAbsProperty z = z == signum z * abs z
+
+-- | Check that @abs@ maps an @EisensteinInteger@ to its associate in first
+-- sextant.
+absProperty :: E.EisensteinInteger -> Bool
+absProperty z = isOrigin || (inFirstSextant && isAssociate)
+  where
+    z'@(x' E.:+ y') = abs z
+    isOrigin = z' == 0 && z == 0
+    -- The First sextant includes the positive real axis, but not the origin
+    -- or the line defined by the linear equation @y = (sqrt 3) * x@ in the
+    -- Cartesian plane.
+    inFirstSextant = x' > y' && y' >= 0
+    isAssociate = z' `elem` map (\e -> z * (1 E.:+ 1) ^ e) [0 .. 5]
+
+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
+  ]

test-suite/Math/NumberTheory/TestUtils.hs

@@ -54,6 +54,7 @@ import Data.Bits
 import GHC.Exts
 import Numeric.Natural
 
+import qualified Math.NumberTheory.EisensteinIntegers as E (EisensteinInteger(..))
 import Math.NumberTheory.GaussianIntegers (GaussianInteger(..))
 import Math.NumberTheory.Moduli.PrimitiveRoot (CyclicGroup(..))
 import qualified Math.NumberTheory.SmoothNumbers as SN
@@ -66,6 +67,13 @@ instance Arbitrary Natural where
   arbitrary = fromInteger <$> (arbitrary `suchThat` (>= 0))
   shrink = map fromInteger . filter (>= 0) . shrink . toInteger
 
+instance Arbitrary E.EisensteinInteger where
+  arbitrary = (E.:+) <$> arbitrary <*> arbitrary
+  shrink (x E.:+ y) = map (x E.:+) (shrink y) ++ map (E.:+ y) (shrink x)
+
+instance Monad m => Serial m E.EisensteinInteger where
+  series = cons2 (E.:+)
+
 instance Arbitrary GaussianInteger where
   arbitrary = (:+) <$> arbitrary <*> arbitrary
   shrink (x :+ y) = map (x :+) (shrink y) ++ map (:+ y) (shrink x)

test-suite/Test.hs

@@ -29,6 +29,8 @@ import qualified Math.NumberTheory.Primes.FactorisationTests as Factorisation
 import qualified Math.NumberTheory.Primes.SieveTests as Sieve
 import qualified Math.NumberTheory.Primes.TestingTests as Testing
 
+import qualified Math.NumberTheory.EisensteinIntegersTests as Eisenstein
+
 import qualified Math.NumberTheory.GaussianIntegersTests as Gaussian
 
 import qualified Math.NumberTheory.ArithmeticFunctionsTests as ArithmeticFunctions
@@ -76,6 +78,7 @@ tests = testGroup "All"
     , Sieve.testSuite
     , Testing.testSuite
     ]
+  , Eisenstein.testSuite
   , Gaussian.testSuite
   , testGroup "ArithmeticFunctions"
     [ ArithmeticFunctions.testSuite