-- |
-- 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 PostfixOperators #-}
{-# 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)
import Data.List.Infinite (Infinite(..), (...))
import qualified Data.List.Infinite as Inf
import Data.List.NonEmpty (NonEmpty(..))
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
ω) *) 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 *) [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 :: Infinite (Prime EisensteinInteger)
primes :: Infinite (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 -> Infinite a -> Infinite a
:< forall a.
(a -> a -> Ordering) -> Infinite a -> Infinite a -> Infinite a
mergeBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing EisensteinInteger -> Integer
norm) Infinite EisensteinInteger
l Infinite EisensteinInteger
r
  where
    leftPrimes, rightPrimes :: Infinite (Prime Integer)
    (Infinite (Prime Integer)
leftPrimes, Infinite (Prime Integer)
rightPrimes) = forall a. (a -> Bool) -> Infinite a -> (Infinite a, Infinite a)
Inf.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' :: Infinite (Prime Integer)
    rightPrimes' :: Infinite (Prime Integer)
rightPrimes' = forall a. (a -> Bool) -> Infinite a -> Infinite a
Inf.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. Infinite a -> Infinite a
Inf.tail Infinite (Prime Integer)
rightPrimes

    l :: Infinite EisensteinInteger
    l :: Infinite EisensteinInteger
l = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Prime Integer
p -> forall a. Prime a -> a
unPrime Prime Integer
p Integer -> Integer -> EisensteinInteger
:+ Integer
0) Infinite (Prime Integer)
leftPrimes

    r :: Infinite EisensteinInteger
    r :: Infinite EisensteinInteger
r = forall a b. (a -> NonEmpty b) -> Infinite a -> Infinite b
Inf.concatMap
        (\Prime Integer
p -> let Integer
x :+ Integer
y = forall a. Prime a -> a
unPrime (Prime Integer -> Prime EisensteinInteger
findPrime Prime Integer
p) in (Integer
x Integer -> Integer -> EisensteinInteger
:+ Integer
y) forall a. a -> [a] -> NonEmpty a
:| [Integer
x Integer -> Integer -> EisensteinInteger
:+ (Integer
x forall a. Num a => a -> a -> a
- Integer
y)])
        Infinite (Prime Integer)
rightPrimes'

-- | [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