Haskell module for Clifford Algebra tensor representations, and a Geometric Algebra module utilizing the tensor representation

Clifford.hs

-- credit goes to sigfpe
{-# LANGUAGE MultiParamTypeClasses
    ,TemplateHaskell
    ,GeneralizedNewtypeDeriving
    ,DeriveFunctor
    ,FunctionalDependencies
    ,FlexibleInstances
    ,UndecidableInstances
    ,FlexibleContexts
    ,Rank2Types #-}

{-# OPTIONS_GHC -fno-warn-name-shadowing -fno-warn-missing-methods #-}

module Clifford where

import Language.Haskell.TH

-- basic components of Clifford algebras

newtype I r = I r deriving (Eq, Ord, Read, Num, Real)
instance (Show r) => Show (I r) where
    show (I r) = ',' :show r

-- Complex numbers form a normed commutative division algebra
data Complex a = C a a deriving (Eq,Show)

-- very useful to have a functor to map over (we can do this via derive)
instance Functor Complex where
    fmap f (C a b) = C (f a) (f b)

-- 4 multiplications
instance Num a => Num (Complex a) where
    fromInteger a = C (fromInteger a) 0
    (C a b)+(C a' b') = C (a+a') (b+b')
    (C a b)-(C a' b') = C (a-a') (b-b')
    (C a b)*(C a' b') = C (a*a'-b*b') (a*b'+b*a')

data Split a = S a a deriving (Eq,Show,Functor)

-- 2 multiplications
instance Num a => Num (Split a) where
    fromInteger a = S (fromInteger a) (fromInteger a)
    (S a b)+(S a' b') = S (a+a') (b+b')
    (S a b)-(S a' b') = S (a-a') (b-b')
    (S a b)*(S a' b') = S (a*a') (b*b')

data Matrix a = M a a a a deriving (Eq,Show,Functor)

-- 8 multiplications
instance Num a => Num (Matrix a) where
    fromInteger a = M (fromInteger a) 0 0 (fromInteger a)
    (M a b c d)+(M a' b' c' d') = M (a+a') (b+b')
                                    (c+c') (d+d')
    (M a b c d)-(M a' b' c' d') = M (a-a') (b-b')
                                    (c-c') (d-d')
    (M a b c d)*(M a' b' c' d') = M (a*a'+b*c') (a*b'+b*d')
                                    (c*a'+d*c') (c*b'+d*d')

-- Quaternions form a normed non-commutative division algebra
-- Quaternions have injective homomorphisms to C(2), R(4)
data Quaternion a = Q a a a a deriving (Eq,Show,Functor)

-- 16 multiplications
instance Num a => Num (Quaternion a) where
    fromInteger a = Q (fromInteger a) 0 0 0
    (Q a b c d)+(Q a' b' c' d') = Q (a+a') (b+b')
                                    (c+c') (d+d')
    (Q a b c d)-(Q a' b' c' d') = Q (a-a') (b-b')
                                    (c-c') (d-d')
    (Q a b c d)*(Q a' b' c' d') = Q (a*a'-b*b'-c*c'-d*d')
                                    (a*b'+b*a'+c*d'-d*c')
                                    (a*c'+c*a'+d*b'-b*d')
                                    (a*d'+d*a'+b*c'-c*b')

-- explicit isomorphisms (can be read in either direction)
-- convert 8 times 2 multiplications into 2 times 8 (trivial mapping)
ms2sm :: (Num a) => Matrix (Split a) -> Split (Matrix a)
ms2sm (M (S a b) (S c d) (S f g) (S h j)) =
    S (M a c f h) (M b d g j)

-- convert 8 times 16 multiplications into 16 times 8 (trivial mapping)
mq2qm :: (Num a) => Matrix (Quaternion a) -> Quaternion (Matrix a)
mq2qm (M (Q a b c d) (Q f g h j) (Q k l m n) (Q p q r s) ) =
    Q (M a f k p)
      (M b g l q)
      (M c h m r)
      (M d j n s)

-- convert 4 times 4 multiplications into 2 times 4!!
-- C (C a b) (C c d) isomorph to 1xC a b + IxC c d
cc2sc :: (Num a) => Complex (Complex a) -> Split (Complex a)
cc2sc (C (C a b) (C c d) ) =
    S (C (a+d) (c-b))
      (C (a-d) (c+b))

-- convert 4 times 16 multiplications into 8 times 4!!
cq2mc :: (Num a) => Complex (Quaternion a) -> Matrix (Complex a)
cq2mc (C (Q a b c d) (Q f g h l)) =
    M (C (a+g) (f-b))
      (C (l+c) (-h-d))
      (C (l-c) (h-d))
      (C (a-g) (f+b))

-- convert 16 times 16 multiplications into 8 times 8!!
qq2mm :: (Num a) => Quaternion (Quaternion a) -> Matrix (Matrix a)
qq2mm (Q (Q a b c d) (Q f g h j) (Q k l m n) (Q p q r s)) =
    M (M (a+g+m+s) (b-f-n+r)  (f-b-n+r)  (a+g-m-s))
      (M (c+j-k-q) (d-h+l-p)  (h-d+l-p)  (c+j+k+q))
      (M (j-c+k-q) (d+h+l+p)  (-d-h+l+p) (j-c-k+q))
      (M (a-g+m-s) (-b-f+n+r) (b+f+n+r)  (a-g-m+s))

bott :: Num a => Quaternion (Matrix (Quaternion (Matrix a))) -> Matrix (Matrix (Matrix (Matrix a)))
bott = qq2mm . fmap mq2qm

-- some educational embeddings from splits, complex and quaternion numbers
-- to the matrix algebra over the real numbers
s2m  :: (Num a) => Split a -> Matrix a
s2m  (S a b) = M (a+b) 0 0 (a-b) -- Jordan form
c2m  :: (Num a) => Complex a -> Matrix a
c2m  (C a b) = M a (-b) b a
q2mm :: (Num a) => Quaternion a -> Matrix (Matrix a)
q2mm (Q a b c d) =
    M (M a b (-b) a)
      (M c d (-d) c)
      (M (-c) d (-d) c)
      (M a (-b) b a)

class Clifford b r | b -> r where
    one :: r -> b
    e_  :: Int -> r -> b -- enumerate grade 0 vectors
    b_  :: Int -> r -> b -- enumerate multi-vectors b0 .. b2^(p+q)-1
    gft :: [r] -> (b,[r])-- fast forward generalized fourier transform
    fmapI :: (r -> r) -> b -> b -- map over the real numbers

instance Clifford (I r) r where
    one a  = I a
    e_ _ _ = error ""
    b_ 0 a = I a
    b_ _ _ = error ""
    gft [] = error ""
    gft (a:as) = (I a, as)
    fmapI f (I r) = I (f r)

instance (Num b,Clifford b r) => Clifford (Complex b) r where
    one a = C (one a) 0
    e_ 1 a = C 0 (one a)
    b_ 0 a = C (one a) 0
    b_ n a
        | hi == 0   = C (b_ lo a) 0
        | otherwise = C 0 (b_ lo a)
        where hi = n `mod` 2
              lo = n `div` 2
    gft as = let (a, as1) = gft as
                 (b, as2) = gft as1
             in (C a b, as2)
    fmapI f (C a b) = C (fmapI f a) (fmapI f b)

instance (Num b,Clifford b r) => Clifford (Split b) r where
    one a = let b = one a in S b b
    e_ 1 a = let b = one a in S b (-b)
    b_ 0 a = let b = one a in S b b
    b_ n a
        | hi == 0   = let b = (b_ lo a) in S b b
        | otherwise = let b = (b_ lo a) in S b (-b)
        where hi = n `mod` 2
              lo = n `div` 2
    gft as = let (a, as1) = gft as
                 (b, as2) = gft as1
             in (S (a+b) (a-b), as2)
    fmapI f (S a b) = S (fmapI f a) (fmapI f b)

instance (Num b,Clifford b r) => Clifford (Quaternion b) r where
    one a = Q (one a) 0 0 0
    e_ 1 a = Q 0 (one a) 0 0
    e_ 2 a = Q 0 0 (one a) 0
    e_ n a = Q 0 0 0 (e_ (n-2) a)
    b_ 0 a = Q (one a) 0 0 0
    b_ n a
        | hi == 0   = Q (b_ lo a) 0 0 0
        | hi == 1   = Q 0 (b_ lo a) 0 0
        | hi == 2   = Q 0 0 (b_ lo a) 0
        | otherwise = Q 0 0 0 (b_ lo a)
        where hi = n `mod` 4
              lo = n `div` 4
    gft as = let (a, as1) = gft as
                 (b, as2) = gft as1
                 (c, as3) = gft as2
                 (d, as4) = gft as3
             in (Q a b c d, as4)
    fmapI f (Q a b c d) = Q (fmapI f a) (fmapI f b) (fmapI f c) (fmapI f d)

instance (Num b,Clifford b r) => Clifford (Matrix b) r where
    one a = let b = one a in M b 0 0 b
    e_ 1 a = let b = one a in M (-b) 0 0 b
    e_ 2 a = let b = one a in M 0 b (-b) 0
    e_ n a = let b = e_ (n-2) a in M 0 b b 0
    b_ 0 a = let b = one a in M b 0 0 b
    b_ n a
        | hi == 0   = let b = (b_ lo a) in M b 0 0 b
        | hi == 1   = let b = (b_ lo a) in M (-b) 0 0 b
        | hi == 2   = let b = (b_ lo a) in M 0 b (-b) 0
        | otherwise = let b = (b_ lo a) in M 0 b b 0
        where hi = n `mod` 4
              lo = n `div` 4
    gft as = let (a, as1) = gft as
                 (b, as2) = gft as1
                 (c, as3) = gft as2
                 (d, as4) = gft as3
             in (M (a-b) (c+d) (d-c) (a+b), as4)
    fmapI f (M a b c d) = M (fmapI f a) (fmapI f b) (fmapI f c) (fmapI f d)

-- a nice application of Template Haskell: declare types at will
cliffpq :: Int -> Int -> Name -> Q Type
cliffpq 0 0 k = do return $ AppT (ConT ''I) (ConT k)
cliffpq 0 1 k = [t| Complex $(cliffpq 0 0 k) |]
cliffpq 1 0 k = [t| Split $(cliffpq 0 0 k) |]

cliffpq p q k
        | p > 0 && q > 0                = [t| Matrix $(cliffpq (p-1) (q-1) k) |]
        | p >= 2 && q == 0              = [t| Matrix $(cliffpq q (p-2) k) |]
        | p == 0 && q >= 2              = [t| Quaternion $(cliffpq (q-2) p k ) |]

-- Geometric Algebra can be defined as
-- one 1.5 :: $(cliffpq 3 0 ''Float)

-- here Haskell starts to 'suck' a bit:
-- I can't do 
-- type GA = $(cliffpq 3 0 ''Float)
-- because the type declaration phase comes before the
-- code splicing phase... (and we have a vicious cycle)

GAlgebra.hs

-- Copyright 2011, Frank Kuehnel
-- Licensed under the Apache License, Version 2.0 (the "License");
-- you may not use this file except in compliance with the License.
--   You may obtain a copy of the License at

--      http://www.apache.org/licenses/LICENSE-2.0

--   Unless required by applicable law or agreed to in writing, software
--   distributed under the License is distributed on an "AS IS" BASIS,
--   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
--   See the License for the specific language governing permissions and
--   limitations under the License

-- Examples
-- e1 = 1.0 `e` [1] :: GA
-- e2 = 1.0 `e` [2] :: GA
-- e3 = 1.0 `e` [3] :: GA
-- x = e1 ^* e2
-- inner product
-- x .*. e1
-- inner product with a scalar
-- x .*. 2.0
-- full inner product
-- x .*. x
-- x .*. (e2 ^* e3)    -> 0
-- (e1*e2) * (e1*e3)   -> -e23
-- (e1*e2) ^* (e1*e3)  -> 0

{-# LANGUAGE MultiParamTypeClasses
    ,FunctionalDependencies
    ,GeneralizedNewtypeDeriving
    ,FlexibleInstances
    ,UndecidableInstances
    ,TypeFamilies #-}

module GAlgebra (GA, cl,(*#),(^*),(.*),(*.),(.*.),(!*),x,e,zero,
    invGFT, fwdGFT, grade, grades, minmaxgrd, maxgrd, mingrd,
    isScalar, isVector, isBivector, rev, invol, norm, normmax, numTrans) where

import Data.List
import qualified Data.IntMap as IM
import qualified Data.Map as Map
import Numeric
import Control.Monad (msum)
import Language.Haskell.TH
import Clifford

-- wish one could somehow parametrize GA with p q?
newtype GA = GA {cl :: $(cliffpq 4 1 ''Double)} deriving (Eq,Num)

-- We need a guide to the Geometric Algebra

-- map from e(-q)..e(-1), e(1)..e(p) to GA (b_ idx _)
gaIndexVec2BasisMap = case (nameBasis (4,1)) of
                        Just a  -> a
                        Nothing -> []

-- map from GA (b_ idx _) grade 1 vectors to e(-q)..e(-1), e(1)..e(p)
gaBasis2IndexVecMap = map (\(a,b) -> (b,a)) gaIndexVec2BasisMap

-- map from GA (b_ idx _) for all multivectors to the index set (frame), i.e. e(-2,-1,3,4)
-- Conceptually, we need to map the power set of the grade 1 vectors to the basis in GA
-- Efficiently, I use a neat trick using the list Monad's tensoring capabilities, similar to:
powerSet f = ffilterM (const [True, False]) f

-- generate a quick lookup method to map from fourier index to index frame set
gaBasis2IndexSetMap = IM.fromList $ zipWith (\(i,s) nl -> (i, (nl, s<0))) (bases ++ [(0,1)]) names
    where
        -- efficiently generate all non 0-grade multivectors
        mvecs = init $ powerSet (*) $ map (\(i,_) -> GA (b_ i 1)) gaBasis2IndexVecMap
        -- let's figure out what fourier index the multivector corresponds to
        bases = map (head . (filter ((/=0) . snd)) . invGFT . cl . head) mvecs
        -- by convention, frame set names for multivectors
        names = powerSet const $ map fst gaIndexVec2BasisMap

-- we could certainly play with graded maps!
gaIndexSet2BasisMap = Map.fromList $ map (\(k,(n,s)) -> (n,(k,s))) $ IM.assocs gaBasis2IndexSetMap

-- a set of grade filters to be applied after invGFT
gaBaseGrades =
    let baseIds = map fst $ invGFT $ cl $ GA (one 1)
    in case mapM (\bi -> gradeOf bi) baseIds of
        Nothing         -> []
        Just grades     -> grades
        where
            gradeOf i = case IM.lookup i gaBasis2IndexSetMap of
                            Nothing             -> Nothing
                            Just ([], sgn)      -> Just 0
                            Just (l@(s:ss),sgn) -> Just $ length l

instance Show GA where
    show (GA cl) = case decode cl gaBasis2IndexSetMap of
                        Nothing     -> "too few basis vectors! Wrong map?"
                        Just ncpoly -> show (NP ncpoly)
instance Read GA where
    readsPrec _ p = case readSigned readFloat p of
                        []          -> getBasis p  1
                        [(a,pn)]    -> getBasis pn a
                    where getBasis p val =
                            case readBasis p of
                                []              -> []
                                [(E [], p)]     -> [(GA (b_ 0 val),p)]
                                [(E (b:bs), p)] -> [(val `e` (b:bs),p)]

instance Fractional GA where
    fromRational r = GA (one $ fromRational r)
    -- general inverse of a multivector may not exist in some cases
    recip g
        | (0,0) == (nm,nx)  = GA (one $ recip $ scalar g)
        | 1 == (nx)         = let ig = invol g -- unique inverse may not exist on a null cone
                              in (recip $ scalar $ g*ig) !* ig
        | isScalar (g*rg)   = let normf = recip $ scalar $ g*rg -- strictly works for versors
                              in normf !* rg
        | isScalar (g*cjg)  = let normf = recip $ scalar $ g*cjg -- sometimes referred to as Lounesto inverse
                              in normf !* rg
        | otherwise         = g -- nonsense result!
        where (nm,nx)   = minmaxgrd g
              rg        = rev g
              cjg       = invol rg  -- try clifford conjugation

-- very faintly inspired by Conal Elliott's vector space
class MultiVectorSpace mv where
    zero    :: mv
    (!*)    :: Double -> mv -> mv    -- multiply scalar with a multivector
    (^*)    :: mv -> mv -> mv       -- multivector exterior product
    (.*.)   :: mv -> mv -> mv       -- multivector symmetric inner product
    (.*)    :: mv -> mv -> mv       -- multivector left contraction
    (*.)    :: mv -> mv -> mv       -- multivector right contraction
    (*#)    :: mv -> mv -> Double    -- multivector scalar product (return grade 0 == scalar)
    x       :: mv -> mv -> mv       -- commutator product
    maxgrd  :: mv -> Int            -- maximum grade of multivector
    mingrd  :: mv -> Int            -- minimum grade of multivector (>=0)
    minmaxgrd :: mv -> (Int, Int)
    grades  :: mv -> [Int]          -- list of grades in multivector
    grade   :: Int -> mv -> mv      -- select grade of multivector
    gradeiGFT :: Int -> [(Int,Double)] -> mv
    scalar  :: mv -> Double          -- extract scalar from multivector
    fwdGFT  :: [(Int,Double)] -> mv  -- forward fast fourier transform
    rev     :: mv -> mv             -- multivector reversion
    invol   :: mv -> mv             -- multivector involution
    e       :: Double -> [Int] -> mv -- convert from index set notation
    -- normed vector space (various norms)
    norm    :: Int -> mv -> Double
    normmax :: mv -> Double
    -- some utility functions
    isScalar    :: mv -> Bool
    isVector    :: mv -> Bool
    isBivector  :: mv -> Bool

instance MultiVectorSpace GA where
    zero            = GA (one 0)
    a !* v          = GA (fmapI (*a) $ cl v)         -- believe this is an efficient operation

    v ^*  w         = foldl' (+) zero outer
                        where (iv,iw) = (invGFT $ cl v, invGFT $ cl w)
                              (gv,gw) = (gradesiGFT iv,gradesiGFT iw)
                              outer = [grade (k+l) (vb*wb)
                                | k<-gv, l<-gw, let vb=gradeiGFT k iv, let wb = gradeiGFT l iw]

    v .*.  w        = foldl' (+) zero outer
                        where (iv,iw) = (invGFT $ cl v, invGFT $ cl w)
                              (gv,gw) = (gradesiGFT iv,gradesiGFT iw)
                              outer = [grade (abs (k-l)) (vb*wb)
                                | k<-gv, l<-gw, let vb=gradeiGFT k iv, let wb = gradeiGFT l iw]

    v .*  w         = foldl' (+) zero outer
                        where (iv,iw) = (invGFT $ cl v, invGFT $ cl w)
                              (gv,gw) = (gradesiGFT iv,gradesiGFT iw)
                              outer = [grade (l-k) (vb*wb)
                                | k<-gv, l<-gw,l>=k, let vb=gradeiGFT k iv, let wb = gradeiGFT l iw]

    v *.  w         = foldl' (+) zero outer
                        where (iv,iw) = (invGFT $ cl v, invGFT $ cl w)
                              (gv,gw) = (gradesiGFT iv, gradesiGFT iw)
                              outer = [grade (k-l) (vb*wb)
                                | k<-gv, l<-gw,k>=l, let vb=gradeiGFT k iv, let wb = gradeiGFT l iw]

    v *# w          = scalar   (v*w)
    v `x` w         = 0.5 !* (v*w - w*v)

    isScalar   g    = 0 == maxgrd g
    isVector   g    = (1,1) == minmaxgrd g
    isBivector g    = (2,2) == minmaxgrd g

    normmax (GA c)  = maximum $ map (abs . snd) (invGFT c)
    norm _ (GA c)   = sum $ map (abs . snd) (invGFT c)

    maxgrd          = snd . minmaxgrd

    mingrd          = fst . minmaxgrd

    minmaxgrd g     = case grades g of
                          []    -> (0,0)
                          l     -> (head l, last l)

    grades (GA c)   = gradesiGFT (invGFT c)

    scalar (GA c)   = (snd . head) (invGFT c)

    grade n (GA c)  = gradeiGFT n (invGFT c)
    gradeiGFT n ic  = let res = zipWith (\g (i,v) -> if (g==n) then (i,v) else (i,0)) gaBaseGrades ic
                      in fwdGFT res

    fwdGFT []       = GA (one 0)
    fwdGFT l@(a:as) = GA (fst $ gft $ map snd l)

    rev (GA c)      = fwdGFT reverted
                        where reverted = zipWith revert gaBaseGrades (invGFT c)
                              revert = (\g (i,v) -> if even (g `div` 2) then (i,v) else (i,(-v)))

    invol (GA c)    = fwdGFT reverted
                        where reverted = zipWith revert gaBaseGrades (invGFT c)
                              revert = (\g (i,v) -> if even g then (i,v) else (i,(-v)))

    v `e` [n]       = case lookup n gaIndexVec2BasisMap of
                            Just i          -> GA (b_ i v)
                            Nothing         -> GA (one v)
    v `e` l@(n:ns)  = case Map.lookup norml gaIndexSet2BasisMap of
                            Just (i,sgn)    -> GA (b_ i (normv sgn))
                            Nothing         -> GA (one v)
                        where (norml,sgn2)  = nubPairB $ sort l
                              normv sgn1    = if sgn1 `xor` sgn2 `xor` (odd $ numTrans l) then (-v) else v
                              a `xor` b     = a && not b || not a && b
-- helper function
gradesiGFT :: [(Int,Double)] -> [Int]
gradesiGFT ic   = nub $ map fst $ sort $ filter ((/=0).snd) res
                    where res = zipWith (\g (i,v) -> (g,v)) gaBaseGrades ic

-- base for multivector (tensor algebra in general) (inspired by David Amos)
newtype Basis = E [Int] deriving (Eq,Ord)
instance Show Basis where
    show (E []) = " "
    show (E i)  = 'e': show i
instance Read Basis where
    readsPrec _ = readBasis

readBasis :: ReadS Basis
readBasis []        = [(E [], [])] -- assume a scalar
readBasis (v:vs)
        | v == ' '  = readBasis vs
        | v == 'e'  = if (not $ null basis) then [(E basis, rest)] else []
        | otherwise = []
    where parsed        = readList vs :: [([Int], String)]
          (basis, rest) = if (not $ null parsed) then head parsed else ([], vs)

newtype NPoly r v = NP [(v,r)]
instance (Show r, Eq v, Show v) => Show (NPoly r v) where
    show (NP []) = "0"
    show (NP ts) =
        let (c:cs) = concatMap showTerm ts
        in if c == '+' then cs else c:cs
        where showTerm (m,a) =
                case show a of
                "1" -> "+" ++ show m
                "-1" -> "-" ++ show m
                cs -> showCoeff cs ++ show m
              showCoeff (c:cs) = if any (`elem` ['+','-']) cs
                                    then "+(" ++ c:cs ++ ")"
                                    else if c == '-' then c:cs else '+':c:cs

-- Clifford Algebra Transformations:
-- invGFT implements a fast inverse generalized fourier transform
--        on the Clifford Algebra. This is rather close to the definition
--        in "A generalized FFT for Clifford algebras", Paul Leopardi

-- (wanted to work with type families here, for specialization
-- between integer and floating point case!
-- My experimentation wasn't successful, type equivalence in superclass isn't
-- implemented in GHC till Milestone 7.2.x!)

class CliffordT b r | b -> r where
    invGFT :: b -> [(Int, r)]

instance CliffordT (I r) r where
    invGFT (I r) = [(0,r)]

instance (Num b,Clifford b r,CliffordT b r) => CliffordT (Complex b) r where
    invGFT (C a b) = do (id1,elem) <- zip [0..] [a,b]
                        (id2, v) <- invGFT elem
                        return (id1+2*id2, v)

instance (Num b, CliffordT b Double) => CliffordT (Split b) Double where
    invGFT (S a b) = do (idx1,elem) <- zip [0..] [a+b,a-b]
                        (idx2, v) <- invGFT elem
                        return (idx1+2*idx2, v / 2.0)

instance (Num b,Clifford b r,CliffordT b r) => CliffordT (Quaternion b) r where
    invGFT (Q a b c d) = do (id1,elem) <- zip [0..] [a,b,c,d]
                            (id2, v) <- invGFT elem
                            return (id1+4*id2, v)

instance (Num b,Clifford b Double,CliffordT b Double) => CliffordT (Matrix b) Double where
    invGFT (M a b c d) = do (idx1,elem) <- zip [0..] [a+d,d-a,b-c,b+c]
                            (idx2, v) <- invGFT elem
                            return (idx1+4*idx2, v / 2.0)

-- make sense of the index value pairs of the fourier transformed
-- Clifford Algebra via an index map (only use non-zero elements)
decode cl idMap = let nonzeros = filter ((/=) 0 . snd) $ invGFT cl
                  in case mapM (\(id,v) -> IM.lookup id idMap) nonzeros of
                      Nothing     -> Nothing
                      Just frames -> Just $ zipWith (\(b,s) (id,v) ->(E b, adjSgn s v)) frames nonzeros
                        where adjSgn :: Num a => Bool -> a -> a
                              adjSgn s v = if s then -v else v

-- We need to be able to work with frame sets, i.e vectors, bi-vectors, ---
-- and need a sensible nomenclature

-- given a permutation of [1..n] (as a list),
-- returns the number of transpositions which led to it
numTrans [] = 0
numTrans xs = numTrans' 0 [] xs where
    numTrans' i ls [r] = i + numTrans (reverse ls)
    numTrans' i ls (r1:r2:rs) =
        if r1 <= r2
        then numTrans' i (r1:ls) (r2:rs)     -- no swap needed
        else numTrans' (i+1) (r2:ls) (r1:rs) -- swap needed

-- eliminate duplicate basis elements pairwise in a sorted list
-- coevaluate signature
nubPairB :: [Int] -> ([Int],Bool)
nubPairB [] = ([],False)
nubPairB [a] = ([a],False)
nubPairB (a:b:bs)
        | a==b && a>0   = nubPairB bs
        | a==b && b<0   = let (r,s) = nubPairB bs
                          in (r,not s)
        | otherwise     = let (r,s) = nubPairB (b:bs)
                          in (a:r,s)

obviousChoices :: Int -> [[Int]]
obviousChoices 0 = []
obviousChoices 1 = [[1]]
obviousChoices 2 = [[1,2],[1,3]]
obviousChoices 3 = [[1,2,7],[1,3,6]]
obviousChoices 4 = [[1,2,7,11],[1,2,7,15],[1,3,6,10],[1,3,6,14]]
obviousChoices 5 = [[1,2,7,11,31],[1,2,7,15,27],[1,3,6,10,30],[1,3,6,14,26]]
-- should be possible to find a sequence generator for this
obviousChoices dim = [map (\i -> getIndex $ GA (e_ i 1)) [1..dim]]
    where getIndex = fst . head . (filter ((/=) 0 . snd)) . invGFT .cl

signature :: Int -> Bool
signature i = let v = GA (b_ i 1)
              in (>0) . snd . head . invGFT . cl $ v*v

-- find a vector basis in GA that has the proper signature
fixBasis :: (Int,Int) -> Maybe [Int]
fixBasis (p,q)
    | p < 0 || q < 0    = Nothing
    | otherwise = let choices = obviousChoices (p+q)
                  in msum [if pTest set then return set else Nothing | set <- choices]
                    where pTest []       = False
                          pTest [x]      = (p == 1) && signature x
                          pTest l@(x:xs) = (==p) . length $ filter id $ fmap signature l

-- use a convenient notation for the vector basis
type Key = Int
nameBasis :: (Int,Int) -> Maybe [(Key,Int)]
nameBasis (p,q) = do basis <- fixBasis (p,q)
                     let orderedB = map snd $ sort $ zip (fmap signature basis) basis
                     return $ zip names orderedB
                         where names = take (p+q) $ [-q..(-1)] ++ [1..]

-- some other utility

-- a functional extension to filterM
ffilterM            :: (Monad m) => (a -> m Bool) -> (a -> a -> a) -> [a] -> m [a]
ffilterM _ _ []     =  return []
ffilterM p f (x:xs) = do
   flg <- p x
   ys  <- ffilterM p f xs
   return (if flg then (if null ys then [x] else (x `f` (head ys)):ys) else ys)