-- |
-- Module:      Math.NumberTheory.GaussianIntegers
-- Copyright:   (c) 2016 Chris Fredrickson, Google Inc.
-- Licence:     MIT
-- Maintainer:  Chris Fredrickson <chris.p.fredrickson@gmail.com>
--
-- This module exports functions for manipulating Gaussian integers, including
-- computing their prime factorisations.
--

{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE PostfixOperators #-}
{-# 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)
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

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 :+
-- |A Gaussian integer is a+bi, where a and b are both integers.
data GaussianInteger = (:+) { GaussianInteger -> Integer
real :: !Integer, GaussianInteger -> Integer
imag :: !Integer }
    deriving (GaussianInteger -> GaussianInteger -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: GaussianInteger -> GaussianInteger -> Bool
$c/= :: GaussianInteger -> GaussianInteger -> Bool
== :: GaussianInteger -> GaussianInteger -> Bool
$c== :: GaussianInteger -> GaussianInteger -> Bool
Eq, Eq GaussianInteger
GaussianInteger -> GaussianInteger -> Bool
GaussianInteger -> GaussianInteger -> Ordering
GaussianInteger -> GaussianInteger -> GaussianInteger
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 :: GaussianInteger -> GaussianInteger -> GaussianInteger
$cmin :: GaussianInteger -> GaussianInteger -> GaussianInteger
max :: GaussianInteger -> GaussianInteger -> GaussianInteger
$cmax :: GaussianInteger -> GaussianInteger -> GaussianInteger
>= :: GaussianInteger -> GaussianInteger -> Bool
$c>= :: GaussianInteger -> GaussianInteger -> Bool
> :: GaussianInteger -> GaussianInteger -> Bool
$c> :: GaussianInteger -> GaussianInteger -> Bool
<= :: GaussianInteger -> GaussianInteger -> Bool
$c<= :: GaussianInteger -> GaussianInteger -> Bool
< :: GaussianInteger -> GaussianInteger -> Bool
$c< :: GaussianInteger -> GaussianInteger -> Bool
compare :: GaussianInteger -> GaussianInteger -> Ordering
$ccompare :: GaussianInteger -> GaussianInteger -> Ordering
Ord, forall x. Rep GaussianInteger x -> GaussianInteger
forall x. GaussianInteger -> Rep GaussianInteger x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cto :: forall x. Rep GaussianInteger x -> GaussianInteger
$cfrom :: forall x. GaussianInteger -> Rep GaussianInteger x
Generic)

instance NFData GaussianInteger

-- |The imaginary unit, where
--
-- > ι .^ 2 == -1
ι :: GaussianInteger
ι :: GaussianInteger
ι = Integer
0 Integer -> Integer -> GaussianInteger
:+ Integer
1

instance Show GaussianInteger where
    show :: GaussianInteger -> 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 GaussianInteger where
    + :: GaussianInteger -> GaussianInteger -> GaussianInteger
(+) (Integer
a :+ Integer
b) (Integer
c :+ Integer
d) = (Integer
a forall a. Num a => a -> a -> a
+ Integer
c) Integer -> Integer -> GaussianInteger
:+ (Integer
b forall a. Num a => a -> a -> a
+ Integer
d)
    * :: GaussianInteger -> GaussianInteger -> GaussianInteger
(*) (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 -> GaussianInteger
:+ (Integer
a forall a. Num a => a -> a -> a
* Integer
d forall a. Num a => a -> a -> a
+ Integer
b forall a. Num a => a -> a -> a
* Integer
c)
    abs :: GaussianInteger -> GaussianInteger
abs = forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. GaussianInteger -> (GaussianInteger, GaussianInteger)
absSignum
    negate :: GaussianInteger -> GaussianInteger
negate (Integer
a :+ Integer
b) = (-Integer
a) Integer -> Integer -> GaussianInteger
:+ (-Integer
b)
    fromInteger :: Integer -> GaussianInteger
fromInteger Integer
n = Integer
n Integer -> Integer -> GaussianInteger
:+ Integer
0
    signum :: GaussianInteger -> GaussianInteger
signum = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. GaussianInteger -> (GaussianInteger, GaussianInteger)
absSignum

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

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

absSignum :: GaussianInteger -> (GaussianInteger, GaussianInteger)
absSignum :: GaussianInteger -> (GaussianInteger, GaussianInteger)
absSignum GaussianInteger
0 = (GaussianInteger
0, GaussianInteger
0)
absSignum z :: GaussianInteger
z@(Integer
a :+ Integer
b)
  -- first quadrant: (0, inf) x [0, inf)i
  | Integer
a forall a. Ord a => a -> a -> Bool
>  Integer
0 Bool -> Bool -> Bool
&& Integer
b forall a. Ord a => a -> a -> Bool
>= Integer
0 = (GaussianInteger
z, GaussianInteger
1)
  -- second quadrant: (-inf, 0] x (0, inf)i
  | Integer
a forall a. Ord a => a -> a -> Bool
<= Integer
0 Bool -> Bool -> Bool
&& Integer
b forall a. Ord a => a -> a -> Bool
>  Integer
0 = (Integer
b Integer -> Integer -> GaussianInteger
:+ (-Integer
a), GaussianInteger
ι)
  -- third quadrant: (-inf, 0) x (-inf, 0]i
  | Integer
a forall a. Ord a => a -> a -> Bool
<  Integer
0 Bool -> Bool -> Bool
&& Integer
b forall a. Ord a => a -> a -> Bool
<= Integer
0 = (-GaussianInteger
z, -GaussianInteger
1)
  -- fourth quadrant: [0, inf) x (-inf, 0)i
  | Bool
otherwise        = ((-Integer
b) Integer -> Integer -> GaussianInteger
:+ Integer
a, -GaussianInteger
ι)

instance GcdDomain GaussianInteger

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

quotRemInt :: GaussianInteger -> Integer -> (GaussianInteger, GaussianInteger)
quotRemInt :: GaussianInteger -> Integer -> (GaussianInteger, GaussianInteger)
quotRemInt GaussianInteger
z   Integer
1  = ( GaussianInteger
z, GaussianInteger
0)
quotRemInt GaussianInteger
z (-1) = (-GaussianInteger
z, GaussianInteger
0)
quotRemInt (Integer
a :+ Integer
b) Integer
c = (Integer
qa Integer -> Integer -> GaussianInteger
:+ Integer
qb, (Integer
ra forall a. Num a => a -> a -> a
- Integer
bumpA) Integer -> Integer -> GaussianInteger
:+ (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 Gaussian integer.
conjugate :: GaussianInteger -> GaussianInteger
conjugate :: GaussianInteger -> GaussianInteger
conjugate (Integer
r :+ Integer
i) = Integer
r Integer -> Integer -> GaussianInteger
:+ (-Integer
i)

-- |The square of the magnitude of a Gaussian integer.
norm :: GaussianInteger -> Integer
norm :: GaussianInteger -> Integer
norm (Integer
x :+ Integer
y) = Integer
x forall a. Num a => a -> a -> a
* Integer
x forall a. Num a => a -> a -> a
+ Integer
y forall a. Num a => a -> a -> a
* Integer
y

-- |Compute whether a given Gaussian integer is prime.
isPrime :: GaussianInteger -> Bool
isPrime :: GaussianInteger -> Bool
isPrime g :: GaussianInteger
g@(Integer
x :+ Integer
y)
    | Integer
x forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
y forall a. Eq a => a -> a -> Bool
/= Integer
0 = forall a. Num a => a -> a
abs Integer
y forall a. Integral a => a -> a -> a
`mod` Integer
4 forall a. Eq a => a -> a -> Bool
== Integer
3 Bool -> Bool -> Bool
&& forall a. Maybe a -> Bool
isJust (forall a. UniqueFactorisation a => a -> Maybe (Prime a)
U.isPrime Integer
y)
    | Integer
y forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
x forall a. Eq a => a -> a -> Bool
/= Integer
0 = forall a. Num a => a -> a
abs Integer
x forall a. Integral a => a -> a -> a
`mod` Integer
4 forall a. Eq a => a -> a -> Bool
== Integer
3 Bool -> Bool -> Bool
&& forall a. Maybe a -> Bool
isJust (forall a. UniqueFactorisation a => a -> Maybe (Prime a)
U.isPrime Integer
x)
    | Bool
otherwise        = forall a. Maybe a -> Bool
isJust forall a b. (a -> b) -> a -> b
$ forall a. UniqueFactorisation a => a -> Maybe (Prime a)
U.isPrime forall a b. (a -> b) -> a -> b
$ GaussianInteger -> Integer
norm GaussianInteger
g

-- |An infinite list of the Gaussian primes. Uses primes in Z to exhaustively
-- generate all Gaussian primes (up to associates), in order of ascending
-- magnitude.
--
-- >>> take 10 primes
-- [Prime 1+ι,Prime 2+ι,Prime 1+2*ι,Prime 3,Prime 3+2*ι,Prime 2+3*ι,Prime 4+ι,Prime 1+4*ι,Prime 5+2*ι,Prime 2+5*ι]
primes :: Infinite (U.Prime GaussianInteger)
primes :: Infinite (Prime GaussianInteger)
primes = coerce :: forall a b. Coercible a b => a -> b
coerce forall a b. (a -> b) -> a -> b
$ (Integer
1 Integer -> Integer -> GaussianInteger
:+ 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 GaussianInteger -> Integer
norm) Infinite GaussianInteger
l Infinite GaussianInteger
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
4 forall a. Eq a => a -> a -> Bool
== Integer
3) (forall a.
(Bits a, Integral a, UniqueFactorisation a) =>
a -> Prime a
U.nextPrime Integer
3 ...)

    l :: Infinite (GaussianInteger)
    l :: Infinite GaussianInteger
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 -> GaussianInteger
:+ Integer
0) Infinite (Prime Integer)
leftPrimes

    r :: Infinite (GaussianInteger)
    r :: Infinite GaussianInteger
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 GaussianInteger
findPrime Prime Integer
p) in (Integer
x Integer -> Integer -> GaussianInteger
:+ Integer
y) forall a. a -> [a] -> NonEmpty a
:| [Integer
y Integer -> Integer -> GaussianInteger
:+ Integer
x])
        Infinite (Prime Integer)
rightPrimes

-- |Find a Gaussian integer whose norm is the given prime number
-- of form 4k + 1 using
-- <http://www.ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/S0025-5718-1972-0314745-6.pdf Hermite-Serret algorithm>.
--
-- >>> import Math.NumberTheory.Primes (nextPrime)
-- >>> findPrime (nextPrime 5)
-- Prime 2+ι
findPrime :: Prime Integer -> U.Prime GaussianInteger
findPrime :: Prime Integer -> Prime GaussianInteger
findPrime Prime Integer
p = case Integer -> Prime Integer -> [Integer]
sqrtsModPrime (-Integer
1) Prime Integer
p of
    []    -> forall a. HasCallStack => String -> a
error String
"findPrime: an argument must be prime p = 4k + 1"
    Integer
z : [Integer]
_ -> forall a. a -> Prime a
Prime forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> GaussianInteger
go (forall a. Prime a -> a
unPrime Prime Integer
p) Integer
z -- Effectively we calculate gcdG' (p :+ 0) (z :+ 1)
    where
        sqrtp :: Integer
        sqrtp :: Integer
sqrtp = forall a. Integral a => a -> a
integerSquareRoot (forall a. Prime a -> a
unPrime Prime Integer
p)

        go :: Integer -> Integer -> GaussianInteger
        go :: Integer -> Integer -> GaussianInteger
go Integer
g Integer
h
            | Integer
g forall a. Ord a => a -> a -> Bool
<= Integer
sqrtp = Integer
g Integer -> Integer -> GaussianInteger
:+ Integer
h
            | Bool
otherwise  = Integer -> Integer -> GaussianInteger
go Integer
h (Integer
g forall a. Integral a => a -> a -> a
`mod` Integer
h)

-- | Compute the prime factorisation of a Gaussian integer. This is unique up to units (+/- 1, +/- i).
-- Unit factors are not included in the result.
factorise :: GaussianInteger -> [(Prime GaussianInteger, Word)]
factorise :: GaussianInteger -> [(Prime GaussianInteger, Word)]
factorise GaussianInteger
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 GaussianInteger
-> (Prime Integer, Word)
-> (GaussianInteger, [(Prime GaussianInteger, Word)])
go GaussianInteger
g (forall a. UniqueFactorisation a => a -> [(Prime a, Word)]
U.factorise forall a b. (a -> b) -> a -> b
$ GaussianInteger -> Integer
norm GaussianInteger
g)
    where
        go :: GaussianInteger -> (Prime Integer, Word) -> (GaussianInteger, [(Prime GaussianInteger, Word)])
        go :: GaussianInteger
-> (Prime Integer, Word)
-> (GaussianInteger, [(Prime GaussianInteger, Word)])
go GaussianInteger
z (Prime Integer
2, Word
e) = (GaussianInteger -> GaussianInteger
divideByTwo GaussianInteger
z, [(forall a. a -> Prime a
Prime (Integer
1 Integer -> Integer -> GaussianInteger
:+ Integer
1), Word
e)])
        go GaussianInteger
z (Prime Integer
p, Word
e)
            | forall a. Prime a -> a
unPrime Prime Integer
p forall a. Integral a => a -> a -> a
`mod` Integer
4 forall a. Eq a => a -> a -> Bool
== Integer
3
            = let e' :: Word
e' = Word
e forall a. Euclidean a => a -> a -> a
`quot` Word
2 in (GaussianInteger
z GaussianInteger -> Integer -> GaussianInteger
`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 -> GaussianInteger
:+ Integer
0), Word
e')])
            | Bool
otherwise
            = (GaussianInteger
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 GaussianInteger
gp, Word
k), (Prime GaussianInteger
gp', Word
k')])
                where
                    gp :: Prime GaussianInteger
gp = Prime Integer -> Prime GaussianInteger
findPrime Prime Integer
p
                    (Word
k, Word
k', GaussianInteger
z') = Prime GaussianInteger
-> Integer
-> Word
-> GaussianInteger
-> (Word, Word, GaussianInteger)
divideByPrime Prime GaussianInteger
gp (forall a. Prime a -> a
unPrime Prime Integer
p) Word
e GaussianInteger
z
                    gp' :: Prime GaussianInteger
gp' = forall a. a -> Prime a
Prime (forall a. Num a => a -> a
abs (GaussianInteger -> GaussianInteger
conjugate (forall a. Prime a -> a
unPrime Prime GaussianInteger
gp)))

-- | Remove all (1:+1) factors from the argument,
-- avoiding complex division.
divideByTwo :: GaussianInteger -> GaussianInteger
divideByTwo :: GaussianInteger -> GaussianInteger
divideByTwo z :: GaussianInteger
z@(Integer
x :+ Integer
y)
    | forall a. Integral a => a -> Bool
even Integer
x, forall a. Integral a => a -> Bool
even Integer
y
    = GaussianInteger -> GaussianInteger
divideByTwo forall a b. (a -> b) -> a -> b
$ GaussianInteger
z GaussianInteger -> Integer -> GaussianInteger
`quotI` Integer
2
    | forall a. Integral a => a -> Bool
odd Integer
x, forall a. Integral a => a -> Bool
odd Integer
y
    = (Integer
x forall a. Num a => a -> a -> a
- Integer
y) forall a. Euclidean a => a -> a -> a
`quot` Integer
2 Integer -> Integer -> GaussianInteger
:+ (Integer
x forall a. Num a => a -> a -> a
+ Integer
y) forall a. Euclidean a => a -> a -> a
`quot` Integer
2
    | Bool
otherwise
    = GaussianInteger
z

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

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

        err :: GaussianInteger
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 GaussianInteger
p, Integer
np, Word
k)

quotI :: GaussianInteger -> Integer -> GaussianInteger
quotI :: GaussianInteger -> Integer -> GaussianInteger
quotI (Integer
x :+ Integer
y) Integer
n = Integer
x forall a. Euclidean a => a -> a -> a
`quot` Integer
n Integer -> Integer -> GaussianInteger
:+ Integer
y forall a. Euclidean a => a -> a -> a
`quot` Integer
n

quotEvenI :: GaussianInteger -> Integer -> Maybe GaussianInteger
quotEvenI :: GaussianInteger -> Integer -> Maybe GaussianInteger
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 -> GaussianInteger
:+ 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

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

instance U.UniqueFactorisation GaussianInteger where
  factorise :: GaussianInteger -> [(Prime GaussianInteger, Word)]
factorise GaussianInteger
0 = []
  factorise GaussianInteger
g = coerce :: forall a b. Coercible a b => a -> b
coerce forall a b. (a -> b) -> a -> b
$ GaussianInteger -> [(Prime GaussianInteger, Word)]
factorise GaussianInteger
g

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