{-|
Module      : Clif.Basis
Copyright   : (c) Matti A. Eskelinen, 2016-2017
License     : MIT
Maintainer  : matti.a.eskelinen@gmail.com
Stability   : experimental
Portability : POSIX

= Motivation
A Clifford algebra is generated by a set of linearly independent (but not necessarily orthogonal) basis vectors and a metric (bilinear form) from the vectors to a field (or unital ring). After the basis and bilinear form are fixed, the Clifford algebra is exactly the quotient of the free algebra of the basis vectors by the equality relation generated by the form (see <https://en.wikipedia.org/wiki/Clifford_algebra>).

= Implementation
Defining a basis is abstracted using the two-parameter typeclass 'Basis'. A 'Basis' /b a/ for a Clifford algebra is defined by implementing a 'metric' that takes two vectors of type /b/ (for /basis/) and returns a number of type /a/, which is then used by 'basisMul' to implement multiplication of vectors. For computational reasons, any vector type must have some canonical ordering (i.e. be an instance of 'Ord'). The constructors for 'Euclidean' and 'Lorentzian' types can be used to map any such type with a corresponding metric.

While not strictly necessary for the abstraction, functions 'canonical' and 'basisMul' are included in the class with respective default implementations ('orderOrthoBasis' and 'orthoMul') in order to provide an easy way to override them using a more efficient implementation for a user-supplied basis. Note that the default implementations are only valid for a diagonal metric; Use 'orderBasis' and 'nonOrthoMul' if you define a non-diagonal metric.

-}
{-# LANGUAGE 
    MultiParamTypeClasses
  , FlexibleInstances
  , DeriveGeneric 
  #-}
module Clif.Basis
    ( -- * Classes
      Basis(..)
      -- * Basis construction
    , Euclidean(..)
    , Lorentzian(..)
      -- * Multiplication of blades
    , orthoMul
    , nonOrthoMul
    , freeMul
      -- * Computation of canonical order 
    , orderOrthoBasis
    , orderBasis
    ) where
import GHC.Generics 


-- | A typeclass for specifying a metric space
-- The default implementation uses 'orthoMul'; Replace it if you require a non-diagonal metric.
class (Ord b, Num a) => Basis b a where
    metric :: b -> b -> a 
    canonical :: ([b], a) -> [([b], a)]
    canonical (b, a) = [orderOrthoBasis [] b a] 
    basisMul :: ([b], a) -> ([b], a) -> [([b], a)]
    basisMul a b = [orthoMul a b]

-- | Data type for constructing Euclidean basis vectors (metric is 1 for matching vectors, 0 otherwise).
newtype Euclidean a = E {getE :: a}
    deriving (Eq, Ord, Generic)

instance Show a => Show (Euclidean a) where
    show (E a) = "E " ++ show a

instance (Ord b, Num a) => Basis (Euclidean b) a where
    metric a b = if a == b then 1 else 0

-- | Data type for Lorentzian (mixed signature) basis vectors.
data Lorentzian a = T a -- ^ Timelike basis vector (-1)
                  | S a -- ^ Spacelike basis vector (+1)
    deriving (Eq, Ord, Generic)

instance Show a => Show (Lorentzian a) where
    show (T a) = "T " ++ show a
    show (S a) = "S " ++ show a

instance (Ord b, Num a) => Basis (Lorentzian b) a where
    metric (T a) (T b) = if a == b then -1 else 0
    metric (S a) (S b) = if a == b then  1 else 0
    metric _ _         = 0

-- | Free multiplication of scaled blades (concatenation of vectors and multiplication of scalars).
freeMul :: Num a => ([b], a) -> ([b], a) -> ([b], a)
freeMul (s, m) (t, n) = (s++t, m * n) 

-- | Multiplies together two blades of orthogonal vectors and simplifies the results to canonical order using 'orderOrthoBasis'.
orthoMul :: Basis b a => ([b], a) -> ([b], a) -> ([b], a)
orthoMul = (uncurry (orderOrthoBasis []) .) . freeMul

-- | Multiplies two blades of (possibly non-orthogonal) basis vectors and simplifies the result to canonical order using 'orderBasis'.
-- Note that this generally results in a direct sum of blades, i.e. a list.
nonOrthoMul :: (Eq a, Basis b a) => ([b], a) -> ([b], a) -> [([b], a)]
nonOrthoMul = (uncurry (orderBasis []) .) . freeMul

-- | Canonical (ordered) form of a scaled basis blade of orthogonal vectors
-- Uses a gnome sort to commute the basis vectors to order and keeps track of the commutations 
-- w.r.t. the blade multiplier, applying the equation 
--
-- /ba = -ab/
--
-- for each consecutive pair /ba/ in wrong order.
-- Note that this relation holds only for a basis with a diagonal bilinear form, for others use the more general (albeit slower) 'orderBasis'.
orderOrthoBasis, backOrderOrtho :: Basis b a => [b] -> [b] -> a -> ([b], a)
orderOrthoBasis pre (a:b:rest) c

    -- Vectors are in order: move on.
    | a < b  = orderOrthoBasis (pre ++ [a]) (b:rest) c

    -- Same vector: simplify using the metric signature for the vector
    | a == b = backOrderOrtho pre rest $ metric a b * c

    -- Wrong order: flip and back up one step. Since the vectors are
    -- orthogonal, we flip the sign of the geometric product.
    | a > b  = backOrderOrtho pre (b:a:rest) (-c)
orderOrthoBasis pre rest c = (pre ++ rest, c)

-- Backwards gnome sort for orthogonal basis
backOrderOrtho []  rest c = orderOrthoBasis []         rest            c
backOrderOrtho pre rest c = orderOrthoBasis (init pre) (last pre:rest) c


-- | Canonical (ordered) form of a scaled basis blade of vectors
-- Uses a gnome sort to commute the basis vectors to order and keeps track of the commutations 
-- applying the equation 
--
-- /ab = 2 B(a,b) - ba/
--
-- for any given bilinear form B between the basis vectors and the number ring.
-- Note that this is slower than using 'orderOrthoBasis' if your bilinear form is diagonal (/B(a,b) = 0/ for /a != b/)
orderBasis, backOrder :: (Eq a, Basis b a) => [([b],a)] -> [b] -> a -> [([b],a)]
orderBasis _ _ 0 = [] -- Optimize zero to empty list
orderBasis pre (a:b:rest) c = case compare a b of
    -- When in order already, add the first element to each sorted blade
    LT -> let add = ([a],1) in orderBasis (if null pre then [add] else map (`freeMul` add) pre) (b:rest) c
    -- Simplify equals away and propagate the sort back
    EQ -> backOrder pre rest $ metric a b * c
    -- When in wrong order, apply the anticommutator relation and branch out
    GT -> backOrder pre (b:a:rest) (-c) ++ backOrder pre rest ((metric a b + metric a b) * c)
orderBasis []  rest m = [(rest, m)]
orderBasis pre rest m = map f pre
    where f (y, n) = (y ++ rest, n * m)

-- Gnome sort backwards for non-orthogonal basis
backOrder _ _ 0 = []
backOrder [] rest c = orderBasis [] rest c 
backOrder pre rest c = do
  ((prev, s), xs) <- selections pre
  case null prev of 
      True  -> orderBasis xs rest (s * c)
      False -> orderBasis ((init prev, s):xs) (last prev:rest) c

-- Given a list, returns all possible tuples containing an element of the original list and the rest of the list.
selections :: [a] -> [(a, [a])]
selections [] = []
selections (x:xs) = (x,xs):[(y,x:ys) | (y,ys)<-selections xs]