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

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

-- |A module defining the algebra of quaternions over an arbitrary field.
--
-- The quaternions are the algebra defined by the basis {1,i,j,k}, where i^2 = j^2 = k^2 = ijk = -1
module Math.Algebras.Quaternions where

import Prelude hiding ( (*>) )

import Math.Core.Field
import Math.Algebras.VectorSpace
import Math.Algebras.TensorProduct
import Math.Algebras.Structures


-- Conway & Smith, On Quaternions and Octonions

-- QUATERNIONS

data HBasis = One | I | J | K deriving (Eq,Ord)

type Quaternion k = Vect k HBasis

instance Show HBasis where
    show One = "1"
    show I = "i"
    show J = "j"
    show K = "k"

instance (Eq k, Num k) => Algebra k HBasis where
    unit x = x *> return One
    mult = linear mult'
         where mult' (One,b) = return b
               mult' (b,One) = return b
               mult' (I,I) = unit (-1)
               mult' (J,J) = unit (-1)
               mult' (K,K) = unit (-1)
               mult' (I,J) = return K
               mult' (J,I) = -1 *> return K
               mult' (J,K) = return I
               mult' (K,J) = -1 *> return I
               mult' (K,I) = return J
               mult' (I,K) = -1 *> return J

-- |The quaternions have {1,i,j,k} as basis, where i^2 = j^2 = k^2 = ijk = -1.
i,j,k :: Num k => Quaternion k
i = return I
j = return J
k = return K


class Algebra k a => HasConjugation k a where
    -- |A conjugation operation is required to satisfy the following laws:
    --
    -- * conj (x+y) = conj x + conj y
    --
    -- * conj (x*y) = conj y * conj x  (note the order-reversal)
    --
    -- * conj (conj x) = x
    --
    -- * conj x = x if and only if x in k
    conj :: Vect k a -> Vect k a
    -- |The squared norm is defined as sqnorm x = x * conj x. It satisfies:
    --
    -- * sqnorm (x*y) = sqnorm x * sqnorm y
    --
    -- * sqnorm (unit k) = k^2, for k a scalar
    sqnorm :: Vect k a -> k

-- |If an algebra has a conjugation operation, then it has multiplicative inverses,
-- via 1\/x = conj x \/ sqnorm x
instance (Eq k, Fractional k, Ord a, Show a, HasConjugation k a) => Fractional (Vect k a) where
    recip 0 = error "recip 0"
    recip x = (1 / sqnorm x) *> conj x
    fromRational q = fromRational q *> 1

-- |The scalar part of the quaternion w+xi+yj+zk is w. Also called the real part.
scalarPart :: (Num k) => Quaternion k -> k
scalarPart = coeff One

-- |The vector part of the quaternion w+xi+yj+zk is xi+yj+zk. Also called the pure part.
vectorPart :: (Eq k, Num k) => Quaternion k -> Quaternion k
vectorPart q = q - scalarPart q *> 1

instance (Eq k, Num k) => HasConjugation k HBasis where
    conj = (>>= conj') where
        conj' One = return One
        conj' imag = -1 *> return imag
    -- ie conj = linear conj', but avoiding unnecessary nf call
    sqnorm x = sum $ map ((^2) . snd) $ terms x
    -- sqnorm x = scalarPart (conj x * x) -- the vector part will be zero anyway
    -- sqnorm x = x <.> x
{-
instance Fractional k => Fractional (Quaternion k) where
    recip 0 = error "Quaternion.recip 0"
    recip x = (1 / sqnorm x) *> conj x
    fromRational q = fromRational q *> 1
-}

x <.> y = scalarPart (conj x * y)
-- x <..> y = 1/2 * (sqnorm x + sqnorm y - sqnorm (x-y))



x^-1 = recip x

-- Conway p40
refl q = \x -> -q * conj x * q

-- Given a linear function f on the quaternions, return the matrix representing it,
-- relative to a given basis. The matrix is considered as acting on the right.
asMatrix f bs = [ let fi = f ei in [ej <.> fi | ej <- bs] | ei <- bs  ]
-- It is possible to write this function using coeff, instead of <.>,
-- but then you have to pass in I,J,K, instead of i,j,k, which is uglier.

-- Conway p24
-- A homomorphism from H\0 to SO3
-- if q is restricted to unit quaternions, this is a double cover of SO3 (since q, -q induce same rotation)
-- The unit quaternions form the group Spin3
reprSO3' q = \x -> q^-1 * x * q

-- |Given a non-zero quaternion q in H, the map x -> q^-1 * x * q defines an action on the 3-dimensional vector space
-- of pure quaternions X (ie linear combinations of i,j,k). It turns out that this action is a rotation of X,
-- and this is a surjective group homomorphism from H* onto SO3. If we restrict q to the group of unit quaternions
-- (those of norm 1), then this homomorphism is 2-to-1 (since q and -q give the same rotation).
-- This shows that the multiplicative group of unit quaternions is isomorphic to Spin3, the double cover of SO3.
--
-- @reprSO3 q@ returns the 3*3 matrix representing this map.
reprSO3 :: (Eq k, Fractional k) => Quaternion k -> [[k]]
reprSO3 q = reprSO3' q `asMatrix` [i,j,k]
-- It's clear from the definition that repr3' q leaves scalars invariant

-- for achiral elts, ie GO3\SO3, we compose the above with conj

-- For unit quaternions, this is a double cover of SO4 (since (l,r), (-l,-r) induce same rotation)
-- Ordered pairs of unit quaternions form the group Spin4
reprSO4' (l,r) = \x -> l^-1 * x * r
-- then (l1,r1) * (l2,r2) -> (l1*l2,r1*r2)
-- having l^-1 is required for this to work

-- |Given a pair of unit quaternions (l,r), the map x -> l^-1 * x * r defines an action on the 4-dimensional space
-- of quaternions. It turns out that this action is a rotation, and this is a surjective group homomorphism
-- onto SO4. The homomorphism is 2-to-1 (since (l,r) and (-l,-r) give the same map).
-- This shows that the multiplicative group of pairs of unit quaternions (with pointwise multiplication)
-- is isomorphic to Spin4, the double cover of SO4.
--
-- @reprSO4 (l,r)@ returns the 4*4 matrix representing this map.
reprSO4 :: (Eq k, Fractional k) => (Quaternion k, Quaternion k) -> [[k]]
reprSO4 (l,r) = reprSO4' (l,r) `asMatrix` [1,i,j,k]
-- could consider checking that l,r are unit length - except that this is hard to achieve working over Q

reprSO4d lr = reprSO4 (p1 lr, p2 lr)

-- for achiral elts, GO4\SO4, we compose the above with conj


-- DUAL SPACE OF QUATERNIONS AS COALGEBRA

one',i',j',k' :: Num k => Vect k (Dual HBasis)
one' = return (Dual One)
i' = return (Dual I)
j' = return (Dual J)
k' = return (Dual K)

-- Coalgebra structure on the dual vector space to the quaternions
-- The comult is the transpose of mult
instance (Eq k, Num k) => Coalgebra k (Dual HBasis) where
    counit = unwrap . linear counit'
        where counit' (Dual One) = return ()
              counit' _ = zerov
    comult = linear comult'
        where comult' (Dual One) = return (Dual One, Dual One) <+>
                  (-1) *> ( return (Dual I, Dual I) <+> return (Dual J, Dual J) <+> return (Dual K, Dual K) )
              comult' (Dual I) = return (Dual One, Dual I) <+> return (Dual I, Dual One) <+>
                  return (Dual J, Dual K) <+> (-1) *> return (Dual K, Dual J)
              comult' (Dual J) = return (Dual One, Dual J) <+> return (Dual J, Dual One) <+>
                  return (Dual K, Dual I) <+> (-1) *> return (Dual I, Dual K)
              comult' (Dual K) = return (Dual One, Dual K) <+> return (Dual K, Dual One) <+>
                  return (Dual I, Dual J) <+> (-1) *> return (Dual J, Dual I)

{-
-- Of course, we can define this coalgebra structure on the quaternions themselves
-- However, it is not compatible with the algebra structure: we don't get a bialgebra
instance Num k => Coalgebra k HBasis where
    counit = unwrap . linear counit'
        where counit' One = return ()
              counit' _ = zero
    comult = linear comult'
        where comult' One = return (One,One) <+> (-1) *> ( return (I,I) <+> return (J,J) <+> return (K,K) )
              comult' I = return (One,I) <+> return (I,One) <+> return (J,K) <+> (-1) *> return (K,J)
              comult' J = return (One,J) <+> return (J,One) <+> return (K,I) <+> (-1) *> return (I,K)
              comult' K = return (One,K) <+> return (K,One) <+> return (I,J) <+> (-1) *> return (J,I)
-}

{-
-- Set coalgebra instance
instance Num k => Coalgebra k HBasis where
    counit (V ts) = sum [x | (m,x) <- ts] -- trace
    comult = fmap (\m -> T m m)           -- diagonal
-}

{-
instance Num k => Coalgebra k HBasis where
    counit (V ts) = sum [x | (One,x) <- ts]
    comult = linear cm
        where cm m = if m == One then return (m,m) else return (m,One) <+> return (One,m)
-}