{- |
module: Factor.Gaussian
description: Gaussian integers
license: MIT

maintainer: Joe Leslie-Hurd <joe@gilith.com>
stability: provisional
portability: portable
-}
module Factor.Gaussian
where

-------------------------------------------------------------------------------
-- A type of Gaussian integers
-------------------------------------------------------------------------------

data Gaussian = Gaussian Integer Integer
  deriving (Gaussian -> Gaussian -> Bool
(Gaussian -> Gaussian -> Bool)
-> (Gaussian -> Gaussian -> Bool) -> Eq Gaussian
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Gaussian -> Gaussian -> Bool
$c/= :: Gaussian -> Gaussian -> Bool
== :: Gaussian -> Gaussian -> Bool
$c== :: Gaussian -> Gaussian -> Bool
Eq,Eq Gaussian
Eq Gaussian
-> (Gaussian -> Gaussian -> Ordering)
-> (Gaussian -> Gaussian -> Bool)
-> (Gaussian -> Gaussian -> Bool)
-> (Gaussian -> Gaussian -> Bool)
-> (Gaussian -> Gaussian -> Bool)
-> (Gaussian -> Gaussian -> Gaussian)
-> (Gaussian -> Gaussian -> Gaussian)
-> Ord Gaussian
Gaussian -> Gaussian -> Bool
Gaussian -> Gaussian -> Ordering
Gaussian -> Gaussian -> Gaussian
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 :: Gaussian -> Gaussian -> Gaussian
$cmin :: Gaussian -> Gaussian -> Gaussian
max :: Gaussian -> Gaussian -> Gaussian
$cmax :: Gaussian -> Gaussian -> Gaussian
>= :: Gaussian -> Gaussian -> Bool
$c>= :: Gaussian -> Gaussian -> Bool
> :: Gaussian -> Gaussian -> Bool
$c> :: Gaussian -> Gaussian -> Bool
<= :: Gaussian -> Gaussian -> Bool
$c<= :: Gaussian -> Gaussian -> Bool
< :: Gaussian -> Gaussian -> Bool
$c< :: Gaussian -> Gaussian -> Bool
compare :: Gaussian -> Gaussian -> Ordering
$ccompare :: Gaussian -> Gaussian -> Ordering
$cp1Ord :: Eq Gaussian
Ord)

instance Show Gaussian where
  show :: Gaussian -> String
show Gaussian
x = case Gaussian
x of
             Gaussian Integer
a Integer
0 -> Integer -> String
forall a. Show a => a -> String
show Integer
a
             Gaussian Integer
0 Integer
b -> Integer -> String
forall a. (Eq a, Num a, Show a) => a -> String
showIm Integer
b
             Gaussian Integer
a Integer
b | Integer
b Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 -> Integer -> String
forall a. Show a => a -> String
show Integer
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" - " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. (Eq a, Num a, Show a) => a -> String
showIm (-Integer
b)
             Gaussian Integer
a Integer
b | Bool
otherwise -> Integer -> String
forall a. Show a => a -> String
show Integer
a String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" + " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. (Eq a, Num a, Show a) => a -> String
showIm Integer
b
    where
      showIm :: a -> String
showIm a
1 = String
"i"
      showIm (-1) = String
"-i"
      showIm a
i = a -> String
forall a. Show a => a -> String
show a
i String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"i"

real :: Integer -> Gaussian
real :: Integer -> Gaussian
real = (Integer -> Integer -> Gaussian) -> Integer -> Integer -> Gaussian
forall a b c. (a -> b -> c) -> b -> a -> c
flip Integer -> Integer -> Gaussian
Gaussian Integer
0

imaginary :: Integer -> Gaussian
imaginary :: Integer -> Gaussian
imaginary = Integer -> Integer -> Gaussian
Gaussian Integer
0

-------------------------------------------------------------------------------
-- Ring operations
-------------------------------------------------------------------------------

zero :: Gaussian
zero :: Gaussian
zero = Integer -> Gaussian
real Integer
0

one :: Gaussian
one :: Gaussian
one = Integer -> Gaussian
real Integer
1

negate :: Gaussian -> Gaussian
negate :: Gaussian -> Gaussian
negate (Gaussian Integer
x Integer
y) = Integer -> Integer -> Gaussian
Gaussian (Integer -> Integer
forall a. Num a => a -> a
Prelude.negate Integer
x) (Integer -> Integer
forall a. Num a => a -> a
Prelude.negate Integer
y)

add :: Gaussian -> Gaussian -> Gaussian
add :: Gaussian -> Gaussian -> Gaussian
add (Gaussian Integer
a Integer
b) (Gaussian Integer
c Integer
d) = Integer -> Integer -> Gaussian
Gaussian (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
c) (Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d)

subtract :: Gaussian -> Gaussian -> Gaussian
subtract :: Gaussian -> Gaussian -> Gaussian
subtract Gaussian
x Gaussian
y = Gaussian -> Gaussian -> Gaussian
add Gaussian
x (Gaussian -> Gaussian
Factor.Gaussian.negate Gaussian
y)

multiply :: Gaussian -> Gaussian -> Gaussian
multiply :: Gaussian -> Gaussian -> Gaussian
multiply (Gaussian Integer
a Integer
b) (Gaussian Integer
c Integer
d) = Integer -> Integer -> Gaussian
Gaussian (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
d) (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
d Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
c)

-------------------------------------------------------------------------------
-- Gaussian integers form a Euclidean domain, which means there exist
-- functions
--
--   norm :: Gaussian -> Integer
--   division :: Gaussian -> Gaussian -> (Gaussian,Gaussian)
--
-- satisfying the following properties:
--
--   1. norm x >= 0
--   2. norm x == 0 <=> x == 0
--   3. y /= 0 ==> norm x <= norm (x*y)
--   4. y /= 0 && (q,r) = division x y ==> x == q*y + r && norm r < norm y
--
-- In fact the division function below satisfies the stronger property
--
--   norm r <= norm y `div` 2
-------------------------------------------------------------------------------

norm :: Gaussian -> Integer
norm :: Gaussian -> Integer
norm (Gaussian Integer
a Integer
b) = Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
b

division :: Gaussian -> Gaussian -> (Gaussian,Gaussian)
division :: Gaussian -> Gaussian -> (Gaussian, Gaussian)
division Gaussian
x Gaussian
y = (Gaussian
q, Gaussian -> Gaussian -> Gaussian
Factor.Gaussian.subtract Gaussian
x (Gaussian -> Gaussian -> Gaussian
multiply Gaussian
q Gaussian
y))
  where
    q :: Gaussian
q = Integer -> Integer -> Gaussian
Gaussian (Integer -> Integer
closest (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
d)) (Integer -> Integer
closest (Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
d))
    closest :: Integer -> Integer
closest Integer
i = (Integer
i Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
n2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
n  -- returns closest integer to i/n
    n2 :: Integer
n2 = Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
2
    n :: Integer
n = Gaussian -> Integer
norm Gaussian
y
    Gaussian Integer
a Integer
b = Gaussian
x
    Gaussian Integer
c Integer
d = Gaussian
y

quotient :: Gaussian -> Gaussian -> Gaussian
quotient :: Gaussian -> Gaussian -> Gaussian
quotient Gaussian
x Gaussian
y = (Gaussian, Gaussian) -> Gaussian
forall a b. (a, b) -> a
fst ((Gaussian, Gaussian) -> Gaussian)
-> (Gaussian, Gaussian) -> Gaussian
forall a b. (a -> b) -> a -> b
$ Gaussian -> Gaussian -> (Gaussian, Gaussian)
division Gaussian
x Gaussian
y

remainder :: Gaussian -> Gaussian -> Gaussian
remainder :: Gaussian -> Gaussian -> Gaussian
remainder Gaussian
x Gaussian
y = (Gaussian, Gaussian) -> Gaussian
forall a b. (a, b) -> b
snd ((Gaussian, Gaussian) -> Gaussian)
-> (Gaussian, Gaussian) -> Gaussian
forall a b. (a -> b) -> a -> b
$ Gaussian -> Gaussian -> (Gaussian, Gaussian)
division Gaussian
x Gaussian
y

divides :: Gaussian -> Gaussian -> Bool
divides :: Gaussian -> Gaussian -> Bool
divides Gaussian
x Gaussian
y = Gaussian
y Gaussian -> Gaussian -> Bool
forall a. Eq a => a -> a -> Bool
== Gaussian
zero Bool -> Bool -> Bool
|| (Gaussian
x Gaussian -> Gaussian -> Bool
forall a. Eq a => a -> a -> Bool
/= Gaussian
zero Bool -> Bool -> Bool
&& Gaussian -> Gaussian -> Gaussian
remainder Gaussian
y Gaussian
x Gaussian -> Gaussian -> Bool
forall a. Eq a => a -> a -> Bool
== Gaussian
zero)

-------------------------------------------------------------------------------
-- Every Euclidean domain a admits the definition of a greatest common
-- divisor function
--
--   egcd :: a -> a -> (a,(a,a))
--
-- such that if (g,(s,t)) = egcd x y then:
--
--   1. divides g x
--   2. divides g y
--   3. s*x + t*y == g
-------------------------------------------------------------------------------

egcd :: Gaussian -> Gaussian -> (Gaussian,(Gaussian,Gaussian))
egcd :: Gaussian -> Gaussian -> (Gaussian, (Gaussian, Gaussian))
egcd Gaussian
x Gaussian
y | Gaussian
y Gaussian -> Gaussian -> Bool
forall a. Eq a => a -> a -> Bool
== Gaussian
zero = (Gaussian
x,(Gaussian
one,Gaussian
zero))
egcd Gaussian
x Gaussian
y | Bool
otherwise =
    (Gaussian
g, (Gaussian
t, Gaussian -> Gaussian -> Gaussian
Factor.Gaussian.subtract Gaussian
s (Gaussian -> Gaussian -> Gaussian
multiply Gaussian
q Gaussian
t)))
  where
    (Gaussian
q,Gaussian
r) = Gaussian -> Gaussian -> (Gaussian, Gaussian)
division Gaussian
x Gaussian
y
    (Gaussian
g,(Gaussian
s,Gaussian
t)) = Gaussian -> Gaussian -> (Gaussian, (Gaussian, Gaussian))
egcd Gaussian
y Gaussian
r

gcd :: Gaussian -> Gaussian -> Gaussian
gcd :: Gaussian -> Gaussian -> Gaussian
gcd Gaussian
x Gaussian
y = (Gaussian, (Gaussian, Gaussian)) -> Gaussian
forall a b. (a, b) -> a
fst ((Gaussian, (Gaussian, Gaussian)) -> Gaussian)
-> (Gaussian, (Gaussian, Gaussian)) -> Gaussian
forall a b. (a -> b) -> a -> b
$ Gaussian -> Gaussian -> (Gaussian, (Gaussian, Gaussian))
egcd Gaussian
x Gaussian
y