    -- 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