Merge pull request #121 from rockbmb/eisenstein-integers

[#118] Add Eisenstein integers module

Math/NumberTheory/GaussianIntegers.hs

@@ -47,6 +47,7 @@ import Math.NumberTheory.Primes.Types (PrimeNat(..))
 import qualified Math.NumberTheory.Primes.Factorisation as Factorisation
 import qualified Math.NumberTheory.Primes.Sieve as Sieve
 import qualified Math.NumberTheory.Primes.Testing as Testing
+import Math.NumberTheory.Utils              (mergeBy)
 import Math.NumberTheory.Utils.FromIntegral (integerToNatural)
 
 infix 6 :+
@@ -76,17 +77,18 @@ instance Show GaussianInteger where
 instance Num GaussianInteger where
     (+) (a :+ b) (c :+ d) = (a + c) :+ (b + d)
     (*) (a :+ b) (c :+ d) = (a * c - b * d) :+ (a * d + b * c)
-    abs z@(a :+ b)
-        | a == 0 && b == 0 =   z             -- origin
-        | a >  0 && b >= 0 =   z             -- first quadrant: (0, inf) x [0, inf)i
-        | a <= 0 && b >  0 =   b  :+ (-a)    -- second quadrant: (-inf, 0] x (0, inf)i
-        | a <  0 && b <= 0 = (-a) :+ (-b)    -- third quadrant: (-inf, 0) x (-inf, 0]i
-        | otherwise        = (-b) :+   a     -- fourth quadrant: [0, inf) x (-inf, 0)i
+    abs = fst . absSignum
     negate (a :+ b) = (-a) :+ (-b)
     fromInteger n = n :+ 0
-    signum z@(a :+ b)
-        | a == 0 && b == 0 = z               -- hole at origin
-        | otherwise        = z `divG` abs z
+    signum = snd . absSignum
+
+absSignum :: GaussianInteger -> (GaussianInteger, GaussianInteger)
+absSignum z@(a :+ b)
+    | a == 0 && b == 0 =   (z, 0)              -- origin
+    | a >  0 && b >= 0 =   (z, 1)              -- first quadrant: (0, inf) x [0, inf)i
+    | a <= 0 && b >  0 =   (b  :+ (-a), ι)     -- second quadrant: (-inf, 0] x (0, inf)i
+    | a <  0 && b <= 0 = ((-a) :+ (-b), -1)    -- third quadrant: (-inf, 0) x (-inf, 0]i
+    | otherwise        = ((-b) :+   a, -ι)     -- fourth quadrant: [0, inf) x (-inf, 0)i
 
 -- |Simultaneous 'quot' and 'rem'.
 quotRemG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)
@@ -148,15 +150,6 @@ primes = (1 :+ 1): mergeBy (comparing norm) l r
         l = [p :+ 0 | p <- leftPrimes]
         r = [g | p <- rightPrimes, let x :+ y = findPrime p, g <- [x :+ y, y :+ x]]
 
-mergeBy :: (a -> a -> Ordering) -> [a] -> [a] -> [a]
-mergeBy cmp = loop
-  where
-    loop [] ys  = ys
-    loop xs []  = xs
-    loop (x:xs) (y:ys)
-      = case cmp x y of
-         GT -> y : loop (x:xs) ys
-         _  -> x : loop xs (y:ys)
 
 -- | Compute the GCD of two Gaussian integers. Result is always
 -- in the first quadrant.

Math/NumberTheory/UniqueFactorisation.hs

@@ -24,6 +24,7 @@ import Data.Coerce
 import qualified Math.NumberTheory.Primes.Factorisation as F (factorise)
 import Math.NumberTheory.Primes.Testing.Probabilistic as T (isPrime)
 import Math.NumberTheory.Primes.Types (Prime, Prm(..), PrimeNat(..))
+import qualified Math.NumberTheory.EisensteinIntegers as E
 import qualified Math.NumberTheory.GaussianIntegers as G
 import Math.NumberTheory.Utils.FromIntegral
 
@@ -76,3 +77,14 @@ instance UniqueFactorisation G.GaussianInteger where
 
   factorise 0 = []
   factorise g = map (coerce *** intToWord) $ G.factorise g
+
+newtype EisensteinPrime = EisensteinPrime { _unEisensteinPrime :: E.EisensteinInteger }
+  deriving (Eq, Show)
+
+type instance Prime E.EisensteinInteger = EisensteinPrime
+
+instance UniqueFactorisation E.EisensteinInteger where
+  unPrime = coerce
+
+  factorise 0 = []
+  factorise e = map (coerce *** intToWord) $ E.factorise e

Math/NumberTheory/Utils.hs

@@ -6,7 +6,7 @@
 -- Stability:   Provisional
 -- Portability: Non-portable (GHC extensions)
 --
--- Some utilities for bit twiddling.
+-- Some utilities, mostly for bit twiddling.
 --
 {-# LANGUAGE CPP, MagicHash, UnboxedTuples, BangPatterns #-}
 {-# OPTIONS_HADDOCK hide #-}
@@ -21,6 +21,8 @@ module Math.NumberTheory.Utils
     , uncheckedShiftR
     , splitOff
     , splitOff#
+
+    , mergeBy
     ) where
 
 #include "MachDeps.h"
@@ -173,3 +175,15 @@ splitOff# p n = go 0# n
       (# q, 0## #) -> go (k +# 1#) q
       _            -> (# k, m #)
 {-# INLINABLE splitOff# #-}
+
+-- | Merges two ordered lists into an ordered list. Checks for neither its
+-- precondition or postcondition.
+mergeBy :: (a -> a -> Ordering) -> [a] -> [a] -> [a]
+mergeBy cmp = loop
+  where
+    loop [] ys  = ys
+    loop xs []  = xs
+    loop (x:xs) (y:ys)
+      = case cmp x y of
+         GT -> y : loop (x:xs) ys
+         _  -> x : loop xs (y:ys)

arithmoi.cabal

@@ -58,6 +58,7 @@ library
     Math.NumberTheory.ArithmeticFunctions.SieveBlock
     Math.NumberTheory.ArithmeticFunctions.Standard
     Math.NumberTheory.Curves.Montgomery
+    Math.NumberTheory.EisensteinIntegers
     Math.NumberTheory.GaussianIntegers
     Math.NumberTheory.GCD
     Math.NumberTheory.GCD.LowLevel
@@ -139,6 +140,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
@@ -191,6 +193,7 @@ benchmark criterion
       semigroups >=0.8
   other-modules:
     Math.NumberTheory.ArithmeticFunctionsBench
+    Math.NumberTheory.EisensteinIntegersBench
     Math.NumberTheory.GaussianIntegersBench
     Math.NumberTheory.GCDBench
     Math.NumberTheory.JacobiBench

benchmark/Bench.hs

@@ -3,6 +3,7 @@ module Main where
 import Gauge.Main
 
 import Math.NumberTheory.ArithmeticFunctionsBench as ArithmeticFunctions
+import Math.NumberTheory.EisensteinIntegersBench as Eisenstein
 import Math.NumberTheory.GaussianIntegersBench as Gaussian
 import Math.NumberTheory.GCDBench as GCD
 import Math.NumberTheory.JacobiBench as Jacobi
@@ -17,6 +18,7 @@ import Math.NumberTheory.SmoothNumbersBench as SmoothNumbers
 main :: IO ()
 main = defaultMain
   [ ArithmeticFunctions.benchSuite
+  , Eisenstein.benchSuite
   , Gaussian.benchSuite
   , GCD.benchSuite
   , Jacobi.benchSuite

benchmark/Math/NumberTheory/EisensteinIntegersBench.hs

@@ -0,0 +1,26 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+{-# OPTIONS_GHC -fno-warn-orphans       #-}
+
+module Math.NumberTheory.EisensteinIntegersBench
+  ( benchSuite
+  ) where
+
+import Control.DeepSeq
+import Gauge.Main
+
+import Math.NumberTheory.ArithmeticFunctions (tau)
+import Math.NumberTheory.EisensteinIntegers
+
+instance NFData EisensteinInteger
+
+benchFindPrime :: Integer -> Benchmark
+benchFindPrime n = bench (show n) $ nf findPrime n
+
+benchTau :: Integer -> Benchmark
+benchTau n = bench (show n) $ nf (\m -> sum [tau (x :+ y) | x <- [1..m], y <- [0..m]] :: Word) n
+
+benchSuite :: Benchmark
+benchSuite = bgroup "Eisenstein"
+  [ bgroup "findPrime" $ map benchFindPrime [1000003, 10000141, 100000039, 1000000021, 10000000033, 100000000003, 1000000000039, 10000000000051]
+  , bgroup "tau" $ map benchTau [10, 20, 40, 80]
+  ]

test-suite/Math/NumberTheory/EisensteinIntegersTests.hs

@@ -0,0 +1,207 @@
+{-# 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 Math.NumberTheory.Primes                       (primes)
+import Test.Tasty                                     (TestTree, testGroup)
+import Test.Tasty.HUnit                               (Assertion, assertEqual,
+                                                      testCase)
+
+import Math.NumberTheory.TestUtils                    (Positive (..),
+                                                       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]
+
+-- | Verify that @divE@ and @modE@ are what `divModE` produces.
+divModProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
+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
+-- 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
+
+-- | Verify that @gcd z1 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)
+
+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
+
+-- | Checks that the @norm@ of the Euclidean domain of Eisenstein integers
+-- is multiplicative i.e.
+-- @forall e1 e2 in Z[ω] . norm(e1 * e2) == norm(e1) * norm(e2)@.
+euclideanDomainProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool
+euclideanDomainProperty1 e1 e2 = E.norm (e1 * e2) == E.norm e1 * E.norm e2
+
+-- | Checks that the numbers produced by @primes@ are actually Eisenstein
+-- primes.
+primesProperty1 :: Positive Int -> Bool
+primesProperty1 (Positive index) = all E.isPrime $ take index $ E.primes
+
+-- | Checks that the infinite list of Eisenstein primes @primes@ is ordered
+-- by the numbers' norm.
+primesProperty2 :: Positive Int -> Bool
+primesProperty2 (Positive index) =
+    let isOrdered :: [E.EisensteinInteger] -> Bool 
+        isOrdered xs = all (\(x,y) -> E.norm x <= E.norm y) . zip xs $ tail xs
+    in isOrdered $ take index E.primes
+
+-- | Checks that the numbers produced by @primes@ are all in the first
+-- sextant.
+primesProperty3 :: Positive Int -> Bool
+primesProperty3 (Positive index) =
+    all (\e -> abs e == e) $ take index $ E.primes
+
+-- | An Eisenstein integer is either zero or associated (i.e. equal up to
+-- multiplication by a unit) to the product of its factors raised to their
+-- respective exponents.
+factoriseProperty1 :: E.EisensteinInteger -> Bool
+factoriseProperty1 g = g == 0 || abs g == abs g'
+  where
+    factors = E.factorise g
+    g' = product $ map (uncurry (^)) factors
+
+-- | Check that there are no factors with exponent @0@ in the factorisation.
+factoriseProperty2 :: E.EisensteinInteger -> Bool
+factoriseProperty2 z = z == 0 || all ((> 0) . snd) (E.factorise z)
+
+-- | Check that no factor produced by @factorise@ is a unit.
+factoriseProperty3 :: E.EisensteinInteger -> Bool
+factoriseProperty3 z = z == 0 || all ((> 1) . E.norm . fst) (E.factorise z)
+
+-- | Check that every prime factor in the factorisation is primary, excluding
+-- @1 - ω@, if it is a factor.
+factoriseProperty4 :: E.EisensteinInteger -> Bool
+factoriseProperty4 z =
+    z == 0 ||
+    (all (\e -> e `E.modE` 3 == 2) $
+     filter (\e -> not $ elem e $ E.associates $ 1 E.:+ (-1)) $
+     map fst $ E.factorise z)
+
+factoriseSpecialCase1 :: Assertion
+factoriseSpecialCase1 = assertEqual "should be equal"
+  [(2 E.:+ 1, 3), (2 E.:+ 3, 1)]
+  (E.factorise (15 E.:+ 12))
+
+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
+  , 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
+    ]
+
+  , 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
+    ]
+  , testSmallAndQuick "The Eisenstein norm function is multiplicative"
+                    euclideanDomainProperty1
+  , testGroup "Primality"
+    [ testSmallAndQuick "Eisenstein primes found by the norm search used in\
+                        \ findPrime are really prime"
+                        findPrimesProperty1
+    , testSmallAndQuick "Eisenstein primes generated by `primes` are actually\
+                        \ primes"
+                        primesProperty1
+    , testSmallAndQuick "The infinite list of Eisenstein primes produced by\
+                        \ `primes` is ordered. "
+                        primesProperty2
+    , testSmallAndQuick "All generated primes are in the first sextant"
+                        primesProperty3
+    ]
+
+    , testGroup "Factorisation"
+      [ testSmallAndQuick "factorise produces correct results"
+                          factoriseProperty1
+      , testSmallAndQuick "factorise produces no factors with exponent 0"
+                          factoriseProperty2
+      , testSmallAndQuick "factorise produces no unit factors"
+                          factoriseProperty3
+      , testSmallAndQuick "factorise only produces primary primes"
+                          factoriseProperty4
+      , testCase          "factorise 15:+12" factoriseSpecialCase1
+      ]
+  ]