{-# LANGUAGE
CPP, BangPatterns, TypeFamilies, DataKinds, KindSignatures, ScopedTypeVariables,
FlexibleContexts
#-}
module Math.Algebra.Polynomial.Monomial.Exterior.Indexed where
import Data.Maybe
import Data.List
import Data.Array.Unboxed
import Data.Ord
import Data.Bits
import Data.Typeable
import GHC.TypeLits
import Data.Proxy
import Data.Foldable as F
import qualified Data.Set as Set ; import Data.Set (Set)
import qualified Data.Map as Map ; import Data.Map (Map)
import Math.Algebra.Polynomial.Class
import Math.Algebra.Polynomial.Pretty
import Math.Algebra.Polynomial.Misc
import Math.Algebra.Polynomial.FreeModule ( PartialMonoid(..) )
newtype Ext (var :: Symbol) (n :: Nat) = Ext Integer deriving (Eq,Ord,Show)
data SgnExt (var :: Symbol) (n :: Nat) = SgnExt !Sign !(Ext var n) deriving (Eq,Ord,Show)
unExt :: Ext v n -> Integer
unExt (Ext k) = k
instance KnownNat n => PartialMonoid (SgnExt var n) where
pmUnit = emptySgnExt
pmAdd = mulSgnExt
instance KnownSymbol var => Pretty (Ext var n) where
pretty monom =
case [ showX i | i <- extToList monom ] of
[] -> "(1)"
xs -> intercalate "/\\" xs
where
v = extVar monom
showX !(Index i) = v ++ show i
instance KnownSymbol var => Pretty (SgnExt var n) where
pretty (SgnExt sgn ext) = case sgn of
Plus -> '+' : pretty ext
Minus -> '-' : pretty ext
extVar :: KnownSymbol var => Ext var n -> String
extVar = symbolVal . varProxy where
varProxy :: Ext var n -> Proxy var
varProxy _ = Proxy
nOfExt :: KnownNat n => Ext var n -> Int
nOfExt = fromInteger . natVal . natProxy where
natProxy :: Ext var n -> Proxy n
natProxy _ = Proxy
nOfMbExt :: KnownNat n => Maybe (Ext var n) -> Int
nOfMbExt = fromInteger . natVal . natProxy where
natProxy :: Maybe (Ext var n) -> Proxy n
natProxy _ = Proxy
nOfSgnExt :: KnownNat n => SgnExt var n -> Int
nOfSgnExt = fromInteger . natVal . natProxy where
natProxy :: SgnExt var n -> Proxy n
natProxy _ = Proxy
nOfMbSgnExt :: KnownNat n => Maybe (SgnExt var n) -> Int
nOfMbSgnExt = fromInteger . natVal . natProxy where
natProxy :: Maybe (SgnExt var n) -> Proxy n
natProxy _ = Proxy
emptyExt :: KnownNat n => Ext v n
emptyExt = Ext 0
emptySgnExt :: KnownNat n => SgnExt v n
emptySgnExt = SgnExt Plus (Ext 0)
isEmptyExt :: Ext v n -> Bool
isEmptyExt (Ext a) = (a==0)
isNormalExt :: KnownNat n => Ext v n -> Bool
isNormalExt ext@(Ext int) = shiftR int n == 0 where n = nOfExt ext
extFromList :: KnownNat n => [Index] -> Maybe (SgnExt v n)
extFromList list = result where
result
| Map.null multiset = Just emptySgnExt
| fst (Map.findMin multiset) < Index 1 = error "extFromList: index out of bounds (too small)"
| fst (Map.findMax multiset) > Index n = Nothing
| any (>1) (Map.elems multiset) = Nothing
| otherwise = Just (SgnExt sgn $ Ext int)
n = nOfMbSgnExt result
multiset = Map.fromListWith (+) [ (idx,1) | idx <- list ]
perm = sortingPermutationAsc list
sgn = if odd (numberOfInversionsMerge perm) then Minus else Plus
int = sum' [ shiftL 1 (j-1) | Index j <- Map.keys multiset ]
extToList :: Ext v n -> [Index]
extToList (Ext int) = go 0 int where
go _ 0 = []
go !i !a = if testBit a 0
then Index (i+1) : go (i+1) (shiftR a 1)
else go (i+1) (shiftR a 1)
extFromSet :: KnownNat n => Set Index -> Maybe (Ext v n)
extFromSet set = result where
n = nOfMbExt result
result = case Set.lookupMax set of
Nothing -> Just emptyExt
Just (Index k) -> if k > n
then Nothing
else if (Set.findMin set < Index 1)
then error "extFromSet: index smaller than 1"
else Just $ Ext $ sum' [ shiftL 1 (j-1) | Index j <- Set.toList set ]
extToSet :: Ext v n -> Set Index
extToSet = Set.fromList . extToList
variableExt :: KnownNat n => Index -> Ext v n
variableExt (Index idx) = Ext $ shiftL 1 (idx-1)
mulExt :: KnownNat n => Ext v n -> Ext v n -> Maybe (SgnExt v n)
mulExt x y = extFromList (extToList x ++ extToList y)
mulExtCoeff :: (Num c, KnownNat n) => Ext v n -> Ext v n -> Maybe (Ext v n, c)
mulExtCoeff x y = case mulExt x y of
Nothing -> Nothing
Just (SgnExt sgn z) -> case sgn of
Plus -> Just (z, 1)
Minus -> Just (z, -1)
mulSgnExt :: KnownNat n => SgnExt v n -> SgnExt v n -> Maybe (SgnExt v n)
mulSgnExt (SgnExt s x) (SgnExt t y) = case mulExt x y of
Nothing -> Nothing
Just (SgnExt u z) -> Just $ SgnExt (s `mappend` t `mappend` u) z
productExt :: (KnownNat n, Foldable f) => f (Ext v n) -> Maybe (SgnExt v n)
productExt t = go Plus (F.toList t) where
go !sgn list = case list of
[] -> Just (SgnExt sgn emptyExt)
[x] -> Just (SgnExt sgn x)
(x:y:rest) -> case mulExt x y of
Nothing -> Nothing
Just (SgnExt s z) -> go (sgn `mappend` s) (z:rest)
[v1,v2,v3,v4,v5,v6] = [ variableExt (Index i) | i<-[1..6] ] :: [Ext "a" 7]
Just a = productExt [v1,v2,v3]
Just b = productExt [v4,v5,v6]
prop_graded_anticomm :: KnownNat n => Ext v n -> Ext v n -> Bool
prop_graded_anticomm x y
| isNothing mb1 && isNothing mb2 = True
| isJust mb1 && isJust mb2 = u == v && maybeFlip s == t
| otherwise = False
where
mb1 = x `mulExt` y
mb2 = y `mulExt` x
Just (SgnExt s u) = mb1
Just (SgnExt t v) = mb2
d1 = degreeExt x
d2 = degreeExt y
maybeFlip = if odd (d1*d2) then oppositeSign else id
prop_graded_anticomm_sgn :: KnownNat n => SgnExt v n -> SgnExt v n -> Bool
prop_graded_anticomm_sgn x y
| isNothing mb1 && isNothing mb2 = True
| isJust mb1 && isJust mb2 = u == v && maybeFlip s == t
| otherwise = False
where
mb1 = x `mulSgnExt` y
mb2 = y `mulSgnExt` x
Just (SgnExt s u) = mb1
Just (SgnExt t v) = mb2
d1 = degreeSgnExt x
d2 = degreeSgnExt y
maybeFlip = if odd (d1*d2) then oppositeSign else id
degreeExt :: Ext v n -> Int
degreeExt (Ext k) = popCount k
degreeSgnExt :: SgnExt v n -> Int
degreeSgnExt (SgnExt sgn ext) = degreeExt ext
newtype Permutation = Permutation (UArray Int Int) deriving (Eq,Ord)
toPermutationUnsafe :: [Int] -> Permutation
toPermutationUnsafe xs = Permutation perm where
n = length xs
perm = listArray (1,n) xs
sortingPermutationAsc :: Ord a => [a] -> Permutation
sortingPermutationAsc xs = toPermutationUnsafe (map fst sorted) where
sorted = sortBy (comparing snd) $ zip [1..] xs
isEvenPermutation :: Permutation -> Bool
isEvenPermutation = even . numberOfInversionsMerge
numberOfInversionsMerge :: Permutation -> Int
numberOfInversionsMerge (Permutation arr) = fst (sortCnt n $ elems arr) where
(_,n) = bounds arr
sortCnt :: Int -> [Int] -> (Int,[Int])
sortCnt 0 _ = (0,[] )
sortCnt 1 [x] = (0,[x])
sortCnt 2 [x,y] = if x>y then (1,[y,x]) else (0,[x,y])
sortCnt n xs = mergeCnt (sortCnt k us) (sortCnt l vs) where
k = div n 2
l = n - k
(us,vs) = splitAt k xs
mergeCnt :: (Int,[Int]) -> (Int,[Int]) -> (Int,[Int])
mergeCnt (!c,us) (!d,vs) = (c+d+e,ws) where
(e,ws) = go 0 us vs
go !k xs [] = ( k*length xs , xs )
go _ [] ys = ( 0 , ys)
go !k xxs@(x:xs) yys@(y:ys) = if x < y
then let (a,zs) = go k xs yys in (a+k, x:zs)
else let (a,zs) = go (k+1) xxs ys in (a , y:zs)