-- |
-- Module:      Math.NumberTheory.EisensteinIntegers
-- Copyright:   (c) 2018 Alexandre Rodrigues Baldé
-- Licence:     MIT
-- Maintainer:  Alexandre Rodrigues Baldé <alexandrer_b@outlook.com>
--
-- This module exports functions for manipulating Eisenstein integers, including
-- computing their prime factorisations.
--

{-# LANGUAGE BangPatterns   #-}
{-# LANGUAGE DeriveGeneric  #-}
{-# LANGUAGE RankNTypes     #-}
{-# LANGUAGE TypeFamilies   #-}

module Math.NumberTheory.Quadratic.EisensteinIntegers
  ( EisensteinInteger(..)
  , ω
  , conjugate
  , norm
  , associates
  , ids

  -- * Primality functions
  , findPrime
  , primes
  ) where

import Prelude hiding (quot, quotRem, gcd)
import Control.DeepSeq
import Data.Coerce
import Data.Euclidean
import Data.List                                       (mapAccumL, partition)
import Data.Maybe
import Data.Ord                                        (comparing)
import qualified Data.Semiring as S
import GHC.Generics                                    (Generic)

import Math.NumberTheory.Moduli.Sqrt
import Math.NumberTheory.Primes.Types
import qualified Math.NumberTheory.Primes as U
import Math.NumberTheory.Utils                          (mergeBy)
import Math.NumberTheory.Utils.FromIntegral

infix 6 :+

-- | An Eisenstein integer is @a + bω@, where @a@ and @b@ are both integers.
data EisensteinInteger = !Integer :+ !Integer
    deriving (EisensteinInteger -> EisensteinInteger -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: EisensteinInteger -> EisensteinInteger -> Bool
$c/= :: EisensteinInteger -> EisensteinInteger -> Bool
== :: EisensteinInteger -> EisensteinInteger -> Bool
$c== :: EisensteinInteger -> EisensteinInteger -> Bool
Eq, Eq EisensteinInteger
EisensteinInteger -> EisensteinInteger -> Bool
EisensteinInteger -> EisensteinInteger -> Ordering
EisensteinInteger -> EisensteinInteger -> EisensteinInteger
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
$cmin :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
max :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
$cmax :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
>= :: EisensteinInteger -> EisensteinInteger -> Bool
$c>= :: EisensteinInteger -> EisensteinInteger -> Bool
> :: EisensteinInteger -> EisensteinInteger -> Bool
$c> :: EisensteinInteger -> EisensteinInteger -> Bool
<= :: EisensteinInteger -> EisensteinInteger -> Bool
$c<= :: EisensteinInteger -> EisensteinInteger -> Bool
< :: EisensteinInteger -> EisensteinInteger -> Bool
$c< :: EisensteinInteger -> EisensteinInteger -> Bool
compare :: EisensteinInteger -> EisensteinInteger -> Ordering
$ccompare :: EisensteinInteger -> EisensteinInteger -> Ordering
Ord, forall x. Rep EisensteinInteger x -> EisensteinInteger
forall x. EisensteinInteger -> Rep EisensteinInteger x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cto :: forall x. Rep EisensteinInteger x -> EisensteinInteger
$cfrom :: forall x. EisensteinInteger -> Rep EisensteinInteger x
Generic)

instance NFData EisensteinInteger

-- | 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
ω :: EisensteinInteger
ω = Integer
0 Integer -> Integer -> EisensteinInteger
:+ Integer
1

instance Show EisensteinInteger where
    show :: EisensteinInteger -> String
show (Integer
a :+ Integer
b)
        | Integer
b forall a. Eq a => a -> a -> Bool
== Integer
0     = forall a. Show a => a -> String
show Integer
a
        | Integer
a forall a. Eq a => a -> a -> Bool
== Integer
0     = String
s forall a. [a] -> [a] -> [a]
++ String
b'
        | Bool
otherwise  = forall a. Show a => a -> String
show Integer
a forall a. [a] -> [a] -> [a]
++ String
op forall a. [a] -> [a] -> [a]
++ String
b'
        where
            b' :: String
b' = if forall a. Num a => a -> a
abs Integer
b forall a. Eq a => a -> a -> Bool
== Integer
1 then String
"ω" else forall a. Show a => a -> String
show (forall a. Num a => a -> a
abs Integer
b) forall a. [a] -> [a] -> [a]
++ String
"*ω"
            op :: String
op = if Integer
b forall a. Ord a => a -> a -> Bool
> Integer
0 then String
"+" else String
"-"
            s :: String
s  = if Integer
b forall a. Ord a => a -> a -> Bool
> Integer
0 then String
"" else String
"-"

instance Num EisensteinInteger where
    + :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
(+) (Integer
a :+ Integer
b) (Integer
c :+ Integer
d) = (Integer
a forall a. Num a => a -> a -> a
+ Integer
c) Integer -> Integer -> EisensteinInteger
:+ (Integer
b forall a. Num a => a -> a -> a
+ Integer
d)
    * :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
(*) (Integer
a :+ Integer
b) (Integer
c :+ Integer
d) = (Integer
a forall a. Num a => a -> a -> a
* Integer
c forall a. Num a => a -> a -> a
- Integer
b forall a. Num a => a -> a -> a
* Integer
d) Integer -> Integer -> EisensteinInteger
:+ (Integer
b forall a. Num a => a -> a -> a
* (Integer
c forall a. Num a => a -> a -> a
- Integer
d) forall a. Num a => a -> a -> a
+ Integer
a forall a. Num a => a -> a -> a
* Integer
d)
    abs :: EisensteinInteger -> EisensteinInteger
abs = forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. EisensteinInteger -> (EisensteinInteger, EisensteinInteger)
absSignum
    negate :: EisensteinInteger -> EisensteinInteger
negate (Integer
a :+ Integer
b) = (-Integer
a) Integer -> Integer -> EisensteinInteger
:+ (-Integer
b)
    fromInteger :: Integer -> EisensteinInteger
fromInteger Integer
n = Integer
n Integer -> Integer -> EisensteinInteger
:+ Integer
0
    signum :: EisensteinInteger -> EisensteinInteger
signum = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. EisensteinInteger -> (EisensteinInteger, EisensteinInteger)
absSignum

instance S.Semiring EisensteinInteger where
    plus :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
plus          = forall a. Num a => a -> a -> a
(+)
    times :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger
times         = forall a. Num a => a -> a -> a
(*)
    zero :: EisensteinInteger
zero          = Integer
0 Integer -> Integer -> EisensteinInteger
:+ Integer
0
    one :: EisensteinInteger
one           = Integer
1 Integer -> Integer -> EisensteinInteger
:+ Integer
0
    fromNatural :: Natural -> EisensteinInteger
fromNatural Natural
n = Natural -> Integer
naturalToInteger Natural
n Integer -> Integer -> EisensteinInteger
:+ Integer
0

instance S.Ring EisensteinInteger where
    negate :: EisensteinInteger -> EisensteinInteger
negate = forall a. Num a => a -> a
negate

-- | Returns an @EisensteinInteger@'s sign, and its associate in the first
-- sextant.
absSignum :: EisensteinInteger -> (EisensteinInteger, EisensteinInteger)
absSignum :: EisensteinInteger -> (EisensteinInteger, EisensteinInteger)
absSignum EisensteinInteger
0 = (EisensteinInteger
0, EisensteinInteger
0)
absSignum z :: EisensteinInteger
z@(Integer
a :+ Integer
b)
  -- first sextant: 0 ≤ Arg(z) < π/3
  | Integer
a forall a. Ord a => a -> a -> Bool
> Integer
b Bool -> Bool -> Bool
&& Integer
b forall a. Ord a => a -> a -> Bool
>= Integer
0 = (EisensteinInteger
z, EisensteinInteger
1)
  -- second sextant: π/3 ≤ Arg(z) < 2π/3
  | Integer
b forall a. Ord a => a -> a -> Bool
>= Integer
a Bool -> Bool -> Bool
&& Integer
a forall a. Ord a => a -> a -> Bool
> Integer
0 = (Integer
b Integer -> Integer -> EisensteinInteger
:+ (Integer
b forall a. Num a => a -> a -> a
- Integer
a), Integer
1 Integer -> Integer -> EisensteinInteger
:+ Integer
1)
  -- third sextant: 2π/3 ≤ Arg(z) < π
  | Integer
b forall a. Ord a => a -> a -> Bool
> Integer
0 Bool -> Bool -> Bool
&& Integer
0 forall a. Ord a => a -> a -> Bool
>= Integer
a = ((Integer
b forall a. Num a => a -> a -> a
- Integer
a) Integer -> Integer -> EisensteinInteger
:+ (-Integer
a), Integer
0 Integer -> Integer -> EisensteinInteger
:+ Integer
1)
  -- fourth sextant: -π ≤ Arg(z) < -2π/3
  | Integer
a forall a. Ord a => a -> a -> Bool
< Integer
b Bool -> Bool -> Bool
&& Integer
b forall a. Ord a => a -> a -> Bool
<= Integer
0 = (-EisensteinInteger
z, -EisensteinInteger
1)
  -- fifth sextant: -2π/3 ≤ Arg(η) < -π/3
  | Integer
b forall a. Ord a => a -> a -> Bool
<= Integer
a Bool -> Bool -> Bool
&& Integer
a forall a. Ord a => a -> a -> Bool
< Integer
0 = ((-Integer
b) Integer -> Integer -> EisensteinInteger
:+ (Integer
a forall a. Num a => a -> a -> a
- Integer
b), (-Integer
1) Integer -> Integer -> EisensteinInteger
:+ (-Integer
1))
  -- sixth sextant: -π/3 ≤ Arg(η) < 0
  | Bool
otherwise       = ((Integer
a forall a. Num a => a -> a -> a
- Integer
b) Integer -> Integer -> EisensteinInteger
:+ Integer
a, Integer
0 Integer -> Integer -> EisensteinInteger
:+ (-Integer
1))

-- | List of all Eisenstein units, counterclockwise across all sextants,
-- starting with @1@.
ids :: [EisensteinInteger]
ids :: [EisensteinInteger]
ids = forall a. Int -> [a] -> [a]
take Int
6 (forall a. (a -> a) -> a -> [a]
iterate ((EisensteinInteger
1 forall a. Num a => a -> a -> a
+ EisensteinInteger
ω) forall a. Num a => a -> a -> a
*) EisensteinInteger
1)

-- | Produce a list of an @EisensteinInteger@'s associates.
associates :: EisensteinInteger -> [EisensteinInteger]
associates :: EisensteinInteger -> [EisensteinInteger]
associates EisensteinInteger
e = forall a b. (a -> b) -> [a] -> [b]
map (EisensteinInteger
e forall a. Num a => a -> a -> a
*) [EisensteinInteger]
ids

instance GcdDomain EisensteinInteger

instance Euclidean EisensteinInteger where
  degree :: EisensteinInteger -> Natural
degree = forall a. Num a => Integer -> a
fromInteger forall b c a. (b -> c) -> (a -> b) -> a -> c
. EisensteinInteger -> Integer
norm
  quotRem :: EisensteinInteger
-> EisensteinInteger -> (EisensteinInteger, EisensteinInteger)
quotRem EisensteinInteger
x (Integer
d :+ Integer
0) = EisensteinInteger
-> Integer -> (EisensteinInteger, EisensteinInteger)
quotRemInt EisensteinInteger
x Integer
d
  quotRem EisensteinInteger
x EisensteinInteger
y = (EisensteinInteger
q, EisensteinInteger
x forall a. Num a => a -> a -> a
- EisensteinInteger
q forall a. Num a => a -> a -> a
* EisensteinInteger
y)
    where
      (EisensteinInteger
q, EisensteinInteger
_) = EisensteinInteger
-> Integer -> (EisensteinInteger, EisensteinInteger)
quotRemInt (EisensteinInteger
x forall a. Num a => a -> a -> a
* EisensteinInteger -> EisensteinInteger
conjugate EisensteinInteger
y) (EisensteinInteger -> Integer
norm EisensteinInteger
y)

quotRemInt :: EisensteinInteger -> Integer -> (EisensteinInteger, EisensteinInteger)
quotRemInt :: EisensteinInteger
-> Integer -> (EisensteinInteger, EisensteinInteger)
quotRemInt EisensteinInteger
z   Integer
1  = ( EisensteinInteger
z, EisensteinInteger
0)
quotRemInt EisensteinInteger
z (-1) = (-EisensteinInteger
z, EisensteinInteger
0)
quotRemInt (Integer
a :+ Integer
b) Integer
c = (Integer
qa Integer -> Integer -> EisensteinInteger
:+ Integer
qb, (Integer
ra forall a. Num a => a -> a -> a
- Integer
bumpA) Integer -> Integer -> EisensteinInteger
:+ (Integer
rb forall a. Num a => a -> a -> a
- Integer
bumpB))
  where
    halfC :: Integer
halfC    = forall a. Num a => a -> a
abs Integer
c forall a. Euclidean a => a -> a -> a
`quot` Integer
2
    bumpA :: Integer
bumpA    = forall a. Num a => a -> a
signum Integer
a forall a. Num a => a -> a -> a
* Integer
halfC
    bumpB :: Integer
bumpB    = forall a. Num a => a -> a
signum Integer
b forall a. Num a => a -> a -> a
* Integer
halfC
    (Integer
qa, Integer
ra) = (Integer
a forall a. Num a => a -> a -> a
+ Integer
bumpA) forall a. Euclidean a => a -> a -> (a, a)
`quotRem` Integer
c
    (Integer
qb, Integer
rb) = (Integer
b forall a. Num a => a -> a -> a
+ Integer
bumpB) forall a. Euclidean a => a -> a -> (a, a)
`quotRem` Integer
c

-- | Conjugate a Eisenstein integer.
conjugate :: EisensteinInteger -> EisensteinInteger
conjugate :: EisensteinInteger -> EisensteinInteger
conjugate (Integer
a :+ Integer
b) = (Integer
a forall a. Num a => a -> a -> a
- Integer
b) Integer -> Integer -> EisensteinInteger
:+ (-Integer
b)

-- | The square of the magnitude of a Eisenstein integer.
norm :: EisensteinInteger -> Integer
norm :: EisensteinInteger -> Integer
norm (Integer
a :+ Integer
b) = Integer
aforall a. Num a => a -> a -> a
*Integer
a forall a. Num a => a -> a -> a
- Integer
a forall a. Num a => a -> a -> a
* Integer
b forall a. Num a => a -> a -> a
+ Integer
bforall a. Num a => a -> a -> a
*Integer
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 <http://thekeep.eiu.edu/theses/2467 Bandara, Sarada, "An Exposition of the Eisenstein Integers" (2016)>,
-- page 12.
isPrime :: EisensteinInteger -> Bool
isPrime :: EisensteinInteger -> Bool
isPrime EisensteinInteger
e | EisensteinInteger
e forall a. Eq a => a -> a -> Bool
== EisensteinInteger
0                     = Bool
False
          -- Special case, @1 - ω@ is the only Eisenstein prime with norm @3@,
          --  and @abs (1 - ω) = 2 + ω@.
          | Integer
a' forall a. Eq a => a -> a -> Bool
== Integer
2 Bool -> Bool -> Bool
&& Integer
b' forall a. Eq a => a -> a -> Bool
== Integer
1         = Bool
True
          | Integer
b' forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
a' forall a. Integral a => a -> a -> a
`mod` Integer
3 forall a. Eq a => a -> a -> Bool
== Integer
2 = forall a. Maybe a -> Bool
isJust forall a b. (a -> b) -> a -> b
$ forall a. UniqueFactorisation a => a -> Maybe (Prime a)
U.isPrime Integer
a'
          | Integer
nE forall a. Integral a => a -> a -> a
`mod` Integer
3 forall a. Eq a => a -> a -> Bool
== Integer
1            = forall a. Maybe a -> Bool
isJust forall a b. (a -> b) -> a -> b
$ forall a. UniqueFactorisation a => a -> Maybe (Prime a)
U.isPrime Integer
nE
          | Bool
otherwise = Bool
False
  where nE :: Integer
nE       = EisensteinInteger -> Integer
norm EisensteinInteger
e
        Integer
a' :+ Integer
b' = forall a. Num a => a -> a
abs EisensteinInteger
e

-- | Remove @1 - ω@ factors from an @EisensteinInteger@, and calculate that
-- prime's multiplicity in the number's factorisation.
divideByThree :: EisensteinInteger -> (Word, EisensteinInteger)
divideByThree :: EisensteinInteger -> (Word, EisensteinInteger)
divideByThree = Word -> EisensteinInteger -> (Word, EisensteinInteger)
go Word
0
  where
    go :: Word -> EisensteinInteger -> (Word, EisensteinInteger)
    go :: Word -> EisensteinInteger -> (Word, EisensteinInteger)
go !Word
n z :: EisensteinInteger
z@(Integer
a :+ Integer
b) | Integer
r1 forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
r2 forall a. Eq a => a -> a -> Bool
== Integer
0 = Word -> EisensteinInteger -> (Word, EisensteinInteger)
go (Word
n forall a. Num a => a -> a -> a
+ Word
1) (Integer
q1 Integer -> Integer -> EisensteinInteger
:+ Integer
q2)
                     | Bool
otherwise          = (Word
n, forall a. Num a => a -> a
abs EisensteinInteger
z)
      where
        -- @(a + a - b) :+ (a + b)@ is @z * (2 :+ 1)@, and @z * (2 :+ 1)/3@
        -- is the same as @z / (1 :+ (-1))@.
        (Integer
q1, Integer
r1) = forall a. Integral a => a -> a -> (a, a)
divMod (Integer
a forall a. Num a => a -> a -> a
+ Integer
a forall a. Num a => a -> a -> a
- Integer
b) Integer
3
        (Integer
q2, Integer
r2) = forall a. Integral a => a -> a -> (a, a)
divMod (Integer
a forall a. Num a => a -> a -> a
+ Integer
b) Integer
3

-- | Find an Eisenstein integer whose norm is the given prime number
-- in the form @3k + 1@.
--
-- >>> import Math.NumberTheory.Primes (nextPrime)
-- >>> findPrime (nextPrime 7)
-- Prime 3+2*ω
findPrime :: Prime Integer -> U.Prime EisensteinInteger
findPrime :: Prime Integer -> Prime EisensteinInteger
findPrime Prime Integer
p = case (Integer
r, Integer -> Prime Integer -> [Integer]
sqrtsModPrime (Integer
9 forall a. Num a => a -> a -> a
* Integer
q forall a. Num a => a -> a -> a
* Integer
q forall a. Num a => a -> a -> a
- Integer
1) Prime Integer
p) of
  (Integer
1, Integer
z : [Integer]
_) -> forall a. a -> Prime a
Prime forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ forall a. GcdDomain a => a -> a -> a
gcd (forall a. Prime a -> a
unPrime Prime Integer
p Integer -> Integer -> EisensteinInteger
:+ Integer
0) ((Integer
z forall a. Num a => a -> a -> a
- Integer
3 forall a. Num a => a -> a -> a
* Integer
q) Integer -> Integer -> EisensteinInteger
:+ Integer
1)
  (Integer, [Integer])
_ -> forall a. HasCallStack => String -> a
error String
"findPrime: argument must be prime p = 6k + 1"
  where
    (Integer
q, Integer
r) = forall a. Prime a -> a
unPrime Prime Integer
p forall a. Euclidean a => a -> a -> (a, a)
`quotRem` Integer
6

-- | An infinite list of Eisenstein primes. Uses primes in @Z@ to exhaustively
-- generate all Eisenstein primes in order of ascending norm.
--
-- * 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.
--
-- >>> take 10 primes
-- [Prime 2+ω,Prime 2,Prime 3+2*ω,Prime 3+ω,Prime 4+3*ω,Prime 4+ω,Prime 5+3*ω,Prime 5+2*ω,Prime 5,Prime 6+5*ω]
primes :: [Prime EisensteinInteger]
primes :: [Prime EisensteinInteger]
primes = coerce :: forall a b. Coercible a b => a -> b
coerce forall a b. (a -> b) -> a -> b
$ (Integer
2 Integer -> Integer -> EisensteinInteger
:+ Integer
1) forall a. a -> [a] -> [a]
: forall a. (a -> a -> Ordering) -> [a] -> [a] -> [a]
mergeBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing EisensteinInteger -> Integer
norm) [EisensteinInteger]
l [EisensteinInteger]
r
  where
    leftPrimes, rightPrimes :: [Prime Integer]
    ([Prime Integer]
leftPrimes, [Prime Integer]
rightPrimes) = forall a. (a -> Bool) -> [a] -> ([a], [a])
partition (\Prime Integer
p -> forall a. Prime a -> a
unPrime Prime Integer
p forall a. Integral a => a -> a -> a
`mod` Integer
3 forall a. Eq a => a -> a -> Bool
== Integer
2) [forall a.
(Bits a, Integral a, UniqueFactorisation a) =>
a -> Prime a
U.nextPrime Integer
2 ..]
    rightPrimes' :: [Prime Integer]
rightPrimes' = forall a. (a -> Bool) -> [a] -> [a]
filter (\Prime Integer
prime -> forall a. Prime a -> a
unPrime Prime Integer
prime forall a. Integral a => a -> a -> a
`mod` Integer
3 forall a. Eq a => a -> a -> Bool
== Integer
1) forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
tail [Prime Integer]
rightPrimes
    l :: [EisensteinInteger]
l = [forall a. Prime a -> a
unPrime Prime Integer
p Integer -> Integer -> EisensteinInteger
:+ Integer
0 | Prime Integer
p <- [Prime Integer]
leftPrimes]
    r :: [EisensteinInteger]
r = [EisensteinInteger
g | Prime Integer
p <- [Prime Integer]
rightPrimes', let Integer
x :+ Integer
y = forall a. Prime a -> a
unPrime (Prime Integer -> Prime EisensteinInteger
findPrime Prime Integer
p), EisensteinInteger
g <- [Integer
x Integer -> Integer -> EisensteinInteger
:+ Integer
y, Integer
x Integer -> Integer -> EisensteinInteger
:+ (Integer
x forall a. Num a => a -> a -> a
- Integer
y)]]

-- | [Implementation notes for factorise function]
--
-- Compute the prime factorisation of a Eisenstein integer.
--
--     1. 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@.
--     2. 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@.
--     3. In the previously mentioned work <http://thekeep.eiu.edu/theses/2467 Bandara, Sarada, "An Exposition of the Eisenstein Integers" (2016)>,
--        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 Eisenstein primes.
--
-- Aplying @norm@ to both sides of the equation from Theorem 8.4:
--
-- 1. @norm μ = norm ( (-1)^a * ω^b * (1 - ω)^c * product [ π_i^a_i | i <- [1..N]] ) ==@
-- 2. @norm μ = norm ((-1)^a) * norm (ω^b) * norm ((1 - ω)^c) * norm (product [ π_i^a_i | i <- [1..N]]) ==@
-- 3. @norm μ = (norm (-1))^a * (norm ω)^b * (norm (1 - ω))^c * product [ norm (π_i^a_i) | i <- [1..N]] ==@
-- 4. @norm μ = (norm (-1))^a * (norm ω)^b * (norm (1 - ω))^c * product [ (norm π_i)^a_i) | i <- [1..N]] ==@
-- 5. @norm μ = 1^a * 1^b * 3^c * product [ (norm π_i)^a_i) | i <- [1..N]] ==@
-- 6. @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 Eisenstein primes @[e_i | i <- [1..M]]@ such that
-- @map norm [e_i | i <- [1..M]] == map norm [π_i | i <- [1..N]]@
-- where @ 1 < N <= M@ are integers and @==@ is equality on sets
-- (i.e.duplicates do not matter).
--
-- NB: The reason @M >= N@ is because the prime factors of an Eisenstein integer
-- may include a prime factor and its conjugate (both have the same norm),
-- meaning the number may have more Eisenstein prime factors than its norm has
-- integer prime factors.
factorise :: EisensteinInteger -> [(Prime EisensteinInteger, Word)]
factorise :: EisensteinInteger -> [(Prime EisensteinInteger, Word)]
factorise EisensteinInteger
g = forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat forall a b. (a -> b) -> a -> b
$
              forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$
              forall (t :: * -> *) s a b.
Traversable t =>
(s -> a -> (s, b)) -> s -> t a -> (s, t b)
mapAccumL EisensteinInteger
-> (Prime Integer, Word)
-> (EisensteinInteger, [(Prime EisensteinInteger, Word)])
go (forall a. Num a => a -> a
abs EisensteinInteger
g) (forall a. UniqueFactorisation a => a -> [(Prime a, Word)]
U.factorise forall a b. (a -> b) -> a -> b
$ EisensteinInteger -> Integer
norm EisensteinInteger
g)
  where
    go :: EisensteinInteger -> (Prime Integer, Word) -> (EisensteinInteger, [(Prime EisensteinInteger, Word)])
    go :: EisensteinInteger
-> (Prime Integer, Word)
-> (EisensteinInteger, [(Prime EisensteinInteger, Word)])
go EisensteinInteger
z (Prime Integer
3, Word
e)
      | Word
e forall a. Eq a => a -> a -> Bool
== Word
n    = (EisensteinInteger
q, [(forall a. a -> Prime a
Prime (Integer
2 Integer -> Integer -> EisensteinInteger
:+ Integer
1), Word
e)])
      | Bool
otherwise = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"3 is a prime factor of the norm of z = " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show EisensteinInteger
z
                          forall a. [a] -> [a] -> [a]
++ String
" with multiplicity " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Word
e
                          forall a. [a] -> [a] -> [a]
++ String
" but (1 - ω) only divides z " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Word
n forall a. [a] -> [a] -> [a]
++ String
"times."
      where
        -- Remove all @1 :+ (-1)@ (which is associated to @2 :+ 1@) factors
        -- from the argument.
        (Word
n, EisensteinInteger
q) = EisensteinInteger -> (Word, EisensteinInteger)
divideByThree EisensteinInteger
z
    go EisensteinInteger
z (Prime Integer
p, Word
e)
      | forall a. Prime a -> a
unPrime Prime Integer
p forall a. Integral a => a -> a -> a
`mod` Integer
3 forall a. Eq a => a -> a -> Bool
== Integer
2
      = let e' :: Word
e' = Word
e forall a. Euclidean a => a -> a -> a
`quot` Word
2 in (EisensteinInteger
z EisensteinInteger -> Integer -> EisensteinInteger
`quotI` (forall a. Prime a -> a
unPrime Prime Integer
p forall a b. (Num a, Integral b) => a -> b -> a
^ Word
e'), [(forall a. a -> Prime a
Prime (forall a. Prime a -> a
unPrime Prime Integer
p Integer -> Integer -> EisensteinInteger
:+ Integer
0), Word
e')])

      -- The @`rem` 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@.
      | Bool
otherwise   = (EisensteinInteger
z', forall a. (a -> Bool) -> [a] -> [a]
filter ((forall a. Ord a => a -> a -> Bool
> Word
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd) [(Prime EisensteinInteger
gp, Word
k), (Prime EisensteinInteger
gp', Word
k')])
      where
        gp :: Prime EisensteinInteger
gp = Prime Integer -> Prime EisensteinInteger
findPrime Prime Integer
p
        Integer
x :+ Integer
y = forall a. Prime a -> a
unPrime Prime EisensteinInteger
gp
        -- @gp'@ is @gp@'s conjugate.
        gp' :: Prime EisensteinInteger
gp' = forall a. a -> Prime a
Prime (Integer
x Integer -> Integer -> EisensteinInteger
:+ (Integer
x forall a. Num a => a -> a -> a
- Integer
y))
        (Word
k, Word
k', EisensteinInteger
z') = Prime EisensteinInteger
-> Prime EisensteinInteger
-> Integer
-> Word
-> EisensteinInteger
-> (Word, Word, EisensteinInteger)
divideByPrime Prime EisensteinInteger
gp Prime EisensteinInteger
gp' (forall a. Prime a -> a
unPrime Prime Integer
p) Word
e EisensteinInteger
z

        quotI :: EisensteinInteger -> Integer -> EisensteinInteger
quotI (Integer
a :+ Integer
b) Integer
n = Integer
a forall a. Euclidean a => a -> a -> a
`quot` Integer
n Integer -> Integer -> EisensteinInteger
:+ Integer
b forall a. Euclidean a => a -> a -> a
`quot` Integer
n

-- | Remove @p@ and @conjugate p@ factors from the argument, where
-- @p@ is an Eisenstein prime.
divideByPrime
    :: Prime EisensteinInteger -- ^ Eisenstein prime @p@
    -> Prime EisensteinInteger -- ^ Conjugate of @p@
    -> Integer                 -- ^ Precomputed norm of @p@, of form @4k + 1@
    -> Word                    -- ^ Expected number of factors (either @p@ or @conjugate p@)
                               --   in Eisenstein integer @z@
    -> EisensteinInteger       -- ^ Eisenstein integer @z@
    -> ( Word                  -- Multiplicity of factor @p@ in @z@
       , Word                  -- Multiplicity of factor @conjigate p@ in @z@
       , EisensteinInteger     -- Remaining Eisenstein integer
       )
divideByPrime :: Prime EisensteinInteger
-> Prime EisensteinInteger
-> Integer
-> Word
-> EisensteinInteger
-> (Word, Word, EisensteinInteger)
divideByPrime Prime EisensteinInteger
p Prime EisensteinInteger
p' Integer
np Word
k = Word
-> Word -> EisensteinInteger -> (Word, Word, EisensteinInteger)
go Word
k Word
0
    where
        go :: Word -> Word -> EisensteinInteger -> (Word, Word, EisensteinInteger)
        go :: Word
-> Word -> EisensteinInteger -> (Word, Word, EisensteinInteger)
go Word
0 Word
d EisensteinInteger
z = (Word
d, Word
d, EisensteinInteger
z)
        go Word
c Word
d EisensteinInteger
z | Word
c forall a. Ord a => a -> a -> Bool
>= Word
2, Just EisensteinInteger
z' <- EisensteinInteger
z EisensteinInteger -> Integer -> Maybe EisensteinInteger
`quotEvenI` Integer
np = Word
-> Word -> EisensteinInteger -> (Word, Word, EisensteinInteger)
go (Word
c forall a. Num a => a -> a -> a
- Word
2) (Word
d forall a. Num a => a -> a -> a
+ Word
1) EisensteinInteger
z'
        go Word
c Word
d EisensteinInteger
z = (Word
d forall a. Num a => a -> a -> a
+ Word
d1, Word
d forall a. Num a => a -> a -> a
+ Word
d2, EisensteinInteger
z'')
            where
                (Word
d1, EisensteinInteger
z') = Word -> Word -> EisensteinInteger -> (Word, EisensteinInteger)
go1 Word
c Word
0 EisensteinInteger
z
                d2 :: Word
d2 = Word
c forall a. Num a => a -> a -> a
- Word
d1
                z'' :: EisensteinInteger
z'' = forall a. (a -> a) -> a -> [a]
iterate (\EisensteinInteger
g -> forall a. a -> Maybe a -> a
fromMaybe EisensteinInteger
err forall a b. (a -> b) -> a -> b
$ (EisensteinInteger
g forall a. Num a => a -> a -> a
* forall a. Prime a -> a
unPrime Prime EisensteinInteger
p) EisensteinInteger -> Integer -> Maybe EisensteinInteger
`quotEvenI` Integer
np) EisensteinInteger
z' forall a. [a] -> Int -> a
!! forall a. Ord a => a -> a -> a
max Int
0 (Word -> Int
wordToInt Word
d2)

        go1 :: Word -> Word -> EisensteinInteger -> (Word, EisensteinInteger)
        go1 :: Word -> Word -> EisensteinInteger -> (Word, EisensteinInteger)
go1 Word
0 Word
d EisensteinInteger
z = (Word
d, EisensteinInteger
z)
        go1 Word
c Word
d EisensteinInteger
z
            | Just EisensteinInteger
z' <- (EisensteinInteger
z forall a. Num a => a -> a -> a
* forall a. Prime a -> a
unPrime Prime EisensteinInteger
p') EisensteinInteger -> Integer -> Maybe EisensteinInteger
`quotEvenI` Integer
np
            = Word -> Word -> EisensteinInteger -> (Word, EisensteinInteger)
go1 (Word
c forall a. Num a => a -> a -> a
- Word
1) (Word
d forall a. Num a => a -> a -> a
+ Word
1) EisensteinInteger
z'
            | Bool
otherwise
            = (Word
d, EisensteinInteger
z)

        err :: EisensteinInteger
err = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"divideByPrime: malformed arguments" forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (Prime EisensteinInteger
p, Integer
np, Word
k)

-- | Divide an Eisenstein integer by an even integer.
quotEvenI :: EisensteinInteger -> Integer -> Maybe EisensteinInteger
quotEvenI :: EisensteinInteger -> Integer -> Maybe EisensteinInteger
quotEvenI (Integer
x :+ Integer
y) Integer
n
    | Integer
xr forall a. Eq a => a -> a -> Bool
== Integer
0 , Integer
yr forall a. Eq a => a -> a -> Bool
== Integer
0 = forall a. a -> Maybe a
Just (Integer
xq Integer -> Integer -> EisensteinInteger
:+ Integer
yq)
    | Bool
otherwise         = forall a. Maybe a
Nothing
  where
    (Integer
xq, Integer
xr) = Integer
x forall a. Euclidean a => a -> a -> (a, a)
`quotRem` Integer
n
    (Integer
yq, Integer
yr) = Integer
y forall a. Euclidean a => a -> a -> (a, a)
`quotRem` Integer
n

-------------------------------------------------------------------------------

-- | See the source code and Haddock comments for the @factorise@ and @isPrime@
-- functions in this module (they are not exported) for implementation
-- details.
instance U.UniqueFactorisation EisensteinInteger where
  factorise :: EisensteinInteger -> [(Prime EisensteinInteger, Word)]
factorise EisensteinInteger
0 = []
  factorise EisensteinInteger
e = coerce :: forall a b. Coercible a b => a -> b
coerce forall a b. (a -> b) -> a -> b
$ EisensteinInteger -> [(Prime EisensteinInteger, Word)]
factorise EisensteinInteger
e

  isPrime :: EisensteinInteger -> Maybe (Prime EisensteinInteger)
isPrime EisensteinInteger
e = if EisensteinInteger -> Bool
isPrime EisensteinInteger
e then forall a. a -> Maybe a
Just (forall a. a -> Prime a
Prime EisensteinInteger
e) else forall a. Maybe a
Nothing