diff --git a/Math/NumberTheory/EisensteinIntegers.hs b/Math/NumberTheory/EisensteinIntegers.hs new file mode 100644 index 000000000..f84bc1d0c --- /dev/null +++ b/Math/NumberTheory/EisensteinIntegers.hs @@ -0,0 +1,364 @@ +-- | +-- Module: Math.NumberTheory.EisensteinIntegers +-- Copyright: (c) 2018 Alexandre Rodrigues Baldé +-- Licence: MIT +-- Maintainer: Alexandre Rodrigues Baldé +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- This module exports functions for manipulating Eisenstein integers, including +-- computing their prime factorisations. +-- + +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE RankNTypes #-} + +module Math.NumberTheory.EisensteinIntegers + ( EisensteinInteger(..) + , ω + , conjugate + , norm + , associates + , ids + + -- * Division and remainder functions + , divE + , divModE + , modE + , quotRemE + , quotE + , remE + + , divideByThree + , gcdE + + -- * Primality functions + , factorise + , findPrime + , isPrime + , primes + ) where + +import Data.List (mapAccumL, partition) +import Data.Maybe (fromMaybe) +import Data.Ord (comparing) +import GHC.Generics (Generic) + +import qualified Math.NumberTheory.Moduli as Moduli +import Math.NumberTheory.Moduli.Sqrt (FieldCharacteristic(..)) +import qualified Math.NumberTheory.Primes.Factorisation as Factorisation +import Math.NumberTheory.Primes.Types (PrimeNat(..)) +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 :+ + +-- | An Eisenstein integer is a + bω, where a and b are both integers. +data EisensteinInteger = (:+) { real :: !Integer, imag :: !Integer } + deriving (Eq, Ord, Generic) + +-- | The imaginary unit for Eisenstein integers, where +-- +-- > ω == (-1/2) + ((sqrt 3)/2)ι == exp(2*pi*ι/3) +-- and ι is the usual imaginary unit with ι² == -1. +ω :: EisensteinInteger +ω = 0 :+ 1 + +instance Show EisensteinInteger where + show (a :+ b) + | b == 0 = show a + | a == 0 = s ++ b' + | otherwise = show a ++ op ++ b' + where + b' = if abs b == 1 then "ω" else show (abs b) ++ "*ω" + op = if b > 0 then "+" else "-" + s = if b > 0 then "" else "-" + +instance Num EisensteinInteger where + (+) (a :+ b) (c :+ d) = (a + c) :+ (b + d) + (*) (a :+ b) (c :+ d) = (a * c - b * d) :+ (b * (c - d) + a * d) + abs = fst . absSignum + negate (a :+ b) = (-a) :+ (-b) + fromInteger n = n :+ 0 + signum = snd . absSignum + +-- | Returns an @EisensteinInteger@'s sign, and its associate in the first +-- sextant. +absSignum :: EisensteinInteger -> (EisensteinInteger, EisensteinInteger) +absSignum z@(a :+ b) + | a == 0 && b == 0 = (z, 0) -- origin + | a > b && b >= 0 = (z, 1) -- first sextant: 0 ≤ Arg(η) < π/3 + | b >= a && a > 0 = ((-ω) * z, 1 + ω) -- second sextant: π/3 ≤ Arg(η) < 2π/3 + | b > 0 && 0 >= a = ((-1 - ω) * z, ω) -- third sextant: 2π/3 ≤ Arg(η) < π + | a < b && b <= 0 = (- z, -1) -- fourth sextant: -π < Arg(η) < -2π/3 or Arg(η) = π + | b <= a && a < 0 = (ω * z, -1 - ω) -- fifth sextant: -2π/3 ≤ Arg(η) < -π/3 + | otherwise = ((1 + ω) * z, -ω) -- sixth sextant: -π/3 ≤ Arg(η) < 0 + +-- | List of all Eisenstein units, counterclockwise across all sextants, +-- starting with @1@. +ids :: [EisensteinInteger] +ids = take 6 (iterate ((1 + ω) *) 1) + +-- | Produce a list of an @EisensteinInteger@'s associates. +associates :: EisensteinInteger -> [EisensteinInteger] +associates e = map (e *) ids + +-- | Takes an Eisenstein prime whose norm is of the form @3k + 1@ with @k@ +-- a nonnegative integer, and return its primary associate. +-- * Does *not* check for this precondition. +-- * @head@ will fail when supplied a number unsatisfying it. +primary :: EisensteinInteger -> EisensteinInteger +primary = head . filter (\p -> p `modE` 3 == 2) . associates + +-- | Simultaneous 'quot' and 'rem'. +quotRemE + :: EisensteinInteger + -> EisensteinInteger + -> (EisensteinInteger, EisensteinInteger) +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. +divModE + :: EisensteinInteger + -> EisensteinInteger + -> (EisensteinInteger, EisensteinInteger) +divModE = divHelper div + +-- | Eisenstein integer division, truncating toward negative infinity. +divE :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger +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 + +-- | 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) + +-- | The square of the magnitude of a Eisenstein integer. +norm :: EisensteinInteger -> Integer +norm (a :+ b) = a*a - a * b + b*b + +-- | Checks if a given @EisensteinInteger@ is prime. @EisensteinInteger@s +-- whose norm is a prime congruent to @0@ or @1@ modulo 3 are prime. +-- See , +-- page 12. +isPrime :: EisensteinInteger -> Bool +isPrime e | e == 0 = False + -- Special case, @1 - ω@ is the only Eisenstein prime with norm @3@, + -- and @abs (1 - ω) = 2 + ω@. + | a' == 2 && b' == 1 = True + | b' == 0 && a' `mod` 3 == 2 = Testing.isPrime a' + | nE `mod` 3 == 1 = Testing.isPrime nE + | otherwise = False + where nE = norm e + a' :+ b' = abs e + +-- | Remove @1 - ω@ factors from an @EisensteinInteger@, and calculate that +-- prime's multiplicity in the number's factorisation. +divideByThree :: EisensteinInteger -> (Int, EisensteinInteger) +divideByThree = go 0 + where + go :: Int -> EisensteinInteger -> (Int, EisensteinInteger) + go !n z@(a :+ b) | r1 == 0 && r2 == 0 = go (n + 1) (q1 :+ q2) + | otherwise = (n, abs z) + where + -- @(a + a - b) :+ (a + b)@ is @z * (2 :+ 1)@, and @z * (2 :+ 1)/3@ + -- is the same as @z / (1 :+ (-1))@. + (q1, r1) = divMod (a + a - b) 3 + (q2, r2) = divMod (a + b) 3 + +-- | 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)) + +-- | Find an Eisenstein integer whose norm is the given prime number +-- in the form @3k + 1@ using a modification of the +-- . +-- +-- The maintainer +-- derived the following: +-- * Each prime of form @3n+1@ is actually of form @6k+1@. +-- * One has @(z+3k)^2 ≡ z^2 + 6kz + 9k^2 ≡ z^2 + (6k+1)z - z + 9k^2 ≡ z^2 - z + 9k^2 (mod 6k+1)@. +-- +-- * The goal is to solve @z^2 - z + 1 ≡ 0 (mod 6k+1)@. One has: +-- @z^2 - z + 9k^2 ≡ 9k^2 - 1 (mod 6k+1)@ +-- @(z+3k)^2 ≡ 9k^2-1 (mod 6k+1)@ +-- @z+3k = sqrtMod(9k^2-1)@ +-- @z = sqrtMod(9k^2-1) - 3k@ +-- +-- * For example, let @p = 7@, then @k = 1@. Square root of @9*1^2-1 modulo 7@ is @1@. +-- * And @z = 1 - 3*1 = -2 ≡ 5 (mod 7)@. +-- * Truly, @norm (5 :+ 1) = 25 - 5 + 1 = 21 ≡ 0 (mod 7)@. +findPrime :: Integer -> EisensteinInteger +findPrime p = case Moduli.sqrtModMaybe (9*k*k - 1) (FieldCharacteristic (PrimeNat . integerToNatural $ p) 1) of + Nothing -> error "findPrime: argument must be prime p = 6k + 1" + Just sqrtMod -> gcdE (p :+ 0) ((sqrtMod - 3 * k) :+ 1) + where + k :: Integer + k = p `div` 6 + +-- | An infinite list of Eisenstein primes. Uses primes in Z to exhaustively +-- generate all Eisenstein primes in order of ascending magnitude. +-- * Every prime is in the first sextant, so the list contains no associates. +-- * Eisenstein primes from the whole complex plane can be generated by +-- applying @associates@ to each prime in this list. +primes :: [EisensteinInteger] +primes = (2 :+ 1) : mergeBy (comparing norm) l r + where (leftPrimes, rightPrimes) = partition (\p -> p `mod` 3 == 2) Sieve.primes + rightPrimes' = filter (\prime -> prime `mod` 3 == 1) $ tail rightPrimes + l = [p :+ 0 | p <- leftPrimes] + r = [g | p <- rightPrimes', let x :+ y = findPrime p, g <- [x :+ y, x :+ (x - y)]] + +-- | Compute the prime factorisation of a Eisenstein integer. This is unique +-- up to units (+/- 1, +/- ω, +/- ω²). +-- * Unit factors are not included in the result. +-- * All prime factors are primary i.e. @e ≡ 2 (modE 3)@, for an Eisenstein +-- prime factor @e@. +-- +-- * This function works by factorising the norm of an Eisenstein integer +-- and then, for each prime factor, finding the Eisenstein prime whose norm +-- is said prime factor with @findPrime@. +-- +-- * This is only possible because the norm function of the Euclidean Domain of +-- Eisenstein integers is multiplicative: @norm (e1 * e2) == norm e1 * norm e2@ +-- for any two @EisensteinInteger@s @e1, e2@. +-- +-- * In the previously mentioned work , +-- in Theorem 8.4 in Chapter 8, a way is given to express any Eisenstein +-- integer @μ@ as @(-1)^a * ω^b * (1 - ω)^c * product [π_i^a_i | i <- [1..N]]@ +-- where @a, b, c, a_i@ are nonnegative integers, @N > 1@ is an integer and +-- @π_i@ are primary primes (for a primary Eisenstein prime @p@, +-- @p ≡ 2 (modE 3)@, see @primary@ above). +-- +-- * Aplying @norm@ to both sides of Theorem 8.4: +-- @norm μ = norm ((-1)^a * ω^b * (1 - ω)^c * product [ π_i^a_i | i <- [1..N]])@ +-- == @norm μ = norm ((-1)^a) * norm (ω^b) * norm ((1 - ω)^c) * norm (product [ π_i^a_i | i <- [1..N]])@ +-- == @norm μ = (norm (-1))^a * (norm ω)^b * (norm (1 - ω))^c * product [ norm (π_i^a_i) | i <- [1..N]]@ +-- == @norm μ = (norm (-1))^a * (norm ω)^b * (norm (1 - ω))^c * product [ (norm π_i)^a_i) | i <- [1..N]]@ +-- == @norm μ = 1^a * 1^b * 3^c * product [ (norm π_i)^a_i) | i <- [1..N]]@ +-- == @norm μ = 3^c * product [ (norm π_i)^a_i) | i <- [1..N]]@ +-- where @a, b, c, a_i@ are nonnegative integers, and @N > 1@ is an integer. +-- +-- * The remainder of the Eisenstein integer factorisation problem is about +-- finding appropriate @[e_i | i <- [1..M]@ such that +-- @(nub . map norm) [e_i | i <- [1..N]] == [π_i | i <- [1..N]]@ +-- where @ 1 < N <= M@ are integers, @nub@ removes duplicates and @==@ +-- is equality on sets. +-- +-- * The reason @M >= N@ is because the prime factors of an Eisenstein integer +-- may include a prime factor and its conjugate, meaning the number may have +-- more Eisenstein prime factors than its norm has integer prime factors. +factorise :: EisensteinInteger -> [(EisensteinInteger, Int)] +factorise g = concat $ + snd $ + mapAccumL go (abs g) (Factorisation.factorise $ norm g) + where + go :: EisensteinInteger -> (Integer, Int) -> (EisensteinInteger, [(EisensteinInteger, Int)]) + go z (3, e) | e == n = (q, [(2 :+ 1, e)]) + | otherwise = error $ "3 is a prime factor of the norm of z\ + \ == " ++ show z ++ " with multiplicity\ + \ " ++ show e ++ " but (1 - ω) only\ + \ divides z " ++ show n ++ "times." + where + -- Remove all @1 :+ (-1)@ (which is associated to @2 :+ 1@) factors + -- from the argument. + (n, q) = divideByThree z + go z (p, e) | p `mod` 3 == 2 = + let e' = e `quot` 2 in (z `quotI` (p ^ e'), [(p :+ 0, e')]) + + -- The @`mod` 3 == 0@ case need not be verified because the + -- only Eisenstein primes whose norm are a multiple of 3 + -- are @1 - ω@ and its associates, which have already been + -- removed by the above @go z (3, e)@ pattern match. + -- This @otherwise@ is mandatorily @`mod` 3 == 1@. + | otherwise = (z', filter ((> 0) . snd) [(gp, k), (gp', k')]) + where + gp@(x :+ y) = primary $ findPrime p + -- @gp'@ is @gp@'s conjugate. + gp' = primary $ abs $ x :+ (x - y) + (k, k', z') = divideByPrime gp gp' p e z + + quotI (a :+ b) n = (a `quot` n :+ b `quot` n) + +-- | Remove @p@ and @conjugate p@ factors from the argument, where +-- @p@ is an Eisenstein prime. +divideByPrime + :: EisensteinInteger -- ^ Eisenstein prime @p@ + -> EisensteinInteger -- ^ Conjugate of @p@ + -> Integer -- ^ Precomputed norm of @p@, of form @4k + 1@ + -> Int -- ^ Expected number of factors (either @p@ or @conjugate p@) + -- in Eisenstein integer @z@ + -> EisensteinInteger -- ^ Eisenstein integer @z@ + -> ( Int -- Multiplicity of factor @p@ in @z@ + , Int -- Multiplicity of factor @conjigate p@ in @z@ + , EisensteinInteger -- Remaining Eisenstein integer + ) +divideByPrime p p' np k = go k 0 + where + go :: Int -> Int -> EisensteinInteger -> (Int, Int, EisensteinInteger) + go 0 d z = (d, d, z) + go c d z | c >= 2, Just z' <- z `quotEvenI` np = go (c - 2) (d + 1) z' + go c d z = (d + d1, d + d2, z'') + where + (d1, z') = go1 c 0 z + d2 = c - d1 + z'' = head $ drop d2 + $ iterate (\g -> fromMaybe err $ (g * p) `quotEvenI` np) z' + + go1 :: Int -> Int -> EisensteinInteger -> (Int, EisensteinInteger) + go1 0 d z = (d, z) + go1 c d z + | Just z' <- (z * p') `quotEvenI` np + = go1 (c - 1) (d + 1) z' + | otherwise + = (d, z) + + err = error $ "divideByPrime: malformed arguments" ++ show (p, np, k) + +-- | Divide an Eisenstein integer by an even integer. +quotEvenI :: EisensteinInteger -> Integer -> Maybe EisensteinInteger +quotEvenI (x :+ y) n + | xr == 0 , yr == 0 = Just (xq :+ yq) + | otherwise = Nothing + where + (xq, xr) = x `quotRem` n + (yq, yr) = y `quotRem` n \ No newline at end of file diff --git a/Math/NumberTheory/GaussianIntegers.hs b/Math/NumberTheory/GaussianIntegers.hs index cb253c24d..295730330 100644 --- a/Math/NumberTheory/GaussianIntegers.hs +++ b/Math/NumberTheory/GaussianIntegers.hs @@ -46,6 +46,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 :+ @@ -73,17 +74,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) @@ -145,15 +147,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. diff --git a/Math/NumberTheory/UniqueFactorisation.hs b/Math/NumberTheory/UniqueFactorisation.hs index 8d3d1cd3f..8fb36ac99 100644 --- a/Math/NumberTheory/UniqueFactorisation.hs +++ b/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 \ No newline at end of file diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs index d4b831365..5f8abb962 100644 --- a/Math/NumberTheory/Utils.hs +++ b/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) \ No newline at end of file diff --git a/arithmoi.cabal b/arithmoi.cabal index 7d8993587..33cb01e44 100644 --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -57,6 +57,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 @@ -137,6 +138,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 @@ -189,6 +191,7 @@ benchmark criterion semigroups >=0.8 other-modules: Math.NumberTheory.ArithmeticFunctionsBench + Math.NumberTheory.EisensteinIntegersBench Math.NumberTheory.GaussianIntegersBench Math.NumberTheory.GCDBench Math.NumberTheory.JacobiBench diff --git a/benchmark/Bench.hs b/benchmark/Bench.hs index 592974145..5ee90857c 100644 --- a/benchmark/Bench.hs +++ b/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 @@ -15,6 +16,7 @@ import Math.NumberTheory.SmoothNumbersBench as SmoothNumbers main = defaultMain [ ArithmeticFunctions.benchSuite + , Eisenstein.benchSuite , Gaussian.benchSuite , GCD.benchSuite , Jacobi.benchSuite diff --git a/benchmark/Math/NumberTheory/EisensteinIntegersBench.hs b/benchmark/Math/NumberTheory/EisensteinIntegersBench.hs new file mode 100644 index 000000000..c942e115d --- /dev/null +++ b/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] + ] diff --git a/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs b/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs new file mode 100644 index 000000000..572efa651 --- /dev/null +++ b/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é 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 + ] + ] diff --git a/test-suite/Math/NumberTheory/TestUtils.hs b/test-suite/Math/NumberTheory/TestUtils.hs index 9a1dc85e2..f42b4d0c2 100644 --- a/test-suite/Math/NumberTheory/TestUtils.hs +++ b/test-suite/Math/NumberTheory/TestUtils.hs @@ -58,6 +58,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 @@ -70,6 +71,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) diff --git a/test-suite/Test.hs b/test-suite/Test.hs index 0cf8ff634..caf69377e 100644 --- a/test-suite/Test.hs +++ b/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 @@ -77,6 +79,7 @@ tests = testGroup "All" , Sieve.testSuite , Testing.testSuite ] + , Eisenstein.testSuite , Gaussian.testSuite , testGroup "ArithmeticFunctions" [ ArithmeticFunctions.testSuite