Skip to content

Instantly share code, notes, and snippets.

@adamnemecek
Forked from fkuehnel/Clifford.hs
Created March 29, 2019 16:55
Show Gist options
  • Save adamnemecek/8d72c494ec2d02ef5a4953a56dda8d95 to your computer and use it in GitHub Desktop.
Save adamnemecek/8d72c494ec2d02ef5a4953a56dda8d95 to your computer and use it in GitHub Desktop.
Haskell module for Clifford Algebra tensor representations, and a Geometric Algebra module utilizing the tensor representation
-- 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)
-- 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment