{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE TypeFamilies #-}
module Math.NumberTheory.Quadratic.GaussianIntegers (
GaussianInteger(..),
ι,
conjugate,
norm,
primes,
findPrime,
) where
import Prelude hiding (quot, quotRem)
import Control.DeepSeq (NFData)
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
import Math.NumberTheory.Moduli.Sqrt
import Math.NumberTheory.Roots (integerSquareRoot)
import Math.NumberTheory.Primes.Types
import qualified Math.NumberTheory.Primes as U
import Math.NumberTheory.Utils (mergeBy)
import Math.NumberTheory.Utils.FromIntegral
infix 6 :+
data GaussianInteger = (:+) { real :: !Integer, imag :: !Integer }
deriving (Eq, Ord, Generic)
instance NFData GaussianInteger
ι :: GaussianInteger
ι = 0 :+ 1
instance Show GaussianInteger 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 GaussianInteger where
(+) (a :+ b) (c :+ d) = (a + c) :+ (b + d)
(*) (a :+ b) (c :+ d) = (a * c - b * d) :+ (a * d + b * c)
abs = fst . absSignum
negate (a :+ b) = (-a) :+ (-b)
fromInteger n = n :+ 0
signum = snd . absSignum
instance S.Semiring GaussianInteger where
plus = (+)
times = (*)
zero = 0 :+ 0
one = 1 :+ 0
fromNatural n = fromIntegral n :+ 0
instance S.Ring GaussianInteger where
negate = negate
absSignum :: GaussianInteger -> (GaussianInteger, GaussianInteger)
absSignum z@(a :+ b)
| a == 0 && b == 0 = (z, 0)
| a > 0 && b >= 0 = (z, 1)
| a <= 0 && b > 0 = (b :+ (-a), ι)
| a < 0 && b <= 0 = ((-a) :+ (-b), -1)
| otherwise = ((-b) :+ a, -ι)
instance GcdDomain GaussianInteger
instance Euclidean GaussianInteger where
degree = fromInteger . norm
quotRem = divHelper
divHelper
:: GaussianInteger
-> GaussianInteger
-> (GaussianInteger, GaussianInteger)
divHelper g h = (q, r)
where
nr :+ ni = g * conjugate h
denom = norm h
q = ((nr + signum nr * denom `quot` 2) `quot` denom) :+ ((ni + signum ni * denom `quot` 2) `quot` denom)
r = g - h * q
conjugate :: GaussianInteger -> GaussianInteger
conjugate (r :+ i) = r :+ (-i)
norm :: GaussianInteger -> Integer
norm (x :+ y) = x * x + y * y
isPrime :: GaussianInteger -> Bool
isPrime g@(x :+ y)
| x == 0 && y /= 0 = abs y `mod` 4 == 3 && isJust (U.isPrime y)
| y == 0 && x /= 0 = abs x `mod` 4 == 3 && isJust (U.isPrime x)
| otherwise = isJust $ U.isPrime $ norm g
primes :: [U.Prime GaussianInteger]
primes = coerce $ (1 :+ 1) : mergeBy (comparing norm) l r
where
leftPrimes, rightPrimes :: [Prime Integer]
(leftPrimes, rightPrimes) = partition (\p -> unPrime p `mod` 4 == 3) [U.nextPrime 3 ..]
l = [unPrime p :+ 0 | p <- leftPrimes]
r = [g | p <- rightPrimes, let Prime (x :+ y) = findPrime p, g <- [x :+ y, y :+ x]]
findPrime :: Prime Integer -> U.Prime GaussianInteger
findPrime p = case sqrtsModPrime (-1) p of
[] -> error "findPrime: an argument must be prime p = 4k + 1"
z : _ -> Prime $ go (unPrime p) z
where
sqrtp :: Integer
sqrtp = integerSquareRoot (unPrime p)
go :: Integer -> Integer -> GaussianInteger
go g h
| g <= sqrtp = g :+ h
| otherwise = go h (g `mod` h)
factorise :: GaussianInteger -> [(Prime GaussianInteger, Word)]
factorise g = concat $ snd $ mapAccumL go g (U.factorise $ norm g)
where
go :: GaussianInteger -> (Prime Integer, Word) -> (GaussianInteger, [(Prime GaussianInteger, Word)])
go z (Prime 2, e) = (divideByTwo z, [(Prime (1 :+ 1), e)])
go z (p, e)
| unPrime p `mod` 4 == 3
= let e' = e `quot` 2 in (z `quotI` (unPrime p ^ e'), [(Prime (unPrime p :+ 0), e')])
| otherwise
= (z', filter ((> 0) . snd) [(gp, k), (gp', k')])
where
gp = findPrime p
(k, k', z') = divideByPrime gp (unPrime p) e z
gp' = Prime (abs (conjugate (unPrime gp)))
divideByTwo :: GaussianInteger -> GaussianInteger
divideByTwo z@(x :+ y)
| even x, even y
= divideByTwo $ z `quotI` 2
| odd x, odd y
= (x - y) `quot` 2 :+ (x + y) `quot` 2
| otherwise
= z
divideByPrime
:: Prime GaussianInteger
-> Integer
-> Word
-> GaussianInteger
-> ( Word
, Word
, GaussianInteger
)
divideByPrime p np k = go k 0
where
go :: Word -> Word -> GaussianInteger -> (Word, Word, GaussianInteger)
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 (wordToInt d2)
$ iterate (\g -> fromMaybe err $ (g * unPrime p) `quotEvenI` np) z'
go1 :: Word -> Word -> GaussianInteger -> (Word, GaussianInteger)
go1 0 d z = (d, z)
go1 c d z
| Just z' <- (z * conjugate (unPrime p)) `quotEvenI` np
= go1 (c - 1) (d + 1) z'
| otherwise
= (d, z)
err = error $ "divideByPrime: malformed arguments" ++ show (p, np, k)
quotI :: GaussianInteger -> Integer -> GaussianInteger
quotI (x :+ y) n = (x `quot` n :+ y `quot` n)
quotEvenI :: GaussianInteger -> Integer -> Maybe GaussianInteger
quotEvenI (x :+ y) n
| xr == 0
, yr == 0
= Just (xq :+ yq)
| otherwise
= Nothing
where
(xq, xr) = x `quotRem` n
(yq, yr) = y `quotRem` n
instance U.UniqueFactorisation GaussianInteger where
factorise 0 = []
factorise g = coerce $ factorise g
isPrime g = if isPrime g then Just (Prime g) else Nothing