-- Copyright (c) 2011-2015, David Amos. All rights reserved.

{-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances, OverlappingInstances #-}

-- |A module for arithmetic in quadratic number fields. A quadratic number field is a field of the form Q(sqrt d),
-- where d is a square-free integer. For example, we can perform the following calculation in Q(sqrt 2):
--
-- > (1 + sqrt 2) / (2 + sqrt 2)
--
-- It is also possible to mix different square roots in the same calculation. For example:
--
-- > (1 + sqrt 2) * (1 + sqrt 3)
--
-- Square roots of negative numbers are also permitted. For example:
--
-- > i * sqrt(-3)
module Math.NumberTheory.QuadraticField where

import Prelude hiding (sqrt, (*>) )

import Data.List as L
import Math.Core.Field
import Math.Core.Utils (powersetdfs)
import Math.Algebras.VectorSpace
import Math.Algebras.TensorProduct
import Math.Algebras.Structures
import Math.NumberTheory.Factor

import Math.Algebra.LinearAlgebra hiding (inverse, (*>) )
import Math.CommutativeAlgebra.Polynomial


-- Q(sqrt n)

-- |A basis for quadratic number fields Q(sqrt d), where d is a square-free integer.
data QNFBasis = One | Sqrt Integer deriving (QNFBasis -> QNFBasis -> Bool
(QNFBasis -> QNFBasis -> Bool)
-> (QNFBasis -> QNFBasis -> Bool) -> Eq QNFBasis
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: QNFBasis -> QNFBasis -> Bool
$c/= :: QNFBasis -> QNFBasis -> Bool
== :: QNFBasis -> QNFBasis -> Bool
$c== :: QNFBasis -> QNFBasis -> Bool
Eq,Eq QNFBasis
Eq QNFBasis =>
(QNFBasis -> QNFBasis -> Ordering)
-> (QNFBasis -> QNFBasis -> Bool)
-> (QNFBasis -> QNFBasis -> Bool)
-> (QNFBasis -> QNFBasis -> Bool)
-> (QNFBasis -> QNFBasis -> Bool)
-> (QNFBasis -> QNFBasis -> QNFBasis)
-> (QNFBasis -> QNFBasis -> QNFBasis)
-> Ord QNFBasis
QNFBasis -> QNFBasis -> Bool
QNFBasis -> QNFBasis -> Ordering
QNFBasis -> QNFBasis -> QNFBasis
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 :: QNFBasis -> QNFBasis -> QNFBasis
$cmin :: QNFBasis -> QNFBasis -> QNFBasis
max :: QNFBasis -> QNFBasis -> QNFBasis
$cmax :: QNFBasis -> QNFBasis -> QNFBasis
>= :: QNFBasis -> QNFBasis -> Bool
$c>= :: QNFBasis -> QNFBasis -> Bool
> :: QNFBasis -> QNFBasis -> Bool
$c> :: QNFBasis -> QNFBasis -> Bool
<= :: QNFBasis -> QNFBasis -> Bool
$c<= :: QNFBasis -> QNFBasis -> Bool
< :: QNFBasis -> QNFBasis -> Bool
$c< :: QNFBasis -> QNFBasis -> Bool
compare :: QNFBasis -> QNFBasis -> Ordering
$ccompare :: QNFBasis -> QNFBasis -> Ordering
$cp1Ord :: Eq QNFBasis
Ord)

instance Show QNFBasis where
    show :: QNFBasis -> String
show One = "1"
    show (Sqrt d :: Integer
d) | Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== -1 = "i"
                  | Bool
otherwise = "sqrt(" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
d String -> ShowS
forall a. [a] -> [a] -> [a]
++ ")"

-- |The type for elements of quadratic number fields
type QNF = Vect Q QNFBasis

-- |Although this has the same name as the Prelude.sqrt function, it should be thought of as more like a constructor
-- for creating elements of quadratic fields.
--
-- Note that for d positive, sqrt d means the positive square root, and sqrt (-d) should be interpreted as the square root
-- with positive imaginary part, that is i * sqrt d. This has the consequence that for example, sqrt (-2) * sqrt (-3) = - sqrt 6.
sqrt :: Integer -> QNF
sqrt :: Integer -> QNF
sqrt d :: Integer
d | Integer
fr Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1   = Integer -> QNF
forall a. Num a => Integer -> a
fromInteger Integer
sq
       | Bool
otherwise = Integer -> QNF
forall a. Num a => Integer -> a
fromInteger Integer
sq QNF -> QNF -> QNF
forall a. Num a => a -> a -> a
* QNFBasis -> QNF
forall (m :: * -> *) a. Monad m => a -> m a
return (Integer -> QNFBasis
Sqrt Integer
fr)
    where (sq :: Integer
sq,fr :: Integer
fr) = Integer -> Integer -> [Integer] -> (Integer, Integer)
forall b. (Eq b, Num b) => b -> b -> [b] -> (b, b)
squaredFree 1 1 (Integer -> [Integer]
pfactors Integer
d)
          squaredFree :: b -> b -> [b] -> (b, b)
squaredFree squared :: b
squared free :: b
free (d1 :: b
d1:d2 :: b
d2:ds :: [b]
ds) =
              if b
d1 b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== b
d2 then b -> b -> [b] -> (b, b)
squaredFree (b
d1b -> b -> b
forall a. Num a => a -> a -> a
*b
squared) b
free [b]
ds else b -> b -> [b] -> (b, b)
squaredFree b
squared (b
d1b -> b -> b
forall a. Num a => a -> a -> a
*b
free) (b
d2b -> [b] -> [b]
forall a. a -> [a] -> [a]
:[b]
ds)
          squaredFree squared :: b
squared free :: b
free ds :: [b]
ds = (b
squared, b
free b -> b -> b
forall a. Num a => a -> a -> a
* [b] -> b
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [b]
ds)

sqrt2 :: QNF
sqrt2 = Integer -> QNF
sqrt 2
sqrt3 :: QNF
sqrt3 = Integer -> QNF
sqrt 3
sqrt5 :: QNF
sqrt5 = Integer -> QNF
sqrt 5
sqrt6 :: QNF
sqrt6 = Integer -> QNF
sqrt 6
sqrt7 :: QNF
sqrt7 = Integer -> QNF
sqrt 7

i :: QNF
i :: QNF
i = Integer -> QNF
sqrt (-1)

instance (Eq k, Num k) => Algebra k QNFBasis where
    unit :: k -> Vect k QNFBasis
unit x :: k
x = k
x k -> Vect k QNFBasis -> Vect k QNFBasis
forall k b. (Eq k, Num k) => k -> Vect k b -> Vect k b
*> QNFBasis -> Vect k QNFBasis
forall (m :: * -> *) a. Monad m => a -> m a
return QNFBasis
One
    mult :: Vect k (Tensor QNFBasis QNFBasis) -> Vect k QNFBasis
mult = (Tensor QNFBasis QNFBasis -> Vect k QNFBasis)
-> Vect k (Tensor QNFBasis QNFBasis) -> Vect k QNFBasis
forall k b a.
(Eq k, Num k, Ord b) =>
(a -> Vect k b) -> Vect k a -> Vect k b
linear Tensor QNFBasis QNFBasis -> Vect k QNFBasis
forall k.
(Num k, Eq k) =>
Tensor QNFBasis QNFBasis -> Vect k QNFBasis
mult'
         where mult' :: Tensor QNFBasis QNFBasis -> Vect k QNFBasis
mult' (One,x :: QNFBasis
x) = QNFBasis -> Vect k QNFBasis
forall (m :: * -> *) a. Monad m => a -> m a
return QNFBasis
x
               mult' (x :: QNFBasis
x,One) = QNFBasis -> Vect k QNFBasis
forall (m :: * -> *) a. Monad m => a -> m a
return QNFBasis
x
               mult' (Sqrt m :: Integer
m, Sqrt n :: Integer
n) | Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
n = k -> Vect k QNFBasis
forall k b. Algebra k b => k -> Vect k b
unit (Integer -> k
forall a. Num a => Integer -> a
fromInteger Integer
m)
                                      | Bool
otherwise = let (i :: Integer
i,d :: Integer
d) = [Integer] -> [Integer] -> Integer -> Integer -> (Integer, Integer)
forall b. (Ord b, Num b) => [b] -> [b] -> b -> b -> (b, b)
interdiff (Integer -> [Integer]
pfactors Integer
m) (Integer -> [Integer]
pfactors Integer
n) 1 1
                                                    in Integer -> k
forall a. Num a => Integer -> a
fromInteger Integer
i k -> Vect k QNFBasis -> Vect k QNFBasis
forall k b. (Eq k, Num k) => k -> Vect k b -> Vect k b
*> QNFBasis -> Vect k QNFBasis
forall (m :: * -> *) a. Monad m => a -> m a
return (Integer -> QNFBasis
Sqrt Integer
d)
               -- if squarefree a == product ps, b == product qs
               -- then sqrt a * sqrt b = product (intersect ps qs) * sqrt (product (symdiff ps qs))
               -- the following calculates these two products
               -- in particular, it correctly handles the case that either or both contain -1
               interdiff :: [b] -> [b] -> b -> b -> (b, b)
interdiff (p :: b
p:ps :: [b]
ps) (q :: b
q:qs :: [b]
qs) i :: b
i d :: b
d =
                   case b -> b -> Ordering
forall a. Ord a => a -> a -> Ordering
compare b
p b
q of
                   LT -> [b] -> [b] -> b -> b -> (b, b)
interdiff [b]
ps (b
qb -> [b] -> [b]
forall a. a -> [a] -> [a]
:[b]
qs) b
i (b
db -> b -> b
forall a. Num a => a -> a -> a
*b
p)
                   EQ -> [b] -> [b] -> b -> b -> (b, b)
interdiff [b]
ps [b]
qs (b
ib -> b -> b
forall a. Num a => a -> a -> a
*b
p) b
d
                   GT -> [b] -> [b] -> b -> b -> (b, b)
interdiff (b
pb -> [b] -> [b]
forall a. a -> [a] -> [a]
:[b]
ps) [b]
qs b
i (b
db -> b -> b
forall a. Num a => a -> a -> a
*b
q)
               interdiff ps :: [b]
ps qs :: [b]
qs i :: b
i d :: b
d = (b
i, b
d b -> b -> b
forall a. Num a => a -> a -> a
* [b] -> b
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ([b]
ps [b] -> [b] -> [b]
forall a. [a] -> [a] -> [a]
++ [b]
qs))

{-
instance HasConjugation Q QNFBasis where
    conj = (>>= conj') where
        conj' One = return One
        conj' sqrt_d = -1 *> return sqrt_d
    -- ie conj = linear conj', but avoiding unnecessary nf call
    sqnorm x = coeff One (x * conj x)
-}

newtype XVar = X Int deriving (XVar -> XVar -> Bool
(XVar -> XVar -> Bool) -> (XVar -> XVar -> Bool) -> Eq XVar
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: XVar -> XVar -> Bool
$c/= :: XVar -> XVar -> Bool
== :: XVar -> XVar -> Bool
$c== :: XVar -> XVar -> Bool
Eq, Eq XVar
Eq XVar =>
(XVar -> XVar -> Ordering)
-> (XVar -> XVar -> Bool)
-> (XVar -> XVar -> Bool)
-> (XVar -> XVar -> Bool)
-> (XVar -> XVar -> Bool)
-> (XVar -> XVar -> XVar)
-> (XVar -> XVar -> XVar)
-> Ord XVar
XVar -> XVar -> Bool
XVar -> XVar -> Ordering
XVar -> XVar -> XVar
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 :: XVar -> XVar -> XVar
$cmin :: XVar -> XVar -> XVar
max :: XVar -> XVar -> XVar
$cmax :: XVar -> XVar -> XVar
>= :: XVar -> XVar -> Bool
$c>= :: XVar -> XVar -> Bool
> :: XVar -> XVar -> Bool
$c> :: XVar -> XVar -> Bool
<= :: XVar -> XVar -> Bool
$c<= :: XVar -> XVar -> Bool
< :: XVar -> XVar -> Bool
$c< :: XVar -> XVar -> Bool
compare :: XVar -> XVar -> Ordering
$ccompare :: XVar -> XVar -> Ordering
$cp1Ord :: Eq XVar
Ord, Int -> XVar -> ShowS
[XVar] -> ShowS
XVar -> String
(Int -> XVar -> ShowS)
-> (XVar -> String) -> ([XVar] -> ShowS) -> Show XVar
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [XVar] -> ShowS
$cshowList :: [XVar] -> ShowS
show :: XVar -> String
$cshow :: XVar -> String
showsPrec :: Int -> XVar -> ShowS
$cshowsPrec :: Int -> XVar -> ShowS
Show)

instance Fractional QNF where
    recip :: QNF -> QNF
recip x :: QNF
x@(V ts :: [(QNFBasis, Q)]
ts) =
        let ds :: [Integer]
ds = [Integer
d | (Sqrt d :: Integer
d, _) <- QNF -> [(QNFBasis, Q)]
forall k b. Vect k b -> [(b, k)]
terms QNF
x]
            fs :: [Integer]
fs = (if (Integer -> Bool) -> [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<0) [Integer]
ds then [-1] else []) [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ Integer -> [Integer]
pfactors ((Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
lcm 1 [Integer]
ds) -- lcm is always positive
            rs :: [QNFBasis]
rs = (Integer -> QNFBasis) -> [Integer] -> [QNFBasis]
forall a b. (a -> b) -> [a] -> [b]
map (\d :: Integer
d -> case Integer
d of 1 -> QNFBasis
One; d' :: Integer
d' -> Integer -> QNFBasis
Sqrt Integer
d') ([Integer] -> [QNFBasis]) -> [Integer] -> [QNFBasis]
forall a b. (a -> b) -> a -> b
$
                 ([Integer] -> Integer) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ([[Integer]] -> [Integer]) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> a -> b
$ [Integer] -> [[Integer]]
forall a. [a] -> [[a]]
powersetdfs ([Integer] -> [[Integer]]) -> [Integer] -> [[Integer]]
forall a b. (a -> b) -> a -> b
$ [Integer]
fs 
            -- for example, for x == sqrt2 + sqrt3, we would have rs == [One, Sqrt 2, Sqrt 3, Sqrt 6]
            n :: Int
n = [QNFBasis] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [QNFBasis]
rs
            y :: Vect (GlexPoly Q XVar) QNFBasis
y = [(QNFBasis, GlexPoly Q XVar)] -> Vect (GlexPoly Q XVar) QNFBasis
forall k b. [(b, k)] -> Vect k b
V ([(QNFBasis, GlexPoly Q XVar)] -> Vect (GlexPoly Q XVar) QNFBasis)
-> [(QNFBasis, GlexPoly Q XVar)] -> Vect (GlexPoly Q XVar) QNFBasis
forall a b. (a -> b) -> a -> b
$ [QNFBasis] -> [GlexPoly Q XVar] -> [(QNFBasis, GlexPoly Q XVar)]
forall a b. [a] -> [b] -> [(a, b)]
zip [QNFBasis]
rs ([GlexPoly Q XVar] -> [(QNFBasis, GlexPoly Q XVar)])
-> [GlexPoly Q XVar] -> [(QNFBasis, GlexPoly Q XVar)]
forall a b. (a -> b) -> a -> b
$ (Int -> GlexPoly Q XVar) -> [Int] -> [GlexPoly Q XVar]
forall a b. (a -> b) -> [a] -> [b]
map (XVar -> GlexPoly Q XVar
forall v. v -> GlexPoly Q v
glexvar (XVar -> GlexPoly Q XVar)
-> (Int -> XVar) -> Int -> GlexPoly Q XVar
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> XVar
X) [1..Int
n] -- x1*1+x2*r2+...+xn*rn
            x' :: Vect (GlexPoly Q XVar) QNFBasis
x' = [(QNFBasis, GlexPoly Q XVar)] -> Vect (GlexPoly Q XVar) QNFBasis
forall k b. [(b, k)] -> Vect k b
V ([(QNFBasis, GlexPoly Q XVar)] -> Vect (GlexPoly Q XVar) QNFBasis)
-> [(QNFBasis, GlexPoly Q XVar)] -> Vect (GlexPoly Q XVar) QNFBasis
forall a b. (a -> b) -> a -> b
$ ((QNFBasis, Q) -> (QNFBasis, GlexPoly Q XVar))
-> [(QNFBasis, Q)] -> [(QNFBasis, GlexPoly Q XVar)]
forall a b. (a -> b) -> [a] -> [b]
map (\(s :: QNFBasis
s,c :: Q
c) -> (QNFBasis
s, Q -> GlexPoly Q XVar
forall k b. Algebra k b => k -> Vect k b
unit Q
c)) [(QNFBasis, Q)]
ts -- lift the coefficients in x into the polynomial algebra
            one :: Vect (GlexPoly Q XVar) QNFBasis
one = Vect (GlexPoly Q XVar) QNFBasis
x' Vect (GlexPoly Q XVar) QNFBasis
-> Vect (GlexPoly Q XVar) QNFBasis
-> Vect (GlexPoly Q XVar) QNFBasis
forall a. Num a => a -> a -> a
* Vect (GlexPoly Q XVar) QNFBasis
y
            m :: [[Q]]
m = [ [Glex XVar -> GlexPoly Q XVar -> Q
forall k b. (Num k, Eq b) => b -> Vect k b -> k
coeff (XVar -> Glex XVar
forall (m :: * -> *) v. MonomialConstructor m => v -> m v
mvar (Int -> XVar
X Int
j)) GlexPoly Q XVar
c | Int
j <- [1..Int
n]] | QNFBasis
i <- [QNFBasis]
rs, let c :: GlexPoly Q XVar
c = QNFBasis -> Vect (GlexPoly Q XVar) QNFBasis -> GlexPoly Q XVar
forall k b. (Num k, Eq b) => b -> Vect k b -> k
coeff QNFBasis
i Vect (GlexPoly Q XVar) QNFBasis
one] -- matrix of the linear system
            b :: [Q]
b = 1 Q -> [Q] -> [Q]
forall a. a -> [a] -> [a]
: Int -> Q -> [Q]
forall a. Int -> a -> [a]
replicate (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) 0
        in case [[Q]] -> [Q] -> Maybe [Q]
forall a. (Eq a, Fractional a) => [[a]] -> [a] -> Maybe [a]
solveLinearSystem [[Q]]
m [Q]
b of -- find v such that m v == b - ie find the values of x1, x2, ... xn
            Just v :: [Q]
v -> QNF -> QNF
forall k b. (Eq k, Num k, Ord b) => Vect k b -> Vect k b
nf (QNF -> QNF) -> QNF -> QNF
forall a b. (a -> b) -> a -> b
$ [(QNFBasis, Q)] -> QNF
forall k b. [(b, k)] -> Vect k b
V ([(QNFBasis, Q)] -> QNF) -> [(QNFBasis, Q)] -> QNF
forall a b. (a -> b) -> a -> b
$ [QNFBasis] -> [Q] -> [(QNFBasis, Q)]
forall a b. [a] -> [b] -> [(a, b)]
zip [QNFBasis]
rs [Q]
v
            Nothing -> String -> QNF
forall a. HasCallStack => String -> a
error "QNF.recip 0"
    fromRational :: Rational -> QNF
fromRational q :: Rational
q = Rational -> Q
forall a. Fractional a => Rational -> a
fromRational Rational
q Q -> QNF -> QNF
forall k b. (Eq k, Num k) => k -> Vect k b -> Vect k b
*> 1