    {-# 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.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
      ( 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
      )

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