```{-# LANGUAGE
MultiParamTypeClasses,
GeneralizedNewtypeDeriving,
BangPatterns,
TypeOperators,
DeriveDataTypeable,
FlexibleInstances,
TypeFamilies,
UndecidableInstances,
ScopedTypeVariables #-}

module Numeric.Coalgebra.Geometric
(
-- * Geometric coalgebra primitives
, Comultivector
-- * Operations over an eigenbasis
, Eigenbasis(..)
, Eigenmetric(..)
, Euclidean(..)
-- * Inversions
, reverse
, cliffordConjugate
-- * Products
, geometric
, outer
-- * Inner products
, contractL
, contractR
, hestenes
, dot
, liftProduct
) where

import Data.Bits
import Data.Functor.Representable.Trie
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
, Additive,Abelian,LeftModule Natural,RightModule Natural,Monoidal
, Multiplicative,Unital,Commutative
, Semiring,Rig
, DecidableZero,DecidableAssociates,DecidableUnits
)

instance HasTrie (BasisCoblade m) where
type BaseTrie (BasisCoblade m) = BaseTrie Word64
embedKey = embedKey . runBasisCoblade
projectKey = BasisCoblade . projectKey

-- 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
, Data,Typeable
, Additive,LeftModule Natural,RightModule Natural,Monoidal,Abelian,LeftModule Integer,RightModule Integer,Group
, Multiplicative,TriviallyInvolutive,InvolutiveMultiplication,InvolutiveSemiring,Unital,Commutative
, Semiring,Rig,Ring
)

instance HasTrie Euclidean where
type BaseTrie Euclidean = BaseTrie Int
embedKey (Euclidean i) = embedKey i
projectKey = Euclidean . projectKey

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 = 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 :: (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
```