combinat 0.2.8.0 → 0.2.8.1
raw patch · 14 files changed
+660/−675 lines, 14 filesdep +combinatdep +test-frameworkdep +test-framework-quickcheck2dep ~QuickCheckdep ~basedep ~containers
Dependencies added: combinat, test-framework, test-framework-quickcheck2
Dependency ranges changed: QuickCheck, base, containers
Files
- Math/Combinat/Groups/Braid.hs +141/−48
- Math/Combinat/Groups/Braid/NF.hs +18/−28
- Math/Combinat/Groups/Free.hs +1/−1
- Math/Combinat/Helper.hs +43/−20
- Math/Combinat/Numbers/Series.hs +4/−280
- Math/Combinat/Partitions/Integer.hs +124/−41
- Math/Combinat/Partitions/Skew.hs +0/−35
- Math/Combinat/Permutations.hs +147/−179
- Math/Combinat/Sign.hs +6/−1
- Math/Combinat/Tableaux/LittlewoodRichardson.hs +1/−1
- Math/Combinat/Tableaux/Skew.hs +0/−32
- Math/Combinat/TypeLevel.hs +117/−0
- combinat.cabal +17/−9
- test/TestSuite.hs +41/−0
Math/Combinat/Groups/Braid.hs view
@@ -17,7 +17,8 @@ CPP, BangPatterns, ScopedTypeVariables, ExistentialQuantification, DataKinds, KindSignatures, Rank2Types,- TypeOperators, TypeFamilies #-}+ TypeOperators, TypeFamilies,+ StandaloneDeriving #-} module Math.Combinat.Groups.Braid where @@ -44,24 +45,22 @@ import Math.Combinat.ASCII import Math.Combinat.Sign import Math.Combinat.Helper+import Math.Combinat.TypeLevel+import Math.Combinat.Numbers.Series import Math.Combinat.Permutations ( Permutation(..) ) import qualified Math.Combinat.Permutations as P -#ifdef QUICKCHECK-import Test.QuickCheck-#endif- -------------------------------------------------------------------------------- -- * Artin generators -- | A standard Artin generator of a braid: @Sigma i@ represents twisting --- the neighbour strands @i@ and @(i+1)@, such that @#i@ goes /under/ @#(i+1)@+-- the neighbour strands @i@ and @(i+1)@, such that strand @i@ goes /under/ strand @(i+1)@. -- -- Note: The strands are numbered @1..n@. data BrGen- = Sigma !Int- | SigmaInv !Int+ = Sigma !Int -- ^ @i@ goes under @(i+1)@+ | SigmaInv !Int -- ^ @i@ goes above @(i+1)@ deriving (Eq,Ord,Show) -- | The strand (more precisely, the first of the two strands) the generator twistes@@ -118,6 +117,24 @@ withSomeBraid :: SomeBraid -> (forall n. KnownNat n => Braid n -> a) -> a withSomeBraid sbraid f = case sbraid of SomeBraid braid -> f braid +mkBraid :: (forall n. KnownNat n => Braid n -> a) -> Int -> [BrGen] -> a+mkBraid f n w = y where+ sb = someBraid n (Braid w)+ y = withSomeBraid sb f++withBraid + :: Int+ -> (forall (n :: Nat). KnownNat n => Braid n)+ -> (forall (n :: Nat). KnownNat n => Braid n -> a) + -> a+withBraid n polyBraid f = + case snat of + SomeNat pxy -> f (asProxyTypeOf1 polyBraid pxy)+ where+ snat = case someNatVal (fromIntegral n :: Integer) of+ Just sn -> sn+ Nothing -> error "withBraid: input is not a natural number"+ -------------------------------------------------------------------------------- braidWord :: Braid n -> [BrGen]@@ -201,7 +218,7 @@ theGarsideBraid :: KnownNat n => Braid n theGarsideBraid = halfTwist --- | The inner automorphism defined by @X -> Delta^-1 X Delta@, +-- | The inner automorphism defined by @tau(X) = Delta^-1 X Delta@, -- where @Delta@ is the positive half-twist. -- -- This sends each generator @sigma_j@ to @sigma_(n-j)@.@@ -212,6 +229,13 @@ f (Sigma i) = Sigma (n-i) f (SigmaInv i) = SigmaInv (n-i) ++-- | The involution @tau@ on permutations (permutation braids)+--+tauPerm :: Permutation -> Permutation+tauPerm (Permutation arr) = Permutation $ listArray (1,n) [ (n+1) - arr!(n-i) | i<-[0..n-1] ] where+ (1,n) = bounds arr+ -------------------------------------------------------------------------------- -- * Group operations @@ -437,7 +461,40 @@ (sgn,k) = brGenSignIdx g s = signValue sgn +--------------------------------------------------------------------------------+-- * Growth +-- | Bronfman's recursive formula for the reciprocial of the growth function +-- of /positive/ braids. It was already known (by Deligne) that these generating functions +-- are reciprocials of polynomials; Bronfman [1] gave a recursive formula for them.+--+-- > let count n l = length $ nub $ [ braidNormalForm w | w <- allPositiveBraidWords n l ]+-- > let convertPoly (1:cs) = zip (map negate cs) [1..]+-- > pseries' (convertPoly $ bronfmanH n) == expandBronfmanH n == [ count n l | l <- [0..] ] +--+-- * [1] Aaron Bronfman: Growth functions of a class of monoids. Preprint, 2001+--+bronfmanH :: Int -> [Int]+bronfmanH n = bronfmanHsList !! n++-- | An infinite list containing the Bronfman polynomials:+--+-- > bronfmanH n = bronfmanHsList !! n+--+bronfmanHsList :: [[Int]]+bronfmanHsList = list where+ list = map go [0..]+ go 0 = [1]+ go n = sumSeries [ sgn i $ replicate (choose2 i) 0 ++ list !! (n-i) | i<-[1..n] ]+ sgn i = if odd i then id else map negate+ choose2 k = div (k*(k-1)) 2++-- | Expands the reciprocial of @H(n)@ into an infinite power series,+-- giving the growth function of the positive braids on @n@ strands.+expandBronfmanH :: Int -> [Int]+expandBronfmanH n = pseries' (convertPoly $ bronfmanH n) where+ convertPoly (1:cs) = zip (map negate cs) [1..]+ -------------------------------------------------------------------------------- -- * ASCII diagram @@ -534,29 +591,49 @@ -} --------------------------------------------------------------------------------+-- * List of all words++-- | All positive braid words of the given length+allPositiveBraidWords :: KnownNat n => Int -> [Braid n]+allPositiveBraidWords l = braids where+ n = numberOfStrands (head braids)+ braids = map Braid $ _allPositiveBraidWords n l ++-- | All braid words of the given length+allBraidWords :: KnownNat n => Int -> [Braid n]+allBraidWords l = braids where+ n = numberOfStrands (head braids)+ braids = map Braid $ _allBraidWords n l ++-- | Untyped version of 'allPositiveBraidWords'+_allPositiveBraidWords :: Int -> Int -> [[BrGen]]+_allPositiveBraidWords n = go where+ go 0 = [[]]+ go k = [ Sigma i : rest | i<-[1..n-1] , rest <- go (k-1) ]++-- | Untyped version of 'allBraidWords'+_allBraidWords :: Int -> Int -> [[BrGen]]+_allBraidWords n = go where+ go 0 = [[]]+ go k = [ gen : rest | gen <- gens , rest <- go (k-1) ]+ gens = concat [ [ Sigma i , SigmaInv i ] | i<-[1..n-1] ]++-------------------------------------------------------------------------------- -- * Random braids -- | Random braid word of the given length randomBraidWord :: (RandomGen g, KnownNat n) => Int -> g -> (Braid n, g)-randomBraidWord len gen = (braid, gen') where- braid = Braid (map sig bjs)- n = numberOfStrands braid- (gen',bjs) = mapAccumL worker gen [1..len]-- worker !g _ = (g'',(b,j)) where- (j, g' ) = randomR (1,n-1) g- (b, g'') = random g'-- sig :: (Bool,Int) -> BrGen- sig (True ,j) = Sigma j- sig (False,j) = SigmaInv j+randomBraidWord len g = (braid, g') where+ braid = Braid w+ n = numberOfStrands braid+ (w,g') = _randomBraidWord n len g -- | Random /positive/ braid word of the given length randomPositiveBraidWord :: (RandomGen g, KnownNat n) => Int -> g -> (Braid n, g)-randomPositiveBraidWord len gen = (braid, gen') where- braid = Braid (map Sigma js)- n = numberOfStrands braid- (gen',js) = mapAccumL (\(!g) _ -> swap (randomR (1,n-1) g)) gen [1..len]+randomPositiveBraidWord len g = (braid, g') where+ braid = Braid w+ n = numberOfStrands braid+ (w,g') = _randomPositiveBraidWord n len g -------------------------------------------------------------------------------- @@ -621,31 +698,47 @@ -------------------------------------------------------------------------------- -#ifdef QUICKCHECK---- | A permutation braid made convenient to use (type-level hackery)-data PermBraid = forall n. KnownNat n => PermBraid Permutation (Braid n)--mkPermBraid :: Permutation -> PermBraid-mkPermBraid perm = - case snat of - SomeNat pxy -> PermBraid perm (asProxyTypeOf1 (permutationBraid perm) pxy)- where- n = P.permutationSize perm- Just snat = someNatVal (fromIntegral n :: Integer)--prop_permBraid_perm :: PermBraid -> Bool-prop_permBraid_perm (PermBraid perm braid) = (braidPermutation braid == perm)--prop_permBraid_valid :: PermBraid -> Bool-prop_permBraid_valid (PermBraid perm braid) = isPermutationBraid braid+-- | This version of 'randomBraidWord' may be convenient to avoid the type level stuff+withRandomBraidWord + :: RandomGen g + => (forall n. KnownNat n => Braid n -> a) + -> Int -- ^ number of strands+ -> Int -- ^ length of the random word+ -> g -> (a, g)+withRandomBraidWord f n len = runRand $ do+ withSelectedM f (rand $ randomBraidWord len) n -prop_braidPerm_comp :: KnownNat n => Braid n -> Braid n -> Bool-prop_braidPerm_comp b1 b2 = (p == q) where- p = braidPermutation (compose b1 b2) - q = braidPermutation b1 `P.multiply` braidPermutation b2+-- | This version of 'randomPositiveBraidWord' may be convenient to avoid the type level stuff+withRandomPositiveBraidWord + :: RandomGen g + => (forall n. KnownNat n => Braid n -> a) + -> Int -- ^ number of strands+ -> Int -- ^ length of the random word+ -> g -> (a, g)+withRandomPositiveBraidWord f n len = runRand $ do+ withSelectedM f (rand $ randomPositiveBraidWord len) n +-- | Untyped version of 'randomBraidWord'+_randomBraidWord + :: (RandomGen g) + => Int -- ^ number of strands+ -> Int -- ^ length of the random word+ -> g -> ([BrGen], g)+_randomBraidWord n len = runRand $ replicateM len $ do+ k <- randChoose (1,n-1)+ s <- randRoll+ return $ case s of+ Plus -> Sigma k+ Minus -> SigmaInv k -#endif+-- | Untyped version of 'randomPositiveBraidWord'+_randomPositiveBraidWord + :: (RandomGen g) + => Int -- ^ number of strands+ -> Int -- ^ length of the random word+ -> g -> ([BrGen], g)+_randomPositiveBraidWord n len = runRand $ replicateM len $ do+ liftM Sigma $ randChoose (1,n-1) --------------------------------------------------------------------------------+
Math/Combinat/Groups/Braid/NF.hs view
@@ -18,12 +18,17 @@ DataKinds, KindSignatures, Rank2Types #-} module Math.Combinat.Groups.Braid.NF - ( BraidNF (..)+ ( -- * Normal form+ BraidNF (..) , nfReprWord , braidNormalForm , braidNormalForm'-#ifdef QUICKCHECK-#endif+ , braidNormalFormNaive'+ -- * Starting and finishing sets+ , permWordStartingSet+ , permWordFinishingSet + , permutationStartingSet+ , permutationFinishingSet ) where @@ -53,10 +58,6 @@ import Math.Combinat.Groups.Braid -#ifdef QUICKCHECK-import Test.QuickCheck-#endif- -------------------------------------------------------------------------------- -- | A unique normal form for braids, called the /left-greedy normal form/.@@ -95,11 +96,19 @@ braidNormalForm' braid@(Braid gens) = BraidNF (dexp+pexp) perms where n = numberOfStrands braid invless = replaceInverses n gens- -- invless = replaceInversesNaive gens (dexp,posxword) = moveDeltasLeft n invless factors = leftGreedyFactors n $ expandPosXWord n posxword (pexp,perms) = normalizePermFactors n $ map (_braidPermutation n) factors +-- | This one uses the naive inverse replacement method. Probably somewhat slower than 'braidNormalForm''.+braidNormalFormNaive' :: KnownNat n => Braid n -> BraidNF n+braidNormalFormNaive' braid@(Braid gens) = BraidNF (dexp+pexp) perms where+ n = numberOfStrands braid+ invless = replaceInversesNaive gens+ (dexp,posxword) = moveDeltasLeft n invless+ factors = leftGreedyFactors n $ expandPosXWord n posxword+ (pexp,perms) = normalizePermFactors n $ map (_braidPermutation n) factors+ -------------------------------------------------------------------------------- -- | Replaces groups of @sigma_i^-1@ generators by @(Delta^-1 * P)@, @@ -421,7 +430,7 @@ | otherwise = (e , [oddTau (e+ep) p] ) oddTau :: Int -> Permutation -> Permutation- oddTau !e p = if even e then p else permTau p+ oddTau !e p = if even e then p else tauPerm p {- checkDelta :: Int -> Permutation -> [Permutation] -> (Int,[Permutation])@@ -431,11 +440,6 @@ | otherwise = let (e',rs) = worker e rest in (e', oddTau e p : rs) -} --- | The involution tau on permutation-permTau :: Permutation -> Permutation-permTau (Permutation arr) = Permutation $ listArray (1,n) [ (n+1) - arr!(n-i) | i<-[0..n-1] ] where- (1,n) = bounds arr- -------------------------------------------------------------------------------- -- | Given a /positive/ word, we apply left-greedy factorization of@@ -528,17 +532,3 @@ -------------------------------------------------------------------------------- -#ifdef QUICKCHECK--prop_braidnf_reduce :: KnownNat n => Braid n -> Bool-prop_braidnf_reduce braid = (braidNormalForm' braid == braidNormalForm braid)--prop_braidnf_reprs :: KnownNat n => Braid n -> Bool-prop_braidnf_reprs braid = (nf == nf') where- nf = braidNormalForm braid - nf' = braidNormalForm braid'- braid' = nfReprWord nf--#endif----------------------------------------------------------------------------------
Math/Combinat/Groups/Free.hs view
@@ -493,7 +493,7 @@ {- --- some basic testing. TODO: QuickCheck tests+-- some basic testing. TODO: real tests import Math.Combinat.Helper import Math.Combinat.Groups.Free
Math/Combinat/Helper.hs view
@@ -1,12 +1,14 @@ -- | Miscellaneous helper functions -{-# LANGUAGE BangPatterns, PolyKinds #-}+{-# LANGUAGE BangPatterns, PolyKinds, GeneralizedNewtypeDeriving #-} module Math.Combinat.Helper where -------------------------------------------------------------------------------- import Control.Monad+import Control.Applicative ( Applicative(..) ) -- required before AMP (before GHC 7.10)+import Data.Functor.Identity import Data.List import Data.Ord@@ -17,6 +19,9 @@ import Debug.Trace +import System.Random+import Control.Monad.Trans.State+ -------------------------------------------------------------------------------- -- * debugging @@ -24,24 +29,6 @@ debug x y = trace ("-- " ++ show x ++ "\n") y ----------------------------------------------------------------------------------- * proxy--proxyUndef :: Proxy a -> a-proxyUndef _ = error "proxyUndef"--proxyOf :: a -> Proxy a-proxyOf _ = Proxy--proxyOf1 :: f a -> Proxy a-proxyOf1 _ = Proxy--proxyOf2 :: g (f a) -> Proxy a-proxyOf2 _ = Proxy--asProxyTypeOf1 :: f a -> Proxy a -> f a -asProxyTypeOf1 y _ = y---------------------------------------------------------------------------------- -- * pairs swap :: (a,b) -> (b,a)@@ -254,4 +241,40 @@ -} --------------------------------------------------------------------------------- +-- * random++-- | A simple random monad to make life suck less+type Rand g = RandT g Identity++runRand :: Rand g a -> g -> (a,g)+runRand action g = runIdentity (runRandT action g)++flipRunRand :: Rand s a -> s -> (s,a)+flipRunRand action g = runIdentity (flipRunRandT action g)+++-- | The Rand monad transformer+newtype RandT g m a = RandT (StateT g m a) deriving (Functor,Applicative,Monad)++runRandT :: RandT g m a -> g -> m (a,g)+runRandT (RandT stuff) = runStateT stuff++-- | This may be occasionally useful+flipRunRandT :: Monad m => RandT s m a -> s -> m (s,a)+flipRunRandT action ini = liftM swap $ runRandT action ini+++-- | Puts a standard-conforming random function into the monad+rand :: (g -> (a,g)) -> Rand g a+rand user = RandT (state user)++randRoll :: (RandomGen g, Random a) => Rand g a+randRoll = rand random++randChoose :: (RandomGen g, Random a) => (a,a) -> Rand g a+randChoose uv = rand (randomR uv)++randProxy1 :: Rand g (f n) -> Proxy n -> Rand g (f n)+randProxy1 action _ = action++--------------------------------------------------------------------------------
Math/Combinat/Numbers/Series.hs view
@@ -20,11 +20,6 @@ import Math.Combinat.Partitions.Integer import Math.Combinat.Helper -#ifdef QUICKCHECK-import System.Random-import Test.QuickCheck-#endif- -------------------------------------------------------------------------------- -- * Trivial series @@ -55,6 +50,10 @@ addSeries :: Num a => [a] -> [a] -> [a] addSeries xs ys = longZipWith 0 0 (+) xs ys +sumSeries :: Num a => [[a]] -> [a]+sumSeries [] = [0]+sumSeries xs = foldl1' addSeries xs+ subSeries :: Num a => [a] -> [a] -> [a] subSeries xs ys = longZipWith 0 0 (-) xs ys @@ -301,25 +300,6 @@ -- Reciprocals of polynomials, without coefficients -#ifdef QUICKCHECK--- | Expansion of @1 / (1-x^k)@. Included for completeness only; --- it equals to @coinSeries [k]@, and for example--- for @k=4@ it is simply--- --- > [1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0...]----pseries1 :: Int -> [Integer]-pseries1 k1 = convolveWithPSeries1 k1 unitSeries ---- | The expansion of @1 / (1-x^k_1-x^k_2)@-pseries2 :: Int -> Int -> [Integer]-pseries2 k1 k2 = convolveWithPSeries2 k1 k2 unitSeries ---- | The expansion of @1 / (1-x^k_1-x^k_2-x^k_3)@-pseries3 :: Int -> Int -> Int -> [Integer]-pseries3 k1 k2 k3 = convolveWithPSeries3 k1 k2 k3 unitSeries-#endif- -- | The power series expansion of -- -- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)@@ -327,33 +307,6 @@ pseries :: [Int] -> [Integer] pseries ks = convolveWithPSeries ks unitSeries -#ifdef QUICKCHECK--- | Convolve with (the expansion of) @1 / (1-x^k1)@-convolveWithPSeries1 :: Int -> [Integer] -> [Integer]-convolveWithPSeries1 k1 series1 = xs where- series = series1 ++ repeat 0 - xs = zipWith (+) series ( replicate k1 0 ++ xs )---- | Convolve with (the expansion of) @1 / (1-x^k1-x^k2)@-convolveWithPSeries2 :: Int -> Int -> [Integer] -> [Integer]-convolveWithPSeries2 k1 k2 series1 = xs where- series = series1 ++ repeat 0 - xs = zipWith3 (\x y z -> x + y + z)- series- ( replicate k1 0 ++ xs )- ( replicate k2 0 ++ xs )- --- | Convolve with (the expansion of) @1 / (1-x^k_1-x^k_2-x^k_3)@-convolveWithPSeries3 :: Int -> Int -> Int -> [Integer] -> [Integer]-convolveWithPSeries3 k1 k2 k3 series1 = xs where- series = series1 ++ repeat 0 - xs = zipWith4 (\x y z w -> x + y + z + w)- series- ( replicate k1 0 ++ xs )- ( replicate k2 0 ++ xs )- ( replicate k3 0 ++ xs )-#endif- -- | Convolve with (the expansion of) -- -- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)@@ -369,24 +322,6 @@ -------------------------------------------------------------------------------- -- Reciprocals of polynomials, with coefficients -#ifdef QUICKCHECK--- | @1 / (1 - a*x^k)@. --- For example, for @a=3@ and @k=2@ it is just--- --- > [1,0,3,0,9,0,27,0,81,0,243,0,729,0,2187,0,6561,0,19683,0...]----pseries1' :: Num a => (a,Int) -> [a]-pseries1' ak1 = convolveWithPSeries1' ak1 unitSeries---- | @1 / (1 - a_1*x^k_1 - a_2*x^k_2)@-pseries2' :: Num a => (a,Int) -> (a,Int) -> [a]-pseries2' ak1 ak2 = convolveWithPSeries2' ak1 ak2 unitSeries---- | @1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3)@-pseries3' :: Num a => (a,Int) -> (a,Int) -> (a,Int) -> [a]-pseries3' ak1 ak2 ak3 = convolveWithPSeries3' ak1 ak2 ak3 unitSeries-#endif- -- | The expansion of -- -- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)@@ -394,35 +329,6 @@ pseries' :: Num a => [(a,Int)] -> [a] pseries' aks = convolveWithPSeries' aks unitSeries -#ifdef QUICKCHECK--- | Convolve with @1 / (1 - a*x^k)@. -convolveWithPSeries1' :: Num a => (a,Int) -> [a] -> [a]-convolveWithPSeries1' (a1,k1) series1 = xs where- series = series1 ++ repeat 0 - xs = zipWith (+)- series- ( replicate k1 0 ++ map (*a1) xs )---- | Convolve with @1 / (1 - a_1*x^k_1 - a_2*x^k_2)@-convolveWithPSeries2' :: Num a => (a,Int) -> (a,Int) -> [a] -> [a]-convolveWithPSeries2' (a1,k1) (a2,k2) series1 = xs where- series = series1 ++ repeat 0 - xs = zipWith3 (\x y z -> x + y + z)- series- ( replicate k1 0 ++ map (*a1) xs )- ( replicate k2 0 ++ map (*a2) xs )- --- | Convolve with @1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3)@-convolveWithPSeries3' :: Num a => (a,Int) -> (a,Int) -> (a,Int) -> [a] -> [a]-convolveWithPSeries3' (a1,k1) (a2,k2) (a3,k3) series1 = xs where- series = series1 ++ repeat 0 - xs = zipWith4 (\x y z w -> x + y + z + w)- series- ( replicate k1 0 ++ map (*a1) xs )- ( replicate k2 0 ++ map (*a2) xs )- ( replicate k3 0 ++ map (*a3) xs )-#endif- -- | Convolve with (the expansion of) -- -- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)@@ -467,186 +373,4 @@ -------------------------------------------------------------------------------- -#ifdef QUICKCHECK---- * Tests--{--swap :: (a,b) -> (b,a)-swap (x,y) = (y,x)--}---- compare the first 1000 elements of the infinite lists-(=!=) :: (Eq a, Num a) => [a] -> [a] -> Bool-(=!=) xs1 ys1 = (take m xs == take m ys) where - m = 1000- xs = xs1 ++ repeat 0- ys = ys1 ++ repeat 0-infix 4 =!=--newtype Nat = Nat { fromNat :: Int } deriving (Eq,Ord,Show,Num,Random)-newtype Ser = Ser { fromSer :: [Integer] } deriving (Eq,Ord,Show)-newtype Exp = Exp { fromExp :: Int } deriving (Eq,Ord,Show,Num,Random)-newtype Exps = Exps { fromExps :: [Int] } deriving (Eq,Ord,Show)-newtype CoeffExp = CoeffExp { fromCoeffExp :: (Integer,Int) } deriving (Eq,Ord,Show)-newtype CoeffExps = CoeffExps { fromCoeffExps :: [(Integer,Int)] } deriving (Eq,Ord,Show)--minSerSize = 0 :: Int-maxSerSize = 1000 :: Int--minSerValue = -10000 :: Integer-maxSerValue = 10000 :: Integer--rndList :: (RandomGen g, Random a) => Int -> (a, a) -> g -> ([a], g)-rndList n minmax g = swap $ mapAccumL f g [1..n] where- f g _ = swap $ randomR minmax g --instance Arbitrary Nat where- arbitrary = choose (Nat 0 , Nat 750)--instance Arbitrary Exp where- arbitrary = choose (Exp 1 , Exp 32)--instance Arbitrary CoeffExp where- arbitrary = do- coeff <- choose (minSerValue, maxSerValue) :: Gen Integer- exp <- arbitrary :: Gen Exp- return $ CoeffExp (coeff,fromExp exp)- -instance Random Ser where- random g = (Ser list, g2) where- (size,g1) = randomR (minSerSize,maxSerSize) g- (list,g2) = rndList size (minSerValue,maxSerValue) g1- randomR _ = random--instance Random Exps where- random g = (Exps list, g2) where- (size,g1) = randomR (0,10) g- (list,g2) = rndList size (1,32) g1- randomR _ = random--instance Random CoeffExps where- random g = (CoeffExps (zip list2 list1), g3) where- (size,g1) = randomR (0,10) g- (list1,g2) = rndList size (1,32) g1- (list2,g3) = rndList size (minSerValue,maxSerValue) g2- randomR _ = random- -instance Arbitrary Ser where- arbitrary = choose undefined--instance Arbitrary Exps where- arbitrary = choose undefined--instance Arbitrary CoeffExps where- arbitrary = choose undefined- --- TODO: quickcheck test properties--checkAll :: IO ()-checkAll = do- let f :: Testable a => a -> IO ()- f = quickCheck-{- - -- these are very slow, because random is slow- putStrLn "leftIdentity" ; f prop_leftIdentity- putStrLn "rightIdentity" ; f prop_rightIdentity- putStrLn "commutativity" ; f prop_commutativity- putStrLn "associativity" ; f prop_associativity--}- putStrLn "convPSeries1 vs generic" ; f prop_conv1_vs_gen- putStrLn "convPSeries2 vs generic" ; f prop_conv2_vs_gen- putStrLn "convPSeries3 vs generic" ; f prop_conv3_vs_gen- putStrLn "convPSeries1' vs generic" ; f prop_conv1_vs_gen'- putStrLn "convPSeries2' vs generic" ; f prop_conv2_vs_gen'- putStrLn "convPSeries3' vs generic" ; f prop_conv3_vs_gen'- putStrLn "convolve_pseries" ; f prop_convolve_pseries - putStrLn "convolve_pseries'" ; f prop_convolve_pseries' - putStrLn "coinSeries vs pseries" ; f prop_coin_vs_pseries- putStrLn "coinSeries vs pseries'" ; f prop_coin_vs_pseries'- -prop_leftIdentity ser = ( xs =!= unitSeries `convolve` xs ) where - xs = fromSer ser --prop_rightIdentity ser = ( unitSeries `convolve` xs =!= xs ) where - xs = fromSer ser --prop_commutativity ser1 ser2 = ( xs `convolve` ys =!= ys `convolve` xs ) where - xs = fromSer ser1- ys = fromSer ser2--prop_associativity ser1 ser2 ser3 = ( one =!= two ) where- one = (xs `convolve` ys) `convolve` zs- two = xs `convolve` (ys `convolve` zs)- xs = fromSer ser1- ys = fromSer ser2- zs = fromSer ser3- -prop_conv1_vs_gen exp1 ser = ( one =!= two ) where- one = convolveWithPSeries1 k1 xs - two = convolveWithPSeries [k1] xs- k1 = fromExp exp1- xs = fromSer ser --prop_conv2_vs_gen exp1 exp2 ser = (one =!= two) where- one = convolveWithPSeries2 k1 k2 xs - two = convolveWithPSeries [k2,k1] xs- k1 = fromExp exp1- k2 = fromExp exp2- xs = fromSer ser --prop_conv3_vs_gen exp1 exp2 exp3 ser = (one =!= two) where- one = convolveWithPSeries3 k1 k2 k3 xs - two = convolveWithPSeries [k2,k3,k1] xs- k1 = fromExp exp1- k2 = fromExp exp2- k3 = fromExp exp3- xs = fromSer ser --prop_conv1_vs_gen' exp1 ser = ( one =!= two ) where- one = convolveWithPSeries1' ak1 xs - two = convolveWithPSeries' [ak1] xs- ak1 = fromCoeffExp exp1- xs = fromSer ser --prop_conv2_vs_gen' exp1 exp2 ser = (one =!= two) where- one = convolveWithPSeries2' ak1 ak2 xs - two = convolveWithPSeries' [ak2,ak1] xs- ak1 = fromCoeffExp exp1- ak2 = fromCoeffExp exp2- xs = fromSer ser --prop_conv3_vs_gen' exp1 exp2 exp3 ser = (one =!= two) where- one = convolveWithPSeries3' ak1 ak2 ak3 xs - two = convolveWithPSeries' [ak2,ak3,ak1] xs- ak1 = fromCoeffExp exp1- ak2 = fromCoeffExp exp2- ak3 = fromCoeffExp exp3- xs = fromSer ser --prop_convolve_pseries exps1 ser = (one =!= two) where- one = convolveWithPSeries ks1 xs - two = xs `convolve` pseries ks1 - ks1 = fromExps exps1- xs = fromSer ser --prop_convolve_pseries' cexps1 ser = (one =!= two) where- one = convolveWithPSeries' aks1 xs - two = xs `convolve` pseries' aks1 - aks1 = fromCoeffExps cexps1- xs = fromSer ser --prop_coin_vs_pseries exps1 = (one =!= two) where- one = coinSeries ks1 - two = convolveMany (map pseries1 ks1)- ks1 = fromExps exps1--prop_coin_vs_pseries' cexps1 = (one =!= two) where- one = coinSeries' aks1 - two = convolveMany (map pseries1' aks1)- aks1 = fromCoeffExps cexps1- -#endif ----------------------------------------------------------------------------------
Math/Combinat/Partitions/Integer.hs view
@@ -17,12 +17,13 @@ -- <<svg/ferrers.svg>> -- -{-# LANGUAGE CPP, BangPatterns #-}+{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables #-} module Math.Combinat.Partitions.Integer where -------------------------------------------------------------------------------- import Data.List+import Control.Monad ( liftM , replicateM ) -- import Data.Map (Map) -- import qualified Data.Map as Map@@ -32,6 +33,9 @@ import Math.Combinat.Numbers (factorial,binomial,multinomial) import Math.Combinat.Helper +import Data.Array+import System.Random+ -------------------------------------------------------------------------------- -- * Type and basic stuff @@ -208,9 +212,69 @@ go _ 0 = [[]] go !h !n = [ a:as | a<-[1..min n h], as <- go a (n-a) ] +-- | Number of partitions of @n@ countPartitions :: Int -> Integer-countPartitions d = countPartitions' (d,d) d+countPartitions n = partitionCountList !! n +-- | This uses 'countPartitions'', and thus is slow+countPartitionsNaive :: Int -> Integer+countPartitionsNaive d = countPartitions' (d,d) d++--------------------------------------------------------------------------------++-- | Infinite list of number of partitions of @0,1,2,...@+--+-- This uses the infinite product formula the generating function of partitions, recursively+-- expanding it; it is quite fast.+--+-- > partitionCountList == map countPartitions [0..]+--+partitionCountList :: [Integer]+partitionCountList = final where++ final = go 1 (1:repeat 0) ++ go !k (x:xs) = x : go (k+1) ys where+ ys = zipWith (+) xs (take k final ++ ys)+ -- explanation:+ -- xs == drop k $ f (k-1)+ -- ys == drop k $ f (k ) ++{-++Full explanation of 'partitionCountList':+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~++let f k = productPSeries $ map (:[]) [1..k]++f 0 = [1,0,0,0,0,0,0,0...]+f 1 = [1,1,1,1,1,1,1,1...]+f 2 = [1,1,2,2,3,3,4,4...]+f 3 = [1,1,2,3,4,5,7,8...]++observe: ++* take (k+1) (f k) == take (k+1) partitionCountList+* f (k+1) == zipWith (+) (f k) (replicate (k+1) 0 ++ f (k+1))++now apply (drop (k+1)) to the second one : ++* drop (k+1) (f (k+1)) == zipWith (+) (drop (k+1) $ f k) (f (k+1))+* f (k+1) = take (k+1) final ++ drop (k+1) (f (k+1))++-}++--------------------------------------------------------------------------------++-- | Naive infinite list of number of partitions of @0,1,2,...@+--+-- > partitionCountListNaive == map countPartitionsNaive [0..]+--+-- This is much slower than the power series expansion above.+--+partitionCountListNaive :: [Integer]+partitionCountListNaive = map countPartitionsNaive [0..]+ -- | All integer partitions up to a given degree (that is, all integer partitions whose sum is less or equal to @d@) allPartitions :: Int -> [Partition] allPartitions d = concat [ partitions i | i <- [0..d] ]@@ -266,7 +330,65 @@ countPartitions' (h,w) d = sum [ countPartitions' (i,w-1) (d-i) | i <- [1..min d h] ] + ---------------------------------------------------------------------------------+-- * Random partitions++-- | Uniformly random partition of the given weight. +--+-- NOTE: This algorithm is effective for small @n@-s (say @n@ up to a few hundred \/ one thousand it should work nicely),+-- and the first time it is executed may be slower (as it needs to build the table 'partitionCountList' first)+--+-- Algorithm of Nijenhuis and Wilf (1975); see+--+-- * Knuth Vol 4A, pre-fascicle 3B, exercise 47;+--+-- * Nijenhuis and Wilf: Combinatorial Algorithms for Computers and Calculators, chapter 10+--+randomPartition :: RandomGen g => Int -> g -> (Partition, g)+randomPartition n g = (p, g') where+ ([p], g') = randomPartitions 1 n g++-- | Generates several uniformly random partitions of @n@ at the same time.+-- Should be a little bit faster then generating them individually.+--+randomPartitions + :: forall g. RandomGen g + => Int -- ^ number of partitions to generate+ -> Int -- ^ the weight of the partitions+ -> g -> ([Partition], g)+randomPartitions howmany n = runRand $ replicateM howmany (worker n []) where++ table = listArray (0,n) $ take (n+1) partitionCountList :: Array Int Integer+ cnt k = table ! k+ + finish :: [(Int,Int)] -> Partition+ finish = mkPartition . concatMap f where f (j,d) = replicate j d++ fi :: Int -> Integer + fi = fromIntegral++ find_jd :: Int -> Integer -> (Int,Int)+ find_jd m capm = go 0 [ (j,d) | j<-[1..n], d<-[1..div m j] ] where+ go :: Integer -> [(Int,Int)] -> (Int,Int)+ go !s [] = (1,1) -- ??+ go !s [jd] = jd -- ??+ go !s (jd@(j,d):rest) = + if s' > capm + then jd + else go s' rest+ where+ s' = s + fi d * cnt (m - j*d)++ worker :: Int -> [(Int,Int)] -> Rand g Partition+ worker 0 acc = return $ finish acc+ worker !m acc = do+ capm <- randChoose (0, (fi m) * cnt m - 1)+ let jd@(!j,!d) = find_jd m capm+ worker (m - j*d) (jd:acc)+++--------------------------------------------------------------------------------- -- * Dominance order -- | @q \`dominates\` p@ returns @True@ if @q >= p@ in the dominance order of partitions@@ -584,42 +706,3 @@ -------------------------------------------------------------------------------- -#ifdef QUICKCHECK---- * Tests---- we need some custom types for quickcheck generating not too large partitions...--newtype PartitionWeight = PartitionWeight Int-data PartitionWeightPair = PartitionWeightPair Int Int -data PartitionPair = PartitionPair Partition Int --prop_partitions_in_bigbox :: PartitionWeight -> Bool-prop_partitions_in_bigbox (PartitionWeight n) = (partitions n == partitions' (n,n) n)--prop_kparts :: PartitionWeightPair -> Bool-prop_kparts (PartitionWeightPair n k) = (partitionsWithKParts k n == [ mu | mu <- partitions n, numberOfParts mu == k ])--prop_odd_partitions :: PartitionWeight -> Bool-prop_odd_partitions (PartitionWeight n) = - (partitionsWithOddParts n == [ mu | mu <- partitions n, and (map odd (fromPartition mu)) ])--prop_distinct_partitions :: PartitionWeight -> Bool-prop_distinct_partitions (PartitionWeight n) = - (partitionsWithDistinctParts n == [ mu | mu <- partitions n, let xs = fromPartition mu, xs == nub xs ])--prop_subparts :: PartitionPair -> Bool-prop_subparts (PartitionPair lam d) = (subPartitions d lam) == sort [ p | p <- partitions d, isSubPartitionOf p lam ]--prop_dual_dual :: Partition -> Bool-prop_dual_dual lam = (lam == dualPartition (dualPartition lam))--prop_dominated_list :: Partition -> Bool-prop_dominated_list lam = (dominatedPartitions lam == [ mu | mu <- partitions (weight lam), lam `dominates` mu ])--prop_dominating_list :: Partition -> Bool-prop_dominating_list mu = (dominatingPartitions mu == [ lam | lam <- partitions (weight mu ), lam `dominates` mu ])--#endif----------------------------------------------------------------------------------
Math/Combinat/Partitions/Skew.hs view
@@ -133,38 +133,3 @@ -------------------------------------------------------------------------------- -#ifdef QUICKCHECK--prop_dual_dual :: SkewPartition -> Bool-prop_dual_dual sp = (dualSkewPartition (dualSkewPartition sp) == sp)--prop_dual_from :: SkewPartition -> Bool-prop_dual_from sp = (p==p' && q==q') where- (p,q) = fromSkewPartition sp- sp' = dualSkewPartition sp- (p',q') = fromSkewPartition sp'--prop_from_to :: SkewPartition -> Bool-prop_from_to sp = (mkSkewPartition (fromSkewPartition sp) == sp)--prop_to_from :: (Partition,Partition) -> Bool-prop_to_from (p,q) = - case mb of- Nothing -> True- Just sp -> fromSkewPartition sp == (p,q)- where- mb = safeSkewPartition (p,q)--prop_from_to_from :: SkewPartition -> Bool-prop_from_to_from sp = (pq == pq') where- pq = fromSkewPartition sp- sp' = mkSkewPartition pq- pq' = fromSkewPartition sp'--prop_weight :: SkewPartition -> Bool-prop_weight sp = (skewPartitionWeight sp == weight p - weight q) where- (p,q) = fromSkewPartition sp--#endif----------------------------------------------------------------------------------
Math/Combinat/Permutations.hs view
@@ -8,7 +8,7 @@ -- are represented internally. Also now they act on the /right/ by default! -- -{-# LANGUAGE CPP, ScopedTypeVariables, GeneralizedNewtypeDeriving, FlexibleContexts #-}+{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, GeneralizedNewtypeDeriving, FlexibleContexts #-} module Math.Combinat.Permutations ( -- * The Permutation type Permutation (..)@@ -45,6 +45,13 @@ , cycleLeft , cycleRight , reversePermutation+ -- * Inversions+ , inversions+ , numberOfInversions+ , numberOfInversionsNaive+ , numberOfInversionsMerge+ , bubbleSort2+ , bubbleSort -- * Permutation groups , identity , inverse@@ -60,7 +67,8 @@ , asciiPermutation , asciiDisjointCycles , twoLineNotation - , twoLineNotation'+ , inverseTwoLineNotation+ , genericTwoLineNotation -- * List of permutations , permutations , _permutations@@ -78,11 +86,6 @@ , permuteMultiset , countPermuteMultiset , fasc2B_algorithm_L--#ifdef QUICKCHECK- -- * QuickCheck - , checkAll-#endif ) where @@ -91,7 +94,8 @@ import Control.Monad import Control.Monad.ST -import Data.List hiding (permutations)+import Data.List hiding ( permutations )+import Data.Ord ( comparing ) import Data.Array (Array) import Data.Array.ST@@ -104,18 +108,15 @@ import Math.Combinat.Classes import Math.Combinat.Helper import Math.Combinat.Sign-import Math.Combinat.Numbers (factorial,binomial)+import Math.Combinat.Numbers ( factorial , binomial ) import System.Random -#ifdef QUICKCHECK-import Test.QuickCheck-#endif- -------------------------------------------------------------------------------- -- * Types --- | A permutation. Internally it is an (unboxed) array of the integers @[1..n]@. +-- | A permutation. Internally it is an (unboxed) array of the integers @[1..n]@, with +-- indexing range also being @(1,n)@. -- -- If this array of integers is @[p1,p2,...,pn]@, then in two-line -- notations, that represents the permutation@@ -186,7 +187,7 @@ -- the zero index is an unidiomatic hack ar = (accumArray (+) 0 (0,n) $ map f xs) :: UArray Int Int f :: Int -> (Int,Int)- f j = if j<1 || j>n then (0,1) else (j,1)+ f !j = if j<1 || j>n then (0,1) else (j,1) -- | Checks whether the input is a permutation of the numbers @[1..n]@. maybePermutation :: [Int] -> Maybe Permutation@@ -233,16 +234,28 @@ asciiDisjointCycles :: DisjointCycles -> ASCII asciiDisjointCycles (DisjointCycles cycles) = final where final = hCatWith VTop (HSepSpaces 1) boxes - boxes = [ twoLineNotation' (f cyc) | cyc <- cycles ]+ boxes = [ genericTwoLineNotation (f cyc) | cyc <- cycles ] f cyc = pairs (cyc ++ [head cyc]) -- | The standard two-line notation, moving the element indexed by the top row into -- the place indexed by the corresponding element in the bottom row. twoLineNotation :: Permutation -> ASCII-twoLineNotation (Permutation arr) = twoLineNotation' $ zip [1..] (elems arr)+twoLineNotation (Permutation arr) = genericTwoLineNotation $ zip [1..] (elems arr) -twoLineNotation' :: [(Int,Int)] -> ASCII-twoLineNotation' xys = asciiFromLines [ topLine, botLine ] where+-- | The inverse two-line notation, where the it\'s the bottom line +-- which is in standard order. The columns of this are a permutation+-- of the columns 'twoLineNotation'.+--+-- Remark: the top row of @inverseTwoLineNotation perm@ is the same +-- as the bottom row of @twoLineNotation (inverse perm)@.+--+inverseTwoLineNotation :: Permutation -> ASCII+inverseTwoLineNotation (Permutation arr) =+ genericTwoLineNotation $ sortBy (comparing snd) $ zip [1..] (elems arr) ++-- | Two-line notation for any set of numbers+genericTwoLineNotation :: [(Int,Int)] -> ASCII+genericTwoLineNotation xys = asciiFromLines [ topLine, botLine ] where topLine = "( " ++ intercalate " " us ++ " )" botLine = "( " ++ intercalate " " vs ++ " )" pairs = [ (show x, show y) | (x,y) <- xys ]@@ -389,6 +402,119 @@ DisjointCycles cycles = permutationToDisjointCycles perm --------------------------------------------------------------------------------+-- * Inversions++-- | An /inversion/ of a permutation @sigma@ is a pair @(i,j)@ such that+-- @i<j@ and @sigma(i) > sigma(j)@.+--+-- This functions returns the inversion of a permutation.+--+inversions :: Permutation -> [(Int,Int)]+inversions (Permutation arr) = list where+ (_,n) = bounds arr+ list = [ (i,j) | i<-[1..n-1], j<-[i+1..n], arr!i > arr!j ]++-- | Returns the number of inversions:+--+-- > numberOfInversions perm = length (inversions perm)+--+-- Synonym for 'numberOfInversionsMerge'+--+numberOfInversions :: Permutation -> Int+numberOfInversions = numberOfInversionsMerge++-- | Returns the number of inversions, using the merge-sort algorithm.+-- This should be @O(n*log(n))@+--+numberOfInversionsMerge :: Permutation -> Int+numberOfInversionsMerge (Permutation arr) = fst (sortCnt n $ elems arr) where+ (_,n) = bounds arr+ + -- | First argument is length of the list.+ -- Returns also the inversion count.+ 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)++-- | Returns the number of inversions, using the definition, thus it's @O(n^2)@.+--+numberOfInversionsNaive :: Permutation -> Int+numberOfInversionsNaive (Permutation arr) = length list where+ (_,n) = bounds arr+ list = [ (0::Int) | i<-[1..n-1], j<-[i+1..n], arr!i > arr!j ]++-- | Bubble sorts breaks a permutation into the product of adjacent transpositions:+--+-- > multiplyMany' n (map (transposition n) $ bubbleSort2 perm) == perm+--+-- Note that while this is not unique, the number of transpositions +-- equals the number of inversions.+--+bubbleSort2 :: Permutation -> [(Int,Int)]+bubbleSort2 = map f . bubbleSort where f i = (i,i+1)++-- | Another version of bubble sort. An entry @i@ in the return sequence means+-- the transposition @(i,i+1)@:+--+-- > multiplyMany' n (map (adjacentTransposition n) $ bubbleSort perm) == perm+--+bubbleSort :: Permutation -> [Int]+bubbleSort perm@(Permutation tgt) = runST action where+ (_,n) = bounds tgt++ action :: forall s. ST s [Int]+ action = do+ fwd <- newArray_ (1,n) :: ST s (STUArray s Int Int)+ inv <- newArray_ (1,n) :: ST s (STUArray s Int Int)+ forM_ [1..n] $ \i -> writeArray fwd i i+ forM_ [1..n] $ \i -> writeArray inv i i++ list <- forM [1..n] $ \x -> do++ let k = tgt ! x -- we take the number which will be at the @x@-th position at the end+ i <- readArray inv k -- number @k@ is at the moment at position @i@+ let j = x -- but the final place is at @x@ ++ let swaps = move i j+ forM_ swaps $ \y -> do++ a <- readArray fwd y+ b <- readArray fwd (y+1)+ writeArray fwd (y+1) a+ writeArray fwd y b++ u <- readArray inv a+ v <- readArray inv b+ writeArray inv b u+ writeArray inv a v++ return swaps+ + return (concat list)++ move :: Int -> Int -> [Int]+ move !i !j+ | j == i = []+ | j > i = [i..j-1]+ | j < i = [i-1,i-2..j]++-------------------------------------------------------------------------------- -- * Some concrete permutations -- | The permutation @[n,n-1,n-2,...,2,1]@. Note that it is the inverse of itself.@@ -710,6 +836,7 @@ -- The order is lexicographic. fasc2B_algorithm_L :: (Eq a, Ord a) => [a] -> [[a]] fasc2B_algorithm_L xs = unfold1 next (sort xs) where+ -- next :: [a] -> Maybe [a] next xs = case findj (reverse xs,[]) of Nothing -> Nothing@@ -725,170 +852,11 @@ findj ( [] , _ ) = Nothing -- inc :: a -> [a] -> ([a],[a]) -> [a]- inc u us ( (x:xs) , yys ) = if u >= x+ inc !u us ( (x:xs) , yys ) = if u >= x then inc u us ( xs , x : yys ) else reverse (x:us) ++ reverse (u:yys) ++ xs inc _ _ ( [] , _ ) = error "permute: should not happen" -------------------------------------------------------------------------------- -#ifdef QUICKCHECK--minPermSize = 1-maxPermSize = 123--newtype Elem = Elem Int deriving Eq-newtype Nat = Nat { fromNat :: Int } deriving (Eq,Ord,Show,Num,Random)--naturalSet :: Permutation -> Array Int Elem-naturalSet perm = listArray (1,n) [ Elem i | i<-[1..n] ] where- n = permutationSize perm--permInternalSet :: Permutation -> Array Int Elem-permInternalSet perm@(Permutation arr) = listArray (1,n) [ Elem (arr!i) | i<-[1..n] ] where- n = permutationSize perm--sameSize :: Permutation -> Permutation -> Bool-sameSize perm1 perm2 = ( permutationSize perm1 == permutationSize perm2)--newtype CyclicPermutation = Cyclic { fromCyclic :: Permutation } deriving Show--data SameSize = SameSize Permutation Permutation deriving Show--instance Random Permutation where- random g = randomPermutation size g1 where- (size,g1) = randomR (minPermSize,maxPermSize) g- randomR _ = random--instance Random CyclicPermutation where- random g = (Cyclic cycl,g2) where- (size,g1) = randomR (minPermSize,maxPermSize) g- (cycl,g2) = randomCyclicPermutation size g1- randomR _ = random--instance Random DisjointCycles where- random g = (disjcyc,g2) where- (size,g1) = randomR (minPermSize,maxPermSize) g- (perm,g2) = randomPermutation size g1- disjcyc = permutationToDisjointCycles perm- randomR _ = random--instance Random SameSize where- random g = (SameSize prm1 prm2, g3) where- (size,g1) = randomR (minPermSize,maxPermSize) g- (prm1,g2) = randomPermutation size g1 - (prm2,g3) = randomPermutation size g2- randomR _ = random--instance Arbitrary Nat where- arbitrary = choose (Nat 0 , Nat 50)- -instance Arbitrary Permutation where arbitrary = choose undefined-instance Arbitrary CyclicPermutation where arbitrary = choose undefined-instance Arbitrary DisjointCycles where arbitrary = choose undefined-instance Arbitrary SameSize where arbitrary = choose undefined---- | Runs all quickCheck tests-checkAll :: IO ()-checkAll = do- let f :: Testable a => a -> IO ()- f = quickCheck-- f prop_disjcyc_1- f prop_disjcyc_2 -- f prop_disjcyc_Mathematica-- f prop_randCyclic- f prop_inverse-- f prop_mulPerm- f prop_mulPermLeft- f prop_mulPermRight-- f prop_perm- f prop_permLeft- f prop_permRight- f prop_permLeftRight-- f prop_cycleLeft- f prop_cycleRight-- f prop_mulSign - f prop_invMul- f prop_cyclSign- f prop_permIsPerm- f prop_isEven- -prop_disjcyc_1 perm = ( perm == disjointCyclesToPermutation n (permutationToDisjointCycles perm) )- where n = permutationSize perm--prop_disjcyc_2 k dcyc = ( dcyc == permutationToDisjointCycles (disjointCyclesToPermutation n dcyc) )- where - n = fromNat k + m - m = case fromDisjointCycles dcyc of- [] -> 1- xxs -> maximum (concat xxs)---- PermutationCycles[ { 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 } ]--- Cycles [ {{1, 12, 18, 19, 14, 21, 23, 13, 22}, {2, 15, 8, 9, 20, 16, 10, 3, 5}, {4, 6, 7, 17}} ]-prop_disjcyc_Mathematica = (permutationToDisjointCycles perm == disjcyc) - && (disjointCyclesToPermutation n disjcyc == perm)- where- n = permutationSize perm- perm = toPermutation [ 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 ]- disjcyc = DisjointCycles [ [1, 12, 18, 19, 14, 21, 23, 13, 22], [2, 15, 8, 9, 20, 16, 10, 3, 5], [4, 6, 7, 17] ]--xperm = toPermutation [ 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 ]-xdisjcyc = DisjointCycles [ [1, 12, 18, 19, 14, 21, 23, 13, 22], [2, 15, 8, 9, 20, 16, 10, 3, 5], [4, 6, 7, 17] ]--prop_randCyclic cycl = ( isCyclicPermutation (fromCyclic cycl) )--prop_inverse perm = ( perm == inverse (inverse perm) ) --prop_mulPerm (SameSize perm1 perm2) = - ( permute perm2 (permute perm1 set) == permute (perm1 `multiply` perm2) set ) - where - set = naturalSet perm1--prop_mulPermRight (SameSize perm1 perm2) = - ( permuteRight perm2 (permuteRight perm1 set) == permuteRight (perm1 `multiply` perm2) set ) - where - set = naturalSet perm1--prop_mulPermLeft (SameSize perm1 perm2) = - ( permuteLeft perm2 (permuteLeft perm1 set) == permuteLeft (perm2 `multiply` perm1) set ) - where - set = naturalSet perm1--prop_perm perm = permute perm (naturalSet perm) == permInternalSet perm-prop_permLeft perm = permuteLeft perm (permInternalSet perm) == naturalSet perm-prop_permRight perm = permuteRight perm (naturalSet perm) == permInternalSet perm-prop_permLeftRight perm = permuteLeft (inverse perm) (naturalSet perm) == permuteRight (perm) (naturalSet perm) --prop_cycleLeft = permuteList (cycleLeft 5) "abcde" == "bcdea"-prop_cycleRight = permuteList (cycleRight 5) "abcde" == "eabcd"--prop_mulSign (SameSize perm1 perm2) = - ( sgn perm1 * sgn perm2 == sgn (perm1 `multiply` perm2) ) - where - sgn = signValue . signOfPermutation :: Permutation -> Int--prop_invMul (SameSize perm1 perm2) = - ( inverse perm2 `multiply` inverse perm1 == inverse (perm1 `multiply` perm2) ) --prop_cyclSign cycl = ( isEvenPermutation perm == odd n ) where- perm = fromCyclic cycl- n = permutationSize perm- -prop_permIsPerm perm = ( isPermutation (fromPermutation perm) ) --prop_isEven perm = ( isEvenPermutation perm == isEvenAlternative perm ) where- isEvenAlternative p = - even $ sum $ map (\x->x-1) $ map length $ fromDisjointCycles $ permutationToDisjointCycles p---#endif----------------------------------------------------------------------------------
Math/Combinat/Sign.hs view
@@ -7,11 +7,12 @@ -------------------------------------------------------------------------------- import Data.Monoid+import System.Random -------------------------------------------------------------------------------- data Sign- = Plus+ = Plus -- hmm, this way @Plus < Minus@, not sure about that | Minus deriving (Eq,Ord,Show,Read) @@ -19,6 +20,10 @@ mempty = Plus mappend = mulSign mconcat = productOfSigns++instance Random Sign where+ random g = let (b,g') = random g in (if b then Plus else Minus, g')+ randomR (u,v) g = let (y,g') = random g in (if u==v then u else y , g') isPlus, isMinus :: Sign -> Bool isPlus s = case s of { Plus -> True ; _ -> False }
Math/Combinat/Tableaux/LittlewoodRichardson.hs view
@@ -328,7 +328,7 @@ -- | Computes the expansion of the product of Schur polynomials @s[mu]*s[nu]@ using the -- Littlewood-Richardson rule. Note: this is symmetric in the two arguments. ----- Based on the wikipedia article <https://en.wikipedia.org/wiki/LittlewoodRichardson_rule>+-- Based on the wikipedia article <https://en.wikipedia.org/wiki/Littlewood-Richardson_rule> -- -- > lrMult mu nu == Map.fromList list where -- > lamw = weight nu + weight mu
Math/Combinat/Tableaux/Skew.hs view
@@ -221,35 +221,3 @@ -------------------------------------------------------------------------------- -#ifdef QUICKCHECK--prop_dual_dual :: SkewTableau Int -> Bool-prop_dual_dual st = (dualSkewTableau (dualSkewTableau st) == st)--prop_rowWord :: SkewTableau Int -> Bool-prop_rowWord st = (fillSkewPartitionWithRowWord shape content == st) where- shape = skewShape st- content = skewTableauRowWord st--prop_columnWord :: SkewTableau Int -> Bool-prop_columnWord st = (fillSkewPartitionWithColumnWord shape content == st) where- shape = skewShape st- content = skewTableauColumnWord st--prop_fill_shape :: SkewPartition -> Bool-prop_fill_shape shape = (shape == shape') where- tableau = fillSkewPartitionWithColumnWord shape [1..]- shape' = skewShape tableau--prop_semistandard :: SkewPartition -> Bool-prop_semistandard shape = and - [ isSemiStandardSkewTableau st - | n <- [1..nn] - , st <- semiStandardSkewTableaux n shape- ]- where- nn = skewPartitionWeight shape--#endif----------------------------------------------------------------------------------
+ Math/Combinat/TypeLevel.hs view
@@ -0,0 +1,117 @@++-- | Type-level hackery.+--+-- This module is used for groups whose parameters are encoded as type-level natural numbers,+-- for example finite cyclic groups, free groups, symmetric groups and braid groups.+--++{-# LANGUAGE PolyKinds, DataKinds, KindSignatures, ScopedTypeVariables, + ExistentialQuantification, Rank2Types #-}++module Math.Combinat.TypeLevel + ( -- * Proxy+ Proxy(..)+ , proxyUndef+ , proxyOf+ , proxyOf1+ , proxyOf2+ , asProxyTypeOf -- defined in Data.Proxy+ , asProxyTypeOf1+ -- * Type-level naturals as type arguments+ , typeArg + , iTypeArg+ -- * Hiding the type parameter+ , Some (..)+ , withSome , withSomeM+ , selectSome , selectSomeM+ , withSelected , withSelectedM+ )+ where++--------------------------------------------------------------------------------++import Data.Proxy ( Proxy(..) , asProxyTypeOf )+import GHC.TypeLits++import Math.Combinat.Numbers.Primes ( isProbablyPrime )++--------------------------------------------------------------------------------+-- * Proxy++proxyUndef :: Proxy a -> a+proxyUndef _ = error "proxyUndef"++proxyOf :: a -> Proxy a+proxyOf _ = Proxy++proxyOf1 :: f a -> Proxy a+proxyOf1 _ = Proxy++proxyOf2 :: g (f a) -> Proxy a+proxyOf2 _ = Proxy++asProxyTypeOf1 :: f a -> Proxy a -> f a +asProxyTypeOf1 y _ = y++--------------------------------------------------------------------------------+-- * Type-level naturals as type arguments++typeArg :: KnownNat n => f (n :: Nat) -> Integer+typeArg = natVal . proxyOf1++iTypeArg :: KnownNat n => f (n :: Nat) -> Int+iTypeArg = fromIntegral . typeArg++--------------------------------------------------------------------------------+-- * Hiding the type parameter++-- | Hide the type parameter of a functor. Example: @Some Braid@+data Some f = forall n. KnownNat n => Some (f n)++-- | Uses the value inside a 'Some'+withSome :: Some f -> (forall n. KnownNat n => f n -> a) -> a+withSome some polyFun = case some of { Some stuff -> polyFun stuff }++-- | Monadic version of 'withSome'+withSomeM :: Monad m => Some f -> (forall n. KnownNat n => f n -> m a) -> m a+withSomeM some polyAct = case some of { Some stuff -> polyAct stuff }++-- | Given a polymorphic value, we select at run time the+-- one specified by the second argument+selectSome :: Integral int => (forall n. KnownNat n => f n) -> int -> Some f+selectSome poly n = case someNatVal (fromIntegral n :: Integer) of+ Nothing -> error "selectSome: not a natural number"+ Just snat -> case snat of+ SomeNat pxy -> Some (asProxyTypeOf1 poly pxy)++-- | Monadic version of 'selectSome'+selectSomeM :: forall m f int. (Integral int, Monad m) => (forall n. KnownNat n => m (f n)) -> int -> m (Some f)+selectSomeM runPoly n = case someNatVal (fromIntegral n :: Integer) of+ Nothing -> error "selectSomeM: not a natural number"+ Just snat -> case snat of+ SomeNat pxy -> do+ poly <- runPoly + return $ Some (asProxyTypeOf1 poly pxy)++-- | Combination of 'selectSome' and 'withSome': we make a temporary structure+-- of the given size, but we immediately consume it.+withSelected + :: Integral int + => (forall n. KnownNat n => f n -> a) + -> (forall n. KnownNat n => f n) + -> int + -> a+withSelected f x n = withSome (selectSome x n) f++-- | (Half-)monadic version of 'withSelected'+withSelectedM + :: forall m f int a. (Integral int, Monad m) + => (forall n. KnownNat n => f n -> a) + -> (forall n. KnownNat n => m (f n)) + -> int + -> m a+withSelectedM f mx n = do + s <- selectSomeM mx n+ return (withSome s f)++--------------------------------------------------------------------------------
combinat.cabal view
@@ -1,5 +1,5 @@ Name: combinat-Version: 0.2.8.0+Version: 0.2.8.1 Synopsis: Generate and manipulate various combinatorial objects. Description: A collection of functions to generate, count, manipulate and visualize all kinds of combinatorial objects like @@ -22,18 +22,11 @@ extra-source-files: svg/*.svg svg/src/gen_figures.hs -Flag withQuickCheck- Description: Compile with the QuickCheck tests. - default: False- + Library Build-Depends: base >= 4 && < 5, array >= 0.5, containers, random, transformers - if flag(withQuickCheck)- Build-Depends: QuickCheck- cpp-options: -DQUICKCHECK- Exposed-Modules: Math.Combinat Math.Combinat.Classes Math.Combinat.Numbers@@ -68,6 +61,7 @@ Math.Combinat.LatticePaths Math.Combinat.ASCII Math.Combinat.Helper+ Math.Combinat.TypeLevel Default-Extensions: CPP, BangPatterns Other-Extensions: MultiParamTypeClasses, ScopedTypeVariables, @@ -80,3 +74,17 @@ ghc-options: -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports ++test-suite combinat-tests+ + type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: TestSuite.hs++ build-depends: base >= 4 && < 5, array >= 0.5, containers, random, transformers,+ combinat,+ QuickCheck >= 2, test-framework, test-framework-quickcheck2++ Default-Language: Haskell2010+ Default-Extensions: CPP, BangPatterns+
+ test/TestSuite.hs view
@@ -0,0 +1,41 @@++module Main where++--------------------------------------------------------------------------------++import Test.Framework+import Test.Framework.Providers.QuickCheck2++import Tests.Permutations ( testgroup_Permutations )+import Tests.Partitions.Integer ( testgroup_IntegerPartitions )+import Tests.Partitions.Skew ( testgroup_SkewPartitions )+import Tests.Braid ( testgroup_Braid + , testgroup_Braid_NF )+import Tests.Series ( testgroup_PowerSeries )+import Tests.SkewTableaux ( testgroup_SkewTableaux )+import Tests.Thompson ( testgroup_ThompsonF )+import Tests.LatticePaths ( testgroup_LatticePaths )++--------------------------------------------------------------------------------++main :: IO ()+main = defaultMain tests++tests :: [Test]+tests = + [ testgroup_Permutations+ , testGroup "Partitions" + [ testgroup_IntegerPartitions+ , testgroup_SkewPartitions+ ]+ , testgroup_SkewTableaux+ , testgroup_ThompsonF+ , testgroup_LatticePaths+ , testGroup "Braids" + [ testgroup_Braid + , testgroup_Braid_NF + ]+ , testgroup_PowerSeries + ]++--------------------------------------------------------------------------------