content: {-# LANGUAGE 
      MultiParamTypeClasses, 
      GeneralizedNewtypeDeriving, 
      BangPatterns,
      TypeOperators,
      DeriveDataTypeable,
      FlexibleInstances,
      TypeFamilies,
      PatternGuards,
      UndecidableInstances,
      ScopedTypeVariables #-}
  
  module Numeric.Coalgebra.Geometric
    ( 
    -- * Geometric coalgebra primitives
      BasisCoblade(..)
    , Comultivector
    -- * Operations over an eigenbasis
    , Eigenbasis(..)
    , Eigenmetric(..)
    , Euclidean(..)
    -- * Grade
    , grade
    , filterGrade
    -- * Inversions
    , reverse
    , gradeInversion
    , cliffordConjugate
    -- * Products
    , geometric
    , outer
    -- * Inner products
    , contractL
    , contractR
    , hestenes
    , dot
    , liftProduct
    ) where
  
  import Control.Monad (mfilter)
  import Data.Bits
  import Data.Word
  import Data.Data
  import Data.Ix
  import Data.Array.Unboxed
  import Numeric.Algebra
  import Prelude hiding ((-),(*),(+),negate,reverse)
  
  -- a basis vector for a simple geometric coalgebra with the Euclidean inner product
  newtype BasisCoblade m = BasisCoblade { runBasisCoblade :: Word64 } deriving 
    ( Eq,Ord,Num,Bits,Enum,Ix,Bounded,Show,Read,Real,Integral
    , Additive,Abelian,LeftModule Natural,RightModule Natural,Monoidal
    , Multiplicative,Unital,Commutative
    , Semiring,Rig
    , DecidableZero,DecidableAssociates,DecidableUnits
    )
  
  -- A metric space over an eigenbasis
  class Eigenbasis m where
    euclidean     :: proxy m -> Bool
    antiEuclidean :: proxy m -> Bool
    v             :: m -> BasisCoblade m
    e             :: Int -> m
  
  -- assuming n /= 0, find the index of the least significant set bit in a basis blade
  lsb :: BasisCoblade m -> Int
  lsb n = fromIntegral $ ix ! shiftR ((n .&. (-n)) * 0x07EDD5E59A4E28C2) 58
    where 
      -- a 64 bit deBruijn multiplication table
      ix :: UArray (BasisCoblade m) Word8
      ix = listArray (0, 63)
        [ 63,  0, 58,  1, 59, 47, 53,  2
        , 60, 39, 48, 27, 54, 33, 42,  3
        , 61, 51, 37, 40, 49, 18, 28, 20
        , 55, 30, 34, 11, 43, 14, 22,  4
        , 62, 57, 46, 52, 38, 26, 32, 41
        , 50, 36, 17, 19, 29, 10, 13, 21
        , 56, 45, 25, 31, 35, 16,  9, 12
        , 44, 24, 15,  8, 23,  7,  6,  5
        ]
  
  class (Ring r, Eigenbasis m) => Eigenmetric r m where
    metric :: m -> r
  
  type Comultivector r m = Covector r (BasisCoblade m)
  
  -- Euclidean basis, we can work with basis vectors for euclidean spaces of up to 64 dimensions without 
  -- expanding the representation of our basis vectors
  newtype Euclidean = Euclidean Int deriving 
    ( Eq,Ord,Show,Read,Num,Ix,Enum,Real,Integral
    , Data,Typeable
    , Additive,LeftModule Natural,RightModule Natural,Monoidal,Abelian,LeftModule Integer,RightModule Integer,Group
    , Multiplicative,TriviallyInvolutive,InvolutiveMultiplication,InvolutiveSemiring,Unital,Commutative
    , Semiring,Rig,Ring
    )
  
  instance Eigenbasis Euclidean where
    euclidean _ = True
    antiEuclidean _ = False
    v n = shiftL 1 (fromIntegral n)
    e = fromIntegral
  
  instance Ring r => Eigenmetric r Euclidean where
    metric _ = one
  
  grade :: BasisCoblade m -> Int
  grade = fromIntegral . count 5 . count 4 . count 3 . count 2 . count 1 . count 0 where 
    count c x = (x .&. mask) + (shiftR x p .&. mask) where 
      p = shiftL 1 c
      mask = (-1) `div` (shiftL 1 p + 1)
  
  m1powTimes :: (Num n, Bits n, Group r) => n -> r -> r
  m1powTimes n r
    | (n .&. 1) == 0 = r
    | otherwise      = negate r
  
  reorder :: Group r => BasisCoblade m -> BasisCoblade m -> r -> r
  reorder a0 b = m1powTimes $ go 0 (shiftR a0 1)
    where
      go !acc 0 = acc
      go acc a = go (acc + grade (a .&. b)) (shiftR a 1)
  
  -- <A>_k
  filterGrade :: Monoidal r => BasisCoblade m -> Int -> Comultivector r m
  filterGrade b k | grade b == k = zero
                  | otherwise    = return b
  
  instance Eigenmetric r m => Coalgebra r (BasisCoblade m) where
    comult f n m = scale (n .&. m) $ reorder n m $ f $ xor n m where
      scale b
        | euclidean n = id
        | otherwise   = (go one b *)
      go :: Eigenmetric r m => r -> BasisCoblade m -> r
      go acc 0 = acc
      go acc n' | b <- lsb n'
                , m' <- metric (e b :: m)
                = go (acc*m') (clearBit n' b)
  
  instance Eigenmetric r m => CounitalCoalgebra r (BasisCoblade m) where
    counit f = f (BasisCoblade zero)
  
  -- instance Group r => InvertibleModule r BasisCoblade where
    
  -- reversion (A~) is an involution for the outer product
  reverse :: Group r => BasisCoblade m -> Comultivector r m
  reverse b = shiftR (g * (g - 1)) 1 `m1powTimes` return b where
    g = grade b
  
  cliffordConjugate :: Group r => BasisCoblade m -> Comultivector r m
  cliffordConjugate b = shiftR (g * (g + 1)) 1 `m1powTimes` return b where
    g = grade b
  
  -- A^
  gradeInversion :: Group r => BasisCoblade m -> Comultivector r m
  gradeInversion b = grade b `m1powTimes` return b
  
  geometric :: Eigenmetric r m => BasisCoblade m -> BasisCoblade m -> Comultivector r m  
  geometric = multM
  
  outer :: Eigenmetric r m => BasisCoblade m -> BasisCoblade m -> Comultivector r m
  outer m n | m .&. n == 0 = geometric m n 
            | otherwise    = zero
  
  -- A _| B
  -- grade (A _| B) = grade B - grade A
  contractL :: Eigenmetric r m => BasisCoblade m -> BasisCoblade m -> Comultivector r m 
  contractL a b 
    | ga Prelude.> gb   = zero
    | otherwise = mfilter (\r -> grade r == gb - ga) (geometric a b)
    where
      ga = grade a
      gb = grade b
  
  -- A |_ B
  -- grade (A |_ B) = grade A - grade B
  contractR :: Eigenmetric r m => BasisCoblade m -> BasisCoblade m -> Comultivector r m
  contractR a b 
    | ga Prelude.< gb   = zero
    | otherwise = mfilter (\r -> grade r == ga - gb) (geometric a b)
    where
      ga = grade a
      gb = grade b
  
  -- the modified Hestenes' product
  dot :: Eigenmetric r m => BasisCoblade m -> BasisCoblade m -> Comultivector r m
  dot a b = mfilter (\r -> grade r == abs(grade a - grade b)) (geometric a b)
  
  -- Hestenes' inner product
  -- if 0 /= grade a <= grade b then 
  -- dot a b = hestenes a b = leftContract a b
  hestenes :: Eigenmetric r m => BasisCoblade m -> BasisCoblade m -> Comultivector r m
  hestenes a b
    | ga == 0 || gb == 0 = zero
    | otherwise = mfilter (\r -> grade r == abs(ga - gb)) (geometric a b)
    where
      ga = grade a
      gb = grade b
  
  liftProduct :: (BasisCoblade m -> BasisCoblade m -> Comultivector r m) -> Comultivector r m -> Comultivector r m -> Comultivector r m
  liftProduct f ma mb = do
    a <- ma
    b <- mb
    f a b
  