module Math.NumberTheory.GaussianIntegers (
GaussianInteger((:+)),
ι,
real,
imag,
conjugate,
norm,
divModG,
divG,
modG,
quotRemG,
quotG,
remG,
(.^),
isPrime,
primes,
gcdG,
gcdG',
findPrime,
findPrime',
factorise,
) where
import qualified Math.NumberTheory.Moduli as Moduli
import qualified Math.NumberTheory.Powers as Powers
import qualified Math.NumberTheory.Primes.Factorisation as Factorisation
import qualified Math.NumberTheory.Primes.Sieve as Sieve
import qualified Math.NumberTheory.Primes.Testing as Testing
infix 6 :+
infixr 8 .^
data GaussianInteger = (:+) { real :: !Integer, imag :: !Integer } deriving (Eq)
ι :: 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 z@(a :+ b)
| a == 0 && b == 0 = z
| a > 0 && b >= 0 = z
| a <= 0 && b > 0 = b :+ (a)
| a < 0 && b <= 0 = (a) :+ (b)
| otherwise = (b) :+ a
negate (a :+ b) = (a) :+ (b)
fromInteger n = n :+ 0
signum z@(a :+ b)
| a == 0 && b == 0 = z
| otherwise = z `divG` abs z
quotRemG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)
quotRemG = divHelper quot
quotG :: GaussianInteger -> GaussianInteger -> GaussianInteger
n `quotG` d = q where (q,_) = quotRemG n d
remG :: GaussianInteger -> GaussianInteger -> GaussianInteger
n `remG` d = r where (_,r) = quotRemG n d
divModG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)
divModG = divHelper div
divG :: GaussianInteger -> GaussianInteger -> GaussianInteger
n `divG` d = q where (q,_) = divModG n d
modG :: GaussianInteger -> GaussianInteger -> GaussianInteger
n `modG` d = r where (_,r) = divModG n d
divHelper :: (Integer -> Integer -> Integer) -> GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)
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 :: 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 && Testing.isPrime y
| y == 0 && x /= 0 = abs x `mod` 4 == 3 && Testing.isPrime x
| otherwise = Testing.isPrime $ norm g
primes :: [GaussianInteger]
primes = [ g
| p <- Sieve.primes
, g <- if p `mod` 4 == 3
then [p :+ 0]
else
if p == 2
then [1 :+ 1]
else let x :+ y = findPrime' p
in [x :+ y, y :+ x]
]
gcdG :: GaussianInteger -> GaussianInteger -> GaussianInteger
gcdG g h = gcdG' (abs g) (abs h)
gcdG' :: GaussianInteger -> GaussianInteger -> GaussianInteger
gcdG' g h
| h == 0 = g
| otherwise = gcdG' h (abs (g `modG` h))
findPrime :: Integer -> GaussianInteger
findPrime p
| p == 2 || (p `mod` 4 == 1 && Testing.isPrime p) = findPrime' p
| otherwise = error "p must be prime, and not congruent to 3 (mod 4)"
findPrime' :: Integer -> GaussianInteger
findPrime' p =
let (Just c) = Moduli.sqrtModP (1) p
k = Powers.integerSquareRoot p
bs = [1 .. k]
asbs = map (\b' -> ((b' * c) `mod` p, b')) bs
(a, b) = head [ (a', b') | (a', b') <- asbs, a' <= k]
in a :+ b
(.^) :: (Integral a) => GaussianInteger -> a -> GaussianInteger
a .^ e
| e < 0 && norm a == 1 =
case a of
1 :+ 0 -> 1
(1) :+ 0 -> if even e then 1 else (1)
0 :+ 1 -> (0 :+ (1)) .^ (abs e `mod` 4)
_ -> (0 :+ 1) .^ (abs e `mod` 4)
| e < 0 = error "Cannot exponentiate non-unit Gaussian Int to negative power"
| a == 0 = 0
| e == 0 = 1
| even e = s * s
| otherwise = a * a .^ (e 1)
where
s = a .^ div e 2
factorise :: GaussianInteger -> [(GaussianInteger, Int)]
factorise g
| g == 0 = error "0 has no prime factorisation"
| g == 1 = []
| otherwise =
let helper :: [(Integer, Int)] -> GaussianInteger -> [(GaussianInteger, Int)] -> [(GaussianInteger, Int)]
helper [] g' fs = (if g' == 1 then [] else [(g', 1)]) ++ fs
helper ((!p, !e) : pt) g' fs
| p `mod` 4 == 3 =
let pow = div e 2
gp = fromInteger p
in helper pt (g' `divG` (gp .^ pow)) ((gp, pow) : fs)
| otherwise =
let gp = findPrime' p
(!gNext, !facs) = trialDivide g' [gp, abs $ conjugate gp] []
in helper pt gNext (facs ++ fs)
in helper (reverse . Factorisation.factorise $ norm g) g []
trialDivide :: GaussianInteger -> [GaussianInteger] -> [(GaussianInteger, Int)] -> (GaussianInteger, [(GaussianInteger, Int)])
trialDivide g [] fs = (g, fs)
trialDivide g (pf : pft) fs
| g `modG` pf == 0 =
let (cnt, g') = countEvenDivisions g pf
in trialDivide g' pft ((pf, cnt) : fs)
| otherwise = trialDivide g pft fs
countEvenDivisions :: GaussianInteger -> GaussianInteger -> (Int, GaussianInteger)
countEvenDivisions g pf = helper g 0
where
helper :: GaussianInteger -> Int -> (Int, GaussianInteger)
helper g' acc
| g' `modG` pf == 0 = helper (g' `divG` pf) (1 + acc)
| otherwise = (acc, g')