diff --git a/Math/Combinat/Groups/Braid.hs b/Math/Combinat/Groups/Braid.hs --- a/Math/Combinat/Groups/Braid.hs +++ b/Math/Combinat/Groups/Braid.hs @@ -48,7 +48,7 @@ import Math.Combinat.TypeLevel import Math.Combinat.Numbers.Series -import Math.Combinat.Permutations ( Permutation(..) ) +import Math.Combinat.Permutations ( Permutation(..) , (!!!) ) import qualified Math.Combinat.Permutations as P -------------------------------------------------------------------------------- @@ -233,8 +233,8 @@ -- | 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 +tauPerm perm = P.toPermutationUnsafeN n [ (n+1) - perm !!! (n-i) | i<-[0..n-1] ] where + n = P.permutationSize perm -------------------------------------------------------------------------------- -- * Group operations @@ -265,7 +265,7 @@ -- | A braid is pure if its permutation is trivial isPureBraid :: KnownNat n => Braid n -> Bool -isPureBraid braid = (braidPermutation braid == P.identity n) where +isPureBraid braid = (braidPermutation braid == P.identityPermutation n) where n = numberOfStrands braid -- | Returns the left-to-right permutation associated to the braid. @@ -282,13 +282,13 @@ -- we got the two-line notation of the permutation. -- braidPermutation :: KnownNat n => Braid n -> Permutation -braidPermutation braid@ (Braid gens) = perm where +braidPermutation braid@(Braid gens) = perm where n = numberOfStrands braid perm = _braidPermutation n (map brGenIdx gens) -- | This is an untyped version of 'braidPermutation' _braidPermutation :: Int -> [Int] -> Permutation -_braidPermutation n idxs = Permutation (runSTUArray action) where +_braidPermutation n idxs = P.uarrayToPermutationUnsafe (runSTUArray action) where action :: forall s. ST s (STUArray s Int Int) action = do @@ -358,8 +358,8 @@ -- | Returns the individual \"phases\" of the a permutation braid realizing the -- given permutation. _permutationBraid' :: Permutation -> [[Int]] -_permutationBraid' perm@(Permutation arr) = runST action where - (1,n) = bounds arr +_permutationBraid' perm = runST action where + n = P.permutationSize perm action :: forall s. ST s [[Int]] action = do @@ -388,7 +388,7 @@ let worker phase | phase >= n = return [] | otherwise = do - let tgt = (arr ! phase) + let tgt = P.lookupPermutation perm phase -- (arr ! phase) src <- readArray cfwd tgt let this = [src-1,src-2..phase] mapM_ doSwap $ this diff --git a/Math/Combinat/Groups/Braid/NF.hs b/Math/Combinat/Groups/Braid/NF.hs --- a/Math/Combinat/Groups/Braid/NF.hs +++ b/Math/Combinat/Groups/Braid/NF.hs @@ -53,7 +53,7 @@ import Math.Combinat.Helper import Math.Combinat.Sign -import Math.Combinat.Permutations ( Permutation(..) , isIdentityPermutation , isReversePermutation ) +import Math.Combinat.Permutations ( Permutation(..) , (!!!) , isIdentityPermutation , isReversePermutation ) import qualified Math.Combinat.Permutations as P import Math.Combinat.Groups.Braid @@ -128,7 +128,7 @@ facs = leftGreedyFactors n $ map brGenIdx gs replaceFac idxs = XDelta (-1) : map XSigma (_permutationBraid perm) where - perm = (P.reversePermutation n) `P.multiply` (P.adjacentTranspositions n idxs) + perm = (P.reversePermutation n) `P.multiplyPermutation` (P.adjacentTranspositions n idxs) -- | Replaces @sigma_i^-1@ generators by @(Delta^-1 * L_i)@. @@ -259,13 +259,13 @@ -- | The permutation @posL k :: Braid n@ is realizing posLPerm :: Int -> Int -> Permutation posLPerm n k - | k>0 && k0 && k Int -> Permutation posRPerm n k - | k>0 && k0 && k permutationStartingSet p == permWordStartingSet n (_permutationBraid p) -- permutationStartingSet :: Permutation -> [Int] -permutationStartingSet = permutationFinishingSet . P.inverse +permutationStartingSet = permutationFinishingSet . P.inversePermutation -- | This satisfies -- -- > permutationFinishingSet p == permWordFinishingSet n (_permutationBraid p) -- permutationFinishingSet :: Permutation -> [Int] -permutationFinishingSet (Permutation arr) - = [ i | i<-[1..n-1] , arr ! i > arr ! (i+1) ] where (1,n) = bounds arr +permutationFinishingSet perm + = [ i | i<-[1..n-1] , perm !!! i > perm !!! (i+1) ] where n = P.permutationSize perm -- | Returns the list of permutations failing Lemma 2.5 in [2] -- (so an empty list means the implementaton is correct) @@ -417,8 +417,8 @@ [] -> let (e',rs) = worker' eq 0 e rest in (e', preal : rs) (j:_) -> worker' (-e) (-e) e (p':q':rest') where s = P.adjacentTransposition n j - p' = P.multiply s preal - q' = P.multiply qreal s + p' = P.multiplyPermutation s preal + q' = P.multiplyPermutation qreal s where preal = oddTau (e+ep) p -- the "real" p qreal = oddTau (e+eq) q -- the "real" q diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs --- a/Math/Combinat/Helper.hs +++ b/Math/Combinat/Helper.hs @@ -1,5 +1,5 @@ --- | Miscellaneous helper functions +-- | Miscellaneous helper functions used internally {-# LANGUAGE BangPatterns, PolyKinds, GeneralizedNewtypeDeriving #-} module Math.Combinat.Helper where @@ -57,6 +57,46 @@ interleave [x] [] = x : [] interleave [] [] = [] interleave _ _ = error "interleave: shouldn't happen" + +evens, odds :: [a] -> [a] +evens (x:xs) = x : odds xs +evens [] = [] +odds (x:xs) = evens xs +odds [] = [] + +-------------------------------------------------------------------------------- +-- * multiplication + +-- | Product of list of integers, but in interleaved order (for a list of big numbers, +-- it should be faster than the linear order) +productInterleaved :: [Integer] -> Integer +productInterleaved = go where + go [] = 1 + go [x] = x + go [x,y] = x*y + go list = go (evens list) * go (odds list) + +-- | Faster implementation of @product [ i | i <- [a+1..b] ]@ +productFromTo :: Integral a => a -> a -> Integer +productFromTo = go where + go !a !b + | dif < 1 = 1 + | dif < 5 = product [ fromIntegral i | i<-[a+1..b] ] + | otherwise = go a half * go half b + where + dif = b - a + half = div (a+b+1) 2 + +-- | Faster implementation of product @[ i | i <- [a+1,a+3,..b] ]@ +productFromToStride2 :: Integral a => a -> a -> Integer +productFromToStride2 = go where + go !a !b + | dif < 1 = 1 + | dif < 9 = product [ fromIntegral i | i<-[a+1,a+3..b] ] + | otherwise = go a half * go half b + where + dif = b - a + half = a + 2*(div dif 4) -------------------------------------------------------------------------------- -- * equality and ordering diff --git a/Math/Combinat/Numbers.hs b/Math/Combinat/Numbers.hs --- a/Math/Combinat/Numbers.hs +++ b/Math/Combinat/Numbers.hs @@ -1,9 +1,12 @@ module Math.Combinat.Numbers - ( module Math.Combinat.Numbers.Sequences - , module Math.Combinat.Numbers.Integers + ( module Math.Combinat.Numbers.Integers + , module Math.Combinat.Numbers.Primes + , module Math.Combinat.Numbers.Sequences ) where -import Math.Combinat.Numbers.Sequences import Math.Combinat.Numbers.Integers +import Math.Combinat.Numbers.Primes +import Math.Combinat.Numbers.Sequences + diff --git a/Math/Combinat/Numbers/Primes.hs b/Math/Combinat/Numbers/Primes.hs --- a/Math/Combinat/Numbers/Primes.hs +++ b/Math/Combinat/Numbers/Primes.hs @@ -2,13 +2,20 @@ -- | Prime numbers and related number theoretical stuff. module Math.Combinat.Numbers.Primes - ( -- * List of prime numbers - primes + ( -- * Elementary number theory + divides + , divisors, squareFreeDivisors, squareFreeDivisors_ + , divisorSum , divisorSum' + , moebiusMu , eulerTotient , liouvilleLambda + -- * List of prime numbers + , primes , primesSimple , primesTMWE -- * Prime factorization - , groupIntegerFactors + , factorize, factorizeNaive + , productOfFactors , integerFactorsTrialDivision + , groupIntegerFactors -- * Modulo @m@ arithmetic , powerMod -- * Prime testing @@ -20,14 +27,93 @@ -------------------------------------------------------------------------------- +import Data.List ( group , sort , foldl' ) + +import Math.Combinat.Sign +import Math.Combinat.Helper import Math.Combinat.Numbers.Integers -import Data.List ( group , sort ) +-- import Math.Combinat.Sets ( sublists ) -- cyclic dependency... +import Math.Combinat.Tuples ( tuples' ) + import Data.Bits import System.Random -------------------------------------------------------------------------------- + +-- | @d `divides` n@ +divides :: Integer -> Integer -> Bool +divides d n = (mod n d == 0) + +{-# SPECIALIZE moebiusMu :: Int -> Int #-} +{-# SPECIALIZE moebiusMu :: Integer -> Integer #-} +-- | The Moebius mu function +moebiusMu :: (Integral a, Num b) => a -> b +moebiusMu n + | any (>1) expos = 0 + | even (length primes) = 1 + | otherwise = -1 + where + factors = groupIntegerFactors $ integerFactorsTrialDivision $ fromIntegral n + (primes,expos) = unzip factors + +{-# SPECIALIZE liouvilleLambda :: Int -> Int #-} +{-# SPECIALIZE liouvilleLambda :: Integer -> Integer #-} +-- | The Liouville lambda function +liouvilleLambda :: (Integral a, Num b) => a -> b +liouvilleLambda n = + if odd (foldl' (+) 0 $ map snd grps) + then -1 + else 1 + where + grps = groupIntegerFactors $ integerFactorsTrialDivision $ fromIntegral n + +-- | Sum ofthe of the divisors +divisorSum :: Integer -> Integer +divisorSum n = foldl' (+) 0 [ d | d <- divisors n] + +-- | Sum of @k@-th powers of the divisors +divisorSum' :: Int -> Integer -> Integer +divisorSum' k n = foldl' (+) 0 [ d^k | d <- divisors n] + +-- | Euler's totient function +eulerTotient :: Integer -> Integer +eulerTotient n = div n prodp * prodp1 where + grps = groupIntegerFactors $ integerFactorsTrialDivision n + ps = map fst grps + prodp = foldl' (*) 1 [ p | p <- ps ] + prodp1 = foldl' (*) 1 [ p-1 | p <- ps ] + +-- | Divisors of @n@ (note: the result is /not/ ordered!) +divisors :: Integer -> [Integer] +divisors n = [ f tup | tup <- tuples' expos ] where + grps = groupIntegerFactors $ integerFactorsTrialDivision n + (ps,expos) = unzip grps + f es = foldl' (*) 1 $ zipWith (^) ps es + +-- | List of square-free divisors together with their Mobius mu value +-- (note: the result is /not/ ordered!) +squareFreeDivisors :: Integer -> [(Integer,Sign)] +squareFreeDivisors n = map f (sublists primes) where + grps = groupIntegerFactors $ integerFactorsTrialDivision n + primes = map fst grps + f ps = ( foldl' (*) 1 ps , if even (length ps) then Plus else Minus) + +-- | List of square-free divisors +-- (note: the result is /not/ ordered!) +squareFreeDivisors_ :: Integer -> [Integer] +squareFreeDivisors_ n = map f (sublists primes) where + grps = groupIntegerFactors $ integerFactorsTrialDivision n + primes = map fst grps + f ps = foldl' (*) 1 ps + +-- | To avoid cyclic dependencies, I made a local copy of this... +sublists :: [a] -> [[a]] +sublists [] = [[]] +sublists (x:xs) = sublists xs ++ map (x:) (sublists xs) + +-------------------------------------------------------------------------------- -- List of prime numbers -- | Infinite list of primes, using the TMWE algorithm. @@ -54,10 +140,10 @@ pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t wheel = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2: 4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel - gaps k ws@(w:t) cs@ ~(c:u) + gaps k ws@(w:t) cs@(~(c:u)) | k==c = gaps (k+w) t u | True = k : gaps (k+w) t cs - roll k ws@(w:t) ps@ ~(p:u) + roll k ws@(w:t) ps@(~(p:u)) | k==p = scanl (\c d->c+p*d) (p*p) ws : roll (k+w) t u | True = roll (k+w) t ps @@ -77,6 +163,22 @@ -------------------------------------------------------------------------------- -- Prime factorization + +factorize :: Integer -> [(Integer,Int)] +factorize = factorizeNaive + +factorizeNaive :: Integer -> [(Integer,Int)] +factorizeNaive = groupIntegerFactors . integerFactorsTrialDivision + +productOfFactors :: [(Integer,Int)] -> Integer +productOfFactors = productInterleaved . map (uncurry pow) where + pow _ 0 = 1 + pow p 1 = p + pow 2 n = shiftL 1 n + pow p 2 = p*p + pow p n = if even n + then (pow p (shiftR n 1))^2 + else p * (pow p (shiftR n 1))^2 -- | Groups integer factors. Example: from [2,2,2,3,3,5] we produce [(2,3),(3,2),(5,1)] groupIntegerFactors :: [Integer] -> [(Integer,Int)] diff --git a/Math/Combinat/Numbers/Sequences.hs b/Math/Combinat/Numbers/Sequences.hs --- a/Math/Combinat/Numbers/Sequences.hs +++ b/Math/Combinat/Numbers/Sequences.hs @@ -4,39 +4,148 @@ -- See the \"On-Line Encyclopedia of Integer Sequences\", -- . +{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} module Math.Combinat.Numbers.Sequences where -------------------------------------------------------------------------------- import Data.Array +import Data.Bits ( shiftL , shiftR , (.&.) ) -import Math.Combinat.Helper ( sum' ) +import Math.Combinat.Helper import Math.Combinat.Sign +import Math.Combinat.Numbers.Primes ( primes , factorize , productOfFactors ) + +import qualified Data.Map.Strict as Map -- used for factorialPrimeExponentsNaive + -------------------------------------------------------------------------------- -- * Factorial --- | A000142. +-- | The factorial function (A000142). factorial :: Integral a => a -> Integer -factorial n - | n < 0 = error "factorial: input should be nonnegative" +factorial = factorialSplit + +-- | Faster implementation of the factorial function +factorialSplit :: Integral a => a -> Integer +factorialSplit n = productFromTo 1 n + +-- | Naive implementation of factorial +factorialNaive :: Integral a => a -> Integer +factorialNaive n + | n < 0 = error "factorialNaive: input should be nonnegative" | n == 0 = 1 | otherwise = product [1..fromIntegral n] --- | A006882. +-- | \"Swing factorial\" algorithm +factorialSwing :: Integral a => a -> Integer +factorialSwing n = productOfFactors (factorialPrimeExponents $ fromIntegral n) where + +-------------------------------------------------------------------------------- + +-- | Prime factorization of the factorial (using the \"swing factorial\" algorithm) +factorialPrimeExponents :: Int -> [(Integer,Int)] +factorialPrimeExponents n = filter cond $ zip primes (factorialPrimeExponents_ n) where + cond (_,!e) = e > 0 + +factorialPrimeExponentsNaive :: forall a. Integral a => a -> [(Integer,Int)] +factorialPrimeExponentsNaive n = result where + fi = fromIntegral :: a -> Integer + result = Map.toList + $ Map.unionsWith (+) + $ map Map.fromList + $ map factorize + $ map fi [1..n] + +factorialPrimeExponents_ :: Int -> [Int] +factorialPrimeExponents_ = go where + go 0 = [] + go 1 = [] + go 2 = [1] + go !n = longAdd half swing where + half = map (flip shiftL 1) $ go (shiftR n 1) + swing = swingFactorialExponents_ n + + longAdd :: [Int] -> [Int] -> [Int] + longAdd xs [] = xs + longAdd [] ys = ys + longAdd (!x:xs) (!y:ys) = (x+y) : longAdd xs ys + +-- | Prime factorizaiton of the \"swing factorial\" function) +swingFactorialExponents_ :: Int -> [Int] +swingFactorialExponents_ = go where + go 0 = [] + go 1 = [] + go 2 = [1] + go n = expo2 : map expo (tail ps) where + + nn = fromIntegral n :: Integer + + ps :: [Integer] + ps = takeWhile (<=nn) primes + + expo2 :: Int + expo2 = go 0 (shiftR n 1) where + go :: Int -> Int -> Int + go !acc !r + | r < 1 = acc + | otherwise = go acc' r' + where + acc' = acc + (r .&. 1) + r' = shiftR r 1 + + expo :: Integer -> Int + expo pp = go 0 (div n p) where + p = fromInteger pp :: Int + go :: Int -> Int -> Int + go !acc !r + | r < 1 = acc + | otherwise = go acc' r' + where + acc' = acc + (r .&. 1) + r' = div r p + +-------------------------------------------------------------------------------- + +-- | The double factorial doubleFactorial :: Integral a => a -> Integer -doubleFactorial n - | n < 0 = error "doubleFactorial: input should be nonnegative" +doubleFactorial = doubleFactorialSplit + +-- | Faster implementation of the double factorial function +doubleFactorialSplit :: Integral a => a -> Integer +doubleFactorialSplit n + | n < 0 = error "doubleFactorialSplit: input should be nonnegative" | n == 0 = 1 + | odd n = productFromToStride2 2 n + | otherwise = let halfn = div n 2 + in shiftL (factorialSplit halfn) (fromIntegral halfn) + +-- | Naive implementation of the double factorial (A006882). +doubleFactorialNaive :: Integral a => a -> Integer +doubleFactorialNaive n + | n < 0 = error "doubleFactorialNaive: input should be nonnegative" + | n == 0 = 1 | odd n = product [1,3..fromIntegral n] | otherwise = product [2,4..fromIntegral n] -------------------------------------------------------------------------------- -- * Binomial and multinomial --- | A007318. Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below. +-- | Binomial numbers (A007318). Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below. binomial :: Integral a => a -> a -> Integer -binomial n k +binomial = binomialSplit + +-- | Faster implementation of binomial +binomialSplit :: Integral a => a -> a -> Integer +binomialSplit n k + | k > n = 0 + | k < 0 = 0 + | k > (n `div` 2) = binomialSplit n (n-k) + | otherwise = (productFromTo (n-k) n) `div` (productFromTo 1 k) + +-- | A007318. Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below. +binomialNaive :: Integral a => a -> a -> Integer +binomialNaive n k | k > n = 0 | k < 0 = 0 | k > (n `div` 2) = binomial n (n-k) diff --git a/Math/Combinat/Partitions/Integer.hs b/Math/Combinat/Partitions/Integer.hs --- a/Math/Combinat/Partitions/Integer.hs +++ b/Math/Combinat/Partitions/Integer.hs @@ -29,6 +29,10 @@ , toPartition , toPartitionUnsafe , isPartition + -- * Conversion to\/from exponent vectors + , toExponentVector + , fromExponentVector + , dropTailingZeros -- * Union and sum , unionOfPartitions , sumOfPartitions @@ -49,8 +53,12 @@ , randomPartition , randomPartitions -- * Dominating \/ dominated partitions + , dominanceCompare , dominatedPartitions , dominatingPartitions + -- * Conjugate lexicographic ordering + , conjugateLexicographicCompare + , ConjLex (..) , fromConjLex -- * Partitions with given number of parts , partitionsWithKParts -- * Partitions with only odd\/distinct parts @@ -83,7 +91,7 @@ import Data.Array import System.Random -import Math.Combinat.Partitions.Integer.Naive +import Math.Combinat.Partitions.Integer.Naive hiding () -- this is for haddock! import Math.Combinat.Partitions.Integer.IntList import Math.Combinat.Partitions.Integer.Count @@ -116,6 +124,41 @@ isPartition (x:xs@(y:_)) = (x >= y) && isPartition xs -------------------------------------------------------------------------------- +-- * Conversion to\/from exponent vectors + +-- | Converts a partition to an exponent vector. +-- +-- For example, +-- +-- > toExponentVector (Partition [4,4,2,2,2,1]) == [1,3,0,2] +-- +-- meaning @(1^1,2^3,3^0,4^2)@. +-- +toExponentVector :: Partition -> [Int] +toExponentVector part = fun 1 $ reverse $ group (fromPartition part) where + fun _ [] = [] + fun !k gs@(this@(i:_):rest) + | k < i = replicate (i-k) 0 ++ fun i gs + | otherwise = length this : fun (k+1) rest + +fromExponentVector :: [Int] -> Partition +fromExponentVector expos = Partition $ concat $ reverse $ zipWith f [1..] expos where + f !i !e = replicate e i + +dropTailingZeros :: [Int] -> [Int] +dropTailingZeros = reverse . dropWhile (==0) . reverse + +{- +-- alternative implementation +toExponentialVector2 :: Partition -> [Int] +toExponentialVector2 p = go 1 (toExponentialForm p) where + go _ [] = [] + go !i ef@((j,e):rest) = if i Partition -> Maybe Ordering +dominanceCompare p q + | p==q = Just EQ + | p `dominates` q = Just GT + | q `dominates` p = Just LT + | otherwise = Nothing + -- | Lists all partitions of the same weight as @lambda@ and also dominated by @lambda@ -- (that is, all partial sums are less or equal): -- @@ -247,6 +298,28 @@ -- dominatingPartitions :: Partition -> [Partition] dominatingPartitions (Partition_ mu) = map Partition_ (_dominatingPartitions mu) + +-------------------------------------------------------------------------------- +-- * Conjugate lexicographic ordering + +conjugateLexicographicCompare :: Partition -> Partition -> Ordering +conjugateLexicographicCompare p q = compare (dualPartition q) (dualPartition p) + +newtype ConjLex = ConjLex Partition deriving (Eq,Show) + +fromConjLex :: ConjLex -> Partition +fromConjLex (ConjLex p) = p + +instance Ord ConjLex where + compare (ConjLex p) (ConjLex q) = conjugateLexicographicCompare p q + +-- {- CONJUGATE LEXICOGRAPHIC ordering is a refinement of dominance partial ordering -} +-- let test n = [ ConjLex p >= ConjLex q | p <- partitions n , q <-partitions n , p `dominates` q ] +-- and (test 20) + +-- {- LEXICOGRAPHIC ordering is a refinement of dominance partial ordering -} +-- let test n = [ p >= q | p <- partitions n , q <-partitions n , p `dominates` q ] +-- and (test 20) -------------------------------------------------------------------------------- -- * Partitions with given number of parts diff --git a/Math/Combinat/Partitions/Integer/Compact.hs b/Math/Combinat/Partitions/Integer/Compact.hs --- a/Math/Combinat/Partitions/Integer/Compact.hs +++ b/Math/Combinat/Partitions/Integer/Compact.hs @@ -3,25 +3,17 @@ Partitions are conceptually nonincreasing sequences of /positive/ integers. -When the partition fits into a 15x15 rectangle, we encode the parts as nibbles in a single 64-bit word. -The most significant nibble is the first element, and the least significant nibble is used to encode -the length. This way equality and comparison of 64-bit words is the same as the corresponding operations -for partitions (lexicographic ordering). - -This will make working with small partitions much more memory efficient (very helpful when -building tables indexed by partitions, for example!) and hopefully quite a bit faster, too. - -When they do not fit into a 15x15 rectangle, but fit into 255x7, 255x15, 255x23 or 255x31, respectively, -then we extend the above to use the bytes of 1, 2, 3 or 4 64-bit words. - -In the general case, we encode the partition as a list of 64-bit words, each encoding 4 16-bit parts. - -Partitions with elements bigger than 65535 are not supported. +This implementation uses the @compact-word-vectors@ library internally to provide +a much more memory-efficient Partition type that the naive lists of integer. +This is very helpful when building large tables indexed by partitions, for example; +and hopefully quite a bit faster, too. Note: This is an internal module, you are not supposed to import it directly. +It is also not fully ready to be used yet... + -} -{-# LANGUAGE BangPatterns, PatternSynonyms, ViewPatterns, ForeignFunctionInterface #-} +{-# LANGUAGE BangPatterns, PatternSynonyms, ViewPatterns #-} module Math.Combinat.Partitions.Integer.Compact where -------------------------------------------------------------------------------- @@ -31,41 +23,30 @@ import Data.Ord import Data.List ( intercalate , group , sort , sortBy , foldl' , scanl' ) -import Math.Combinat.Compositions ( compositions' ) +import Data.Vector.Compact.WordVec ( WordVec , Shape(..) ) +import qualified Data.Vector.Compact.WordVec as V +import Math.Combinat.Compositions ( compositions' ) -------------------------------------------------------------------------------- -- * The compact partition data type -data Partition - = Nibble {-# UNPACK #-} !Word64 - | Medium1 {-# UNPACK #-} !Word64 - | Medium2 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 - | Medium3 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 - | Medium4 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 - | WordList {-# UNPACK #-} !Int ![Word64] - deriving (Eq,Show) - --------------------------------------------------------------------------------- - --- | for debugging -partitionPrefixChar :: Partition -> Char -partitionPrefixChar p = case p of - Nibble {} -> 'N' - Medium1 {} -> '1' - Medium2 {} -> '2' - Medium3 {} -> '3' - Medium4 {} -> '4' - WordList {} -> 'L' +newtype Partition + = Partition WordVec + deriving Eq -{- instance Show Partition where - show compact = partitionPrefixChar compact - : '<' : intercalate "," (map show $ toList compact) ++ ">" --} + showsPrec = showsPrecPartition +showsPrecPartition :: Int -> Partition -> ShowS +showsPrecPartition prec (Partition vec) + = showParen (prec > 10) + $ showString "Partition" + . showChar ' ' + . shows (V.toList vec) + instance Ord Partition where - compare = cmp + compare = cmpLexico -------------------------------------------------------------------------------- -- * Pattern synonyms @@ -97,200 +78,60 @@ -- * Lexicographic comparison -- | The lexicographic ordering -cmp :: Partition -> Partition -> Ordering -cmp (Nibble a) (Nibble b) = compare a b -cmp (Medium1 a1) (Medium1 b1) = compare a1 b1 -cmp (Medium2 a1 a2) (Medium2 b1 b2) = compare (a1,a2) (b1,b2) -cmp (Medium3 a1 a2 a3) (Medium3 b1 b2 b3) = compare (a1,a2,a3) (b1,b2,b3) -cmp (Medium4 a1 a2 a3 a4) (Medium4 b1 b2 b3 b4) = compare (a1,a2,a3,a4) (b1,b2,b3,b4) -cmp (WordList _ as) (WordList _ bs) = compare as bs -cmp p q = compare (toList p) (toList q) - +cmpLexico :: Partition -> Partition -> Ordering +cmpLexico (Partition vec1) (Partition vec2) = compare (V.toList vec1) (V.toList vec2) + -------------------------------------------------------------------------------- -- * Basic (de)constructrion empty :: Partition -empty = Nibble 0 +empty = Partition (V.empty) isEmpty :: Partition -> Bool -isEmpty compact = case compact of - Nibble x -> (x == 0) - _ -> False +isEmpty (Partition vec) = V.null vec -------------------------------------------------------------------------------- singleton :: Int -> Partition -singleton x - | x == 0 = Nibble 0 - | x <= 15 = Nibble $ shiftL (i2w x) 60 + 1 - | x <= 255 = Medium1 $ shiftL (i2w x) 56 + 1 - | x <= 65535 = WordList 1 [ shiftL (i2w x) 48 ] - | otherwise = error "singleton: partitions with elements bigger than 65535 are not supported" +singleton x + | x > 0 = Partition (V.singleton $ i2w x) + | x == 0 = empty + | otherwise = error "Parittion/singleton: negative input" -------------------------------------------------------------------------------- uncons :: Partition -> Maybe (Int,Partition) -uncons compact = case compact of - Nibble 0 -> Nothing - Nibble w -> Just ( w2i (shiftR w 60) , Nibble $ shiftL (w .&. 0x0ffffffffffffff0) 4 + ((w .&. 15) - 1) ) - Medium1 w1 -> Just ( w2i (shiftR w1 56) , partitionTail compact ) - Medium2 w1 w2 -> Just ( w2i (shiftR w1 56) , partitionTail compact ) - Medium3 w1 w2 w3 -> Just ( w2i (shiftR w1 56) , partitionTail compact ) - Medium4 w1 w2 w3 w4 -> Just ( w2i (shiftR w1 56) , partitionTail compact ) - WordList n (w:rest) -> Just ( w2i (shiftR w 48) , partitionTail compact ) - --------------------------------------------------------------------------------- +uncons (Partition vec) = case V.uncons vec of + Nothing -> Nothing + Just (h,tl) -> Just (w2i h, Partition tl) -- | @partitionTail p == snd (uncons p)@ partitionTail :: Partition -> Partition -partitionTail compact = case compact of - - Nibble 0 -> Nibble 0 - Nibble w -> Nibble $ shiftL (w .&. 0x0ffffffffffffff0) 4 + ((w .&. 15) - 1) - - Medium1 w1 -> - let !y = (shiftR w1 48) .&. 255 -- next element - !n = w1 .&. 15 - in if y <= 15 - then makeNibble (w2i $ n-1) $ safeTail $ toList compact - else Medium1 $ shiftL (w1 .&. 0x00ffffffffffff00) 8 + (n-1) - - Medium2 w1 w2 -> - let !y = (shiftR w1 48) .&. 255 - !n = w2 .&. 255 - in if y <= 15 - then makeNibble (w2i $ n-1) $ safeTail $ toList compact - else if n <= 8 - then Medium1 $ shiftL (w1 .&. 0x00ffffffffffffff) 8 + shiftL (shiftR w2 56) 8 + (n-1) - - else Medium2 ( shiftL w1 8 + shiftR w2 56 ) - ( shiftL (w2 .&. 0x00ffffffffffff00) 8 + (n-1) ) - - Medium3 w1 w2 w3 -> - let !y = (shiftR w1 48) .&. 255 - !n = w3 .&. 255 - in if y <= 15 && n <= 16 - then makeNibble (w2i $ n-1) $ safeTail $ toList compact - else if n <= 16 - then Medium2 ( shiftL w1 8 + shiftR w2 56 ) - ( shiftL w2 8 + shiftR w3 56 + shiftL (shiftR w3 56) 8 + (n-1) ) - - else Medium3 ( shiftL w1 8 + shiftR w2 56 ) - ( shiftL w2 8 + shiftR w3 56 ) - ( shiftL (w3 .&. 0x00ffffffffffff00) 8 + (n-1) ) - - _ -> - let n = width compact - in fromDescList' (n-1) $ safeTail $ toList compact +partitionTail (Partition vec) = Partition (V.tail vec) --------------------------------------------------------------------------------- +------------------------------------------------------------------------------- -- | We assume that @x >= partitionHeight p@! cons :: Int -> Partition -> Partition -cons !x !compact = case compact of - - Nibble 0 -> singleton x - - Nibble word - | x <= 15 && n < 15 -> Nibble $ shiftR word 4 + shiftL xw 60 + (n+1) - | x <= 255 -> makeMedium (w2i $ n+1) (x : toList compact) - | otherwise -> makeWordList (w2i $ n+1) (x : toList compact) - where - n = word .&. 15 - xw = i2w x - - Medium1 w1 - | x <= 255 && n < 7 -> Medium1 (shiftR w1 8 + shiftL xw 56 + (n+1)) - | x <= 255 -> Medium2 (shiftR w1 8 + shiftL xw 56 ) 8 - | otherwise -> makeWordList (w2i $ n+1) (x : toList compact) - where - n = w1 .&. 255 - xw = i2w x - - Medium2 w1 w2 - | x <= 255 && n < 15 -> Medium2 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56 + (n+1)) - | x <= 255 -> Medium3 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56 ) 16 - | otherwise -> makeWordList (w2i $ n+1) (x : toList compact) - where - n = w2 .&. 255 - xw = i2w x - - Medium3 w1 w2 w3 - | x <= 255 && n < 23 -> Medium3 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56) (shiftR w3 8 + shiftL (w2 .&. 255) 56 + (n+1)) - | x <= 255 -> Medium4 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56) (shiftR w3 8 + shiftL (w2 .&. 255) 56 ) 24 - | otherwise -> makeWordList (w2i $ n+1) (x : toList compact) - where - n = w3 .&. 255 - xw = i2w x - - Medium4 w1 w2 w3 w4 - | x <= 255 && n < 31 -> Medium4 (shiftR w1 8 + shiftL xw 56) - (shiftR w2 8 + shiftL (w1 .&. 255) 56) - (shiftR w3 8 + shiftL (w2 .&. 255) 56) - (shiftR w4 8 + shiftL (w3 .&. 255) 56 + (n+1)) - | otherwise -> makeWordList (w2i $ n+1) (x : toList compact) - where - n = w4 .&. 255 - xw = i2w x - - _ -> - let n = width compact - in fromDescList' (n+1) (x : toList compact) +cons !x (Partition !vec) + | V.null vec = Partition (if x > 0 then V.singleton y else V.empty) + | y >= h = Partition (V.cons y vec) + | otherwise = error "Partition/cons: invalid element to cons" + where + y = i2w x + h = V.head vec -------------------------------------------------------------------------------- -- | We assume that the element is not bigger than the last element! snoc :: Partition -> Int -> Partition -snoc !compact 0 = compact -snoc !compact !x = case compact of - - Nibble 0 -> singleton x - - Nibble word - | n < 15 -> Nibble $ (word + 1) .|. shiftL (i2w x) ((15-n)*4) - | otherwise -> makeMedium (n+1) (toList compact ++ [x]) - where - n = w2i (word .&. 15) - - Medium1 w1 - | n < 7 -> Medium1 $ (w1 + 1) .|. shiftL (i2w x) ((7-n)*8) - | otherwise -> Medium2 ((w1 .&. 0xffffffffffffff00) + i2w x) 8 - where - n = w2i (w1 .&. 255) - - Medium2 w1 w2 - | n < 15 -> Medium2 w1 $ (w2 + 1) .|. shiftL (i2w x) ((15-n)*8) - | otherwise -> Medium3 w1 ((w2 .&. 0xffffffffffffff00) + i2w x) 16 - where - n = w2i (w2 .&. 255) - - Medium3 w1 w2 w3 - | n < 23 -> Medium3 w1 w2 $ (w3 + 1) .|. shiftL (i2w x) ((23-n)*8) - | otherwise -> Medium4 w1 w2 ((w3 .&. 0xffffffffffffff00) + i2w x) 24 - where - n = w2i (w3 .&. 255) - - Medium4 w1 w2 w3 w4 - | n < 31 -> Medium4 w1 w2 w3 $ (w4 + 1) .|. shiftL (i2w x) ((31-n)*8) - | otherwise -> makeWordList (n + 1) (toList compact ++ [x]) - where - n = w2i (w4 .&. 255) - - WordList n list -> WordList (n+1) (go list) where - go :: [Word64] -> [Word64] - go (w:[]) = case mod n 4 of - 0 -> w : shiftL (i2w x) 48 : [] - 1 -> w + shiftL (i2w x) 32 : [] - 2 -> w + shiftL (i2w x) 16 : [] - 3 -> w + (i2w x) : [] - go (w:ws) = w : go ws - go [] = shiftL (i2w x) 48 : [] - -{- - _ -> - let n = width compact - in makeWordList (n+1) (toList compact ++ [x]) --} +snoc (Partition !vec) !x + | x == 0 = Partition vec + | V.null vec = Partition (V.singleton y) + | y <= V.last vec = Partition (V.snoc vec y) + | otherwise = error "Partition/snoc: invalid element to snoc" + where + y = i2w x -------------------------------------------------------------------------------- -- * exponential form @@ -308,33 +149,15 @@ -- | Width, or the number of parts width :: Partition -> Int -width compact = case compact of - Nibble word -> w2i (word .&. 15) - Medium1 word -> w2i (word .&. 255) - Medium2 _ word -> w2i (word .&. 255) - Medium3 _ _ word -> w2i (word .&. 255) - Medium4 _ _ _ word -> w2i (word .&. 255) - WordList n _ -> n +width (Partition vec) = V.vecLen vec -- | Height, or the first (that is, the largest) element height :: Partition -> Int -height compact = case compact of - Nibble word -> w2i (shiftR word 60) - Medium1 word -> w2i (shiftR word 56) - Medium2 word _ -> w2i (shiftR word 56) - Medium3 word _ _ -> w2i (shiftR word 56) - Medium4 word _ _ _ -> w2i (shiftR word 56) - WordList _ (word:_) -> w2i (shiftR word 48) +height (Partition vec) = w2i (V.head vec) -- | Width and height widthHeight :: Partition -> (Int,Int) -widthHeight compact = case compact of - Nibble word -> ( w2i (word .&. 15) , w2i (shiftR word 60) ) - Medium1 word -> ( w2i (word .&. 255) , w2i (shiftR word 56) ) - Medium2 word1 word2 -> ( w2i (word2 .&. 255) , w2i (shiftR word1 56) ) - Medium3 word1 _ word3 -> ( w2i (word3 .&. 255) , w2i (shiftR word1 56) ) - Medium4 word1 _ _ word4 -> ( w2i (word4 .&. 255) , w2i (shiftR word1 56) ) - WordList n (word:_) -> ( n , w2i (shiftR word 48) ) +widthHeight (Partition vec) = (V.vecLen vec , w2i (V.head vec)) -------------------------------------------------------------------------------- -- * Differential sequence @@ -342,132 +165,35 @@ -- | From a non-increasing sequence @[a1,a2,..,an]@ this computes the sequence of differences -- @[a1-a2,a2-a3,...,an-0]@ diffSequence :: Partition -> [Int] -diffSequence compact = case compact of - - Nibble 0 -> [] - - Nibble w -> - let !nw = (w .&. 15) - !w' = w - nw - !n = w2i nw - in [ w2i $ (shiftR w (60 - i*4) - shiftR w' (56 - i*4)) .&. 15 | i<-[0..n-1] ] - - Medium1 w -> - let !nw = (w .&. 255) - !w' = w - nw - !n = w2i nw - in [ w2i $ (shiftR w (56 - i*8) - shiftR w' (48 - i*8)) .&. 255 | i<-[0..n-1] ] - - Medium2 w1 w2 -> - let !nw = (w2 .&. 255) - !w2' = w2 - nw - !n = w2i nw - in [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1 (48 - i*8)) .&. 255 | i<-[0..6] ] ++ - ( w2i $ ( w1 - shiftR w2 56 ) .&. 255 ) : - [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2' (48 - i*8)) .&. 255 | i<-[0..n-9] ] - - Medium3 w1 w2 w3 -> - let !nw = (w3 .&. 255) - !w3' = w3 - nw - !n = w2i nw - in [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1 (48 - i*8)) .&. 255 | i<-[0..6] ] ++ - ( w2i $ ( w1 - shiftR w2 56 ) .&. 255 ) : - [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2 (48 - i*8)) .&. 255 | i<-[0..6] ] ++ - ( w2i $ ( w2 - shiftR w3 56 ) .&. 255 ) : - [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3' (48 - i*8)) .&. 255 | i<-[0..n-17] ] - - Medium4 w1 w2 w3 w4 -> - let !nw = (w4 .&. 255) - !w4' = w4 - nw - !n = w2i nw - in [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1 (48 - i*8)) .&. 255 | i<-[0..6] ] ++ - ( w2i $ ( w1 - shiftR w2 56 ) .&. 255 ) : - [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2 (48 - i*8)) .&. 255 | i<-[0..6] ] ++ - ( w2i $ ( w2 - shiftR w3 56 ) .&. 255 ) : - [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3 (48 - i*8)) .&. 255 | i<-[0..6] ] ++ - ( w2i $ ( w3 - shiftR w4 56 ) .&. 255 ) : - [ w2i $ (shiftR w4 (56 - i*8) - shiftR w4' (48 - i*8)) .&. 255 | i<-[0..n-25] ] - - WordList {} -> go (toList compact) where - go (x:ys@(y:_)) = (x-y) : go ys - go [x] = [x] - go [] = [] +diffSequence = go . toDescList where + go (x:ys@(y:_)) = (x-y) : go ys + go [x] = [x] + go [] = [] ---------------------------------------- -- | From a non-increasing sequence @[a1,a2,..,an]@ this computes the reversed sequence of differences -- @[ a[n]-0 , a[n-1]-a[n] , ... , a[2]-a[3] , a[1]-a[2] ] @ reverseDiffSequence :: Partition -> [Int] -reverseDiffSequence compact = case compact of - - Nibble 0 -> [] - - Nibble w -> - let !nw = (w .&. 15) - !w' = w - nw - !n = w2i nw - in [ w2i $ (shiftR w (60 - i*4) - shiftR w' (56 - i*4)) .&. 15 | i<-toZero (n-1) ] - - Medium1 w -> - let !nw = (w .&. 255) - !w' = w - nw - !n = w2i nw - in [ w2i $ (shiftR w (56 - i*8) - shiftR w' (48 - i*8)) .&. 255 | i<-toZero (n-1) ] - - Medium2 w1 w2 -> - let !nw = (w2 .&. 255) - !w2' = w2 - nw - !n = w2i nw - in [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2' (48 - i*8)) .&. 255 | i<-toZero (n-9) ] ++ - ( w2i $ ( w1 - shiftR w2 56 ) .&. 255 ) : - [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1 (48 - i*8)) .&. 255 | i<-toZero 6 ] - - Medium3 w1 w2 w3 -> - let !nw = (w3 .&. 255) - !w3' = w3 - nw - !n = w2i nw - in [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3' (48 - i*8)) .&. 255 | i<-toZero (n-17) ] ++ - ( w2i $ ( w2 - shiftR w3 56 ) .&. 255 ) : - [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2 (48 - i*8)) .&. 255 | i<-toZero 6 ] ++ - ( w2i $ ( w1 - shiftR w2 56 ) .&. 255 ) : - [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1 (48 - i*8)) .&. 255 | i<-toZero 6 ] - - Medium4 w1 w2 w3 w4 -> - let !nw = (w4 .&. 255) - !w4' = w4 - nw - !n = w2i nw - in [ w2i $ (shiftR w4 (56 - i*8) - shiftR w4' (48 - i*8)) .&. 255 | i<-toZero (n-25) ] ++ - ( w2i $ ( w3 - shiftR w4 56 ) .&. 255 ) : - [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3 (48 - i*8)) .&. 255 | i<-toZero 6 ] ++ - ( w2i $ ( w2 - shiftR w3 56 ) .&. 255 ) : - [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2 (48 - i*8)) .&. 255 | i<-toZero 6 ] ++ - ( w2i $ ( w1 - shiftR w2 56 ) .&. 255 ) : - [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1 (48 - i*8)) .&. 255 | i<-toZero 6 ] - - WordList {} -> (h : go asclist) where - asclist@(h:_) = toAscList compact - go (x:ys@(y:_)) = (y-x) : go ys - go [_] = [] - go [] = [] +reverseDiffSequence p = go (0 : toAscList p) where + go (x:ys@(y:_)) = (y-x) : go ys + go [x] = [] + go [] = [] -------------------------------------------------------------------------------- -- * Dual partition -foreign import ccall unsafe "c_dual_nibble" c_dual_nibble :: Word64 -> Word64 - dualPartition :: Partition -> Partition -dualPartition compact = case compact of - - Nibble 0 -> Nibble 0 - Nibble w -> Nibble (c_dual_nibble w) - _ -> if (w <= 255 && h <= 31) - then makeMedium h dualList - else makeWordList h dualList +dualPartition compact@(Partition vec) + | V.null vec = Partition V.empty + | otherwise = Partition (V.fromList' shape $ map i2w dual) where - (w,h) = widthHeight compact - dualList = concat + height = V.head vec + len = V.vecLen vec + shape = Shape (w2i height) (V.bitsNeededFor $ i2w len) + dual = concat [ replicate d j - | (j,d) <- zip (toOne w) (reverseDiffSequence compact) + | (j,d) <- zip (descendToOne len) (reverseDiffSequence compact) ] -------------------------------------------------------------------------------- @@ -478,106 +204,11 @@ -- | returns a descending (non-increasing) list toDescList :: Partition -> [Int] -toDescList compact = case compact of - - Nibble 0 -> [] - - Nibble word -> - let !n = w2i (word .&. 15) - in [ w2i (shiftR word (60 - i*4) .&. 15 ) | i<-[0..n-1] ] - - Medium1 word1 -> - let !n = w2i (word1 .&. 255) - in [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..n-1] ] - - Medium2 word1 word2 -> - let !n = w2i (word2 .&. 255) - in [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..7] ] ++ - [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-[0..n-9] ] - - Medium3 word1 word2 word3 -> - let !n = w2i (word3 .&. 255) - in [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..7] ] ++ - [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-[0..7] ] ++ - [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-[0..n-17] ] - - Medium4 word1 word2 word3 word4 -> - let !n = w2i (word4 .&. 255) - in [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..7] ] ++ - [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-[0..7] ] ++ - [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-[0..7] ] ++ - [ w2i (shiftR word4 (56 - i*8) .&. 255) | i<-[0..n-25] ] - - WordList _ list -> go list where - go :: [Word64] -> [Int] - go !wlist = case wlist of - (!w):(!ws) -> case ws of - (_:_) -> w2i (shiftR w 48 ) : - w2i (shiftR w 32 .&. 65535) : - w2i (shiftR w 16 .&. 65535) : - w2i ( w .&. 65535) : go ws - [] -> takeWhile (/=0) (fromWord w) - [] -> [] - - fromWord :: Word64 -> [Int] - fromWord !word = - [ w2i (shiftR word 48 ) - , w2i (shiftR word 32 .&. 65535) - , w2i (shiftR word 16 .&. 65535) - , w2i ( word .&. 65535) - ] - ----------------------------------------- +toDescList (Partition vec) = map w2i (V.toList vec) -- | Returns a reversed (ascending; non-decreasing) list toAscList :: Partition -> [Int] -toAscList compact = case compact of - - Nibble 0 -> [] - - Nibble word -> - let !n = w2i (word .&. 15) - in [ w2i (shiftR word (60 - i*4) .&. 15 ) | i<-toZero (n-1) ] - - Medium1 word1 -> - let !n = w2i (word1 .&. 255) - in [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero (n-1) ] - - Medium2 word1 word2 -> - let !n = w2i (word2 .&. 255) - in [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-toZero (n-9) ] ++ - [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero 7 ] - - Medium3 word1 word2 word3 -> - let !n = w2i (word3 .&. 255) - in [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-toZero (n-17) ] ++ - [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-toZero 7 ] ++ - [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero 7 ] - - Medium4 word1 word2 word3 word4 -> - let !n = w2i (word4 .&. 255) - in [ w2i (shiftR word4 (56 - i*8) .&. 255) | i<-toZero (n-25) ] ++ - [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-toZero 7 ] ++ - [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-toZero 7 ] ++ - [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero 7 ] - - WordList _ list -> dropWhile (==0) $ go (reverse list) where - go :: [Word64] -> [Int] - go !wlist = case wlist of - (!w):ws -> w2i ( w .&. 65535) : - w2i (shiftR w 16 .&. 65535) : - w2i (shiftR w 32 .&. 65535) : - w2i (shiftR w 48 ) : go ws - [] -> [] - -{- - go :: [Word64] -> [Int] - go (w:[]) = fromWord w - go (w:ws) = fromWord w ++ go ws - go [] = [] - fromWord :: Word64 -> [Int] - fromWord word = [ w2i (shiftR word (48 - i*16) .&. 65535) | i<-toZero 3 ] --} +toAscList (Partition vec) = map w2i (V.toRevList vec) -------------------------------------------------------------------------------- -- * Conversion from list @@ -590,129 +221,33 @@ :: Int -- ^ length -> [Int] -- ^ the list -> Partition -fromDescList' !n !list = - case list of - [] -> empty - (h:_) | h <= 0 -> empty - | h <= 15 && n <= 15 -> makeNibble n list - | h > 65535 -> error "partitions with elements bigger than 65535 are not supported" - | h > 255 || n > 31 -> makeWordList n list - | otherwise -> makeMedium n list - -makeNibble :: Int -> [Int] -> Partition -makeNibble !n list = Nibble $ go (i2w n) 60 list where - go !acc !k (x:xs) = go (acc + shiftL (i2w x) k) (k-4) xs - go !acc _ [] = acc -{- -makeNibble :: Int -> [Int] -> Partition -makeNibble !n list = Nibble - $ sum' [ shiftL (i2w x) (60 - 4*i) | (i,x) <- zip [0..] list ] - + i2w n --} - -makeMedium :: Int -> [Int] -> Partition -makeMedium !n list - | n <= 7 = makeMedium1 n list - | n <= 15 = makeMedium2 n list - | n <= 23 = makeMedium3 n list - | n <= 31 = makeMedium4 n list - | otherwise = error "makeMedium: input list too big (should be smaller than 32)" - -makeMedium1 :: Int -> [Int] -> Partition -makeMedium1 !n list = Medium1 - $ sum' [ shiftL (fromIntegral x) (56 - 8*i) | (i,x) <- zip [0..] list ] - + fromIntegral n - -makeMedium2 :: Int -> [Int] -> Partition -makeMedium2 !n list = Medium2 word1 word2 where - (list1,list2) = splitAt 8 list - word1 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list1 ] - word2 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list2 ] - + fromIntegral n - -makeMedium3 :: Int -> [Int] -> Partition -makeMedium3 !n list = Medium3 word1 word2 word3 where - (list1,tmp ) = splitAt 8 list - (list2,list3) = splitAt 8 tmp - word1 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list1 ] - word2 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list2 ] - word3 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list3 ] - + i2w n - -makeMedium4 :: Int -> [Int] -> Partition -makeMedium4 !n list = Medium4 word1 word2 word3 word4 where - (list1,tmp1 ) = splitAt 8 list - (list2,tmp2 ) = splitAt 8 tmp1 - (list3,list4) = splitAt 8 tmp2 - word1 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list1 ] - word2 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list2 ] - word3 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list3 ] - word4 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list4 ] - + i2w n - -makeWordList :: Int -> [Int] -> Partition -makeWordList !n list = WordList n (go list) where - go :: [Int] -> [Word64] - go !xs = case xs of - (x:y:z:w:rest) -> makeWord x y z w : go rest - (x:y:z: [] ) -> makeWord x y z 0 : [] - (x:y: [] ) -> makeWord x y 0 0 : [] - (x: [] ) -> makeWord x 0 0 0 : [] - [] -> [] - makeWord !x !y !z !w = shiftL (i2w x) 48 - + shiftL (i2w y) 32 - + shiftL (i2w z) 16 - + (i2w w) -{- - go [] = [] - go xs = case splitAt 4 xs of - (this,rest) -> case rest of - [] -> makeWord (take 4 $ this ++ repeat 0) : [] - _ -> makeWord this : go rest - makeWord [x,y,z,w] = shiftL (i2w x) 48 - + shiftL (i2w y) 32 - + shiftL (i2w z) 16 - + (i2w w) --} +fromDescList' !len !list = Partition (V.fromList' (Shape len bits) $ map i2w list) where + bits = case list of + [] -> 4 + (x:xs) -> V.bitsNeededFor (i2w x) -------------------------------------------------------------------------------- -- * Partial orderings +-- @ |p `isSubPartitionOf` q@ isSubPartitionOf :: Partition -> Partition -> Bool -isSubPartitionOf p q = case (p,q) of - - (Nibble 0 , _ ) -> True - - (Nibble u , Nibble v) -> let !n = w2i (u .&. 15) - in and [ (shiftR u (60 - i*4) .&. 15) - <= (shiftR v (60 - i*4) .&. 15) - | i<-[0..n-1] - ] - - _ -> and $ zipWith (<=) (toList p) (toList q ++ repeat 0) +isSubPartitionOf p q = and $ zipWith (<=) (toList p) (toList q ++ repeat 0) +-- | @q `dominates` p@ dominates :: Partition -> Partition -> Bool -dominates q p = case (q,p) of - - (_ , Nibble 0 ) -> True - - (Nibble v , Nibble u ) -> go 60 0 0 where - n = u .&. 15 - klimit = w2i (4*(15-n)) - go !k !b !a = if k <= klimit - then True - else let !b' = b + (shiftR v k .&. 15) - !a' = a + (shiftR u k .&. 15) - in if b' < a' - then False - else go (k-4) b' a' - - _ -> and $ zipWith (>=) (sums $ toList q ++ repeat 0) (sums $ toList p) where - sums = tail . scanl' (+) 0 +dominates (Partition vec_q) (Partition vec_p) = and $ zipWith (>=) (sums (qs ++ repeat 0)) (sums ps) where + sums = tail . scanl' (+) 0 + ps = V.toList vec_p + qs = V.toList vec_q -------------------------------------------------------------------------------- -- * Pieri rule +-- | Expands to product @s[lambda]*h[k]@ as a sum of @s[mu]@-s. See +pieriRule :: Partition -> Int -> [Partition] +pieriRule = error "Partitions/Integer/Compact: pieriRule not implemented yet" + +{- -- | Expands to product @s[lambda]*h[1] = s[lambda]*e[1]@ as a sum of @s[mu]@-s. See pieriRuleSingleBox :: Partition -> [Partition] pieriRuleSingleBox !compact = case compact of @@ -780,37 +315,38 @@ go !k [c] = if c==0 then 0 else shiftL (i2w c) k + 1 go !k (c:cs) = shiftL (i2w c) k + go (k-4) cs go !k [] = error "Compact/pieriRule/encode: shouldn't happen" +-} -------------------------------------------------------------------------------- -- * local (internally used) utility functions {-# INLINE i2w #-} -i2w :: Int -> Word64 +i2w :: Int -> Word i2w = fromIntegral {-# INLINE w2i #-} -w2i :: Word64 -> Int +w2i :: Word -> Int w2i = fromIntegral {-# INLINE sum' #-} -sum' :: [Word64] -> Word64 +sum' :: [Word] -> Word sum' = foldl' (+) 0 {-# INLINE safeTail #-} safeTail :: [Int] -> [Int] safeTail xs = case xs of { [] -> [] ; _ -> tail xs } -{-# INLINE toZero #-} -toZero :: Int -> [Int] -toZero !n - | n > 0 = n : toZero (n-1) +{-# INLINE descendToZero #-} +descendToZero :: Int -> [Int] +descendToZero !n + | n > 0 = n : descendToZero (n-1) | n == 0 = [0] | n < 0 = [] -{-# INLINE toOne #-} -toOne :: Int -> [Int] -toOne !n - | n > 1 = n : toOne (n-1) +{-# INLINE descendToOne #-} +descendToOne :: Int -> [Int] +descendToOne !n + | n > 1 = n : descendToOne (n-1) | n == 1 = [1] | n < 1 = [] diff --git a/Math/Combinat/Partitions/Integer/Naive.hs b/Math/Combinat/Partitions/Integer/Naive.hs --- a/Math/Combinat/Partitions/Integer/Naive.hs +++ b/Math/Combinat/Partitions/Integer/Naive.hs @@ -41,6 +41,18 @@ --------------------------------------------------------------------------------- +toList :: Partition -> [Int] +toList (Partition xs) = xs + +fromList :: [Int] -> Partition +fromList = mkPartition where + mkPartition xs = Partition $ sortBy (reverseCompare) $ filter (>0) xs + +fromListUnsafe :: [Int] -> Partition +fromListUnsafe = Partition + +--------------------------------------------------------------------------------- + isEmptyPartition :: Partition -> Bool isEmptyPartition (Partition p) = null p diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs --- a/Math/Combinat/Permutations.hs +++ b/Math/Combinat/Permutations.hs @@ -13,6 +13,7 @@ ( -- * The Permutation type Permutation (..) , fromPermutation + , lookupPermutation , (!!!) , permutationArray , permutationUArray , uarrayToPermutationUnsafe @@ -20,6 +21,7 @@ , maybePermutation , toPermutation , toPermutationUnsafe + , toPermutationUnsafeN , permutationSize -- * Disjoint cycles , DisjointCycles (..) @@ -54,16 +56,16 @@ , bubbleSort2 , bubbleSort -- * Permutation groups - , identity - , inverse - , multiply - , multiplyMany - , multiplyMany' + , identityPermutation + , inversePermutation + , multiplyPermutation + , productOfPermutations + , productOfPermutations' -- * Action of the permutation group - , permute + , permuteArray , permuteList - , permuteLeft , permuteRight - , permuteLeftList , permuteRightList + , permuteArrayLeft , permuteArrayRight + , permuteListLeft , permuteListRight -- * Sorting , sortingPermutationAsc , sortingPermutationDesc @@ -108,6 +110,9 @@ import Data.Array.MArray import Data.Array.Unsafe +import Data.Vector.Compact.WordVec ( WordVec ) +import qualified Data.Vector.Compact.WordVec as V + import Math.Combinat.ASCII import Math.Combinat.Classes import Math.Combinat.Helper @@ -117,10 +122,58 @@ import System.Random -------------------------------------------------------------------------------- +-- WordVec helpers + +toUArray :: WordVec -> UArray Int Int +toUArray vec = listArray (1,n) (map fromIntegral $ V.toList vec) where n = V.vecLen vec + +fromUArray :: UArray Int Int -> WordVec +fromUArray arr = fromPermListN n (map fromIntegral $ elems arr) where + (1,n) = bounds arr + +-- | maximum = length +fromPermListN :: Int -> [Int] -> WordVec +fromPermListN n perm = V.fromList' shape (map fromIntegral perm) where + shape = V.Shape n bits + bits = V.bitsNeededFor (fromIntegral n :: Word) + +fromPermList :: [Int] -> WordVec +fromPermList perm = V.fromList (map fromIntegral perm) + +(.!) :: WordVec -> Int -> Int +(.!) vec idx = fromIntegral (V.unsafeIndex (idx-1) vec) + +_elems :: WordVec -> [Int] +_elems = map fromIntegral . V.toList + +_assocs :: WordVec -> [(Int,Int)] +_assocs vec = zip [1..] (_elems vec) + +_bound :: WordVec -> Int +_bound = V.vecLen + +{- +-- the old internal representation (UArray Int Int) + +_elems :: UArray Int Int -> [Int] +_elems = elems + +_assocs :: UArray Int Int -> [(Int,Int)] +_assocs = elems + +_bound :: UArray Int Int -> Int +_bound = snd . bounds +-} + + +toPermN :: Int -> [Int] -> Permutation +toPermN n xs = Permutation (fromPermListN n xs) + +-------------------------------------------------------------------------------- -- * Types --- | A permutation. Internally it is an (unboxed) array of the integers @[1..n]@, with --- indexing range also being @(1,n)@. +-- | A permutation. Internally it is an (compact) vector +-- of the integers @[1..n]@. -- -- If this array of integers is @[p1,p2,...,pn]@, then in two-line -- notations, that represents the permutation @@ -136,12 +189,12 @@ -- -- (NOTE: this changed at version 0.2.8.0!) -- -newtype Permutation = Permutation (UArray Int Int) deriving (Eq,Ord) -- ,Show,Read) +newtype Permutation = Permutation WordVec deriving (Eq,Ord) -- ,Show,Read) instance Show Permutation where showsPrec d (Permutation arr) = showParen (d > 10) - $ showString "toPermutation " . showsPrec 11 (elems arr) -- app_prec = 10 + $ showString "toPermutation " . showsPrec 11 (_elems arr) -- app_prec = 10 instance Read Permutation where readsPrec d r = readParen (d > 10) fun r where @@ -164,25 +217,26 @@ newtype DisjointCycles = DisjointCycles [[Int]] deriving (Eq,Ord,Show,Read) fromPermutation :: Permutation -> [Int] -fromPermutation (Permutation ar) = elems ar +fromPermutation (Permutation ar) = _elems ar permutationUArray :: Permutation -> UArray Int Int -permutationUArray (Permutation ar) = ar +permutationUArray (Permutation ar) = toUArray ar --- | Note: this is slower than 'permutationUArray' permutationArray :: Permutation -> Array Int Int -permutationArray (Permutation ar) = listArray (1,n) (elems ar) where - (1,n) = bounds ar +permutationArray (Permutation ar) = listArray (1,n) (_elems ar) where + n = _bound ar -- | Assumes that the input is a permutation of the numbers @[1..n]@. toPermutationUnsafe :: [Int] -> Permutation -toPermutationUnsafe xs = Permutation perm where - n = length xs - perm = listArray (1,n) xs +toPermutationUnsafe xs = Permutation (fromPermList xs) +-- | This is faster than 'toPermutationUnsafe', but you need to supply @n@. +toPermutationUnsafeN :: Int -> [Int] -> Permutation +toPermutationUnsafeN n xs = Permutation (fromPermListN n xs) + -- | Note: Indexing starts from 1. uarrayToPermutationUnsafe :: UArray Int Int -> Permutation -uarrayToPermutationUnsafe = Permutation +uarrayToPermutationUnsafe = Permutation . fromUArray -- | Checks whether the input is a permutation of the numbers @[1..n]@. isPermutation :: [Int] -> Bool @@ -200,7 +254,7 @@ action :: forall s. ST s (Maybe Permutation) action = do ar <- newArray (1,n) 0 :: ST s (STUArray s Int Int) - let go [] = return $ Just (Permutation $ listArray (1,n) input) + let go [] = return $ Just (toPermutationUnsafe input) go (j:js) = if j<1 || j>n then return Nothing else do @@ -218,15 +272,27 @@ -- | Returns @n@, where the input is a permutation of the numbers @[1..n]@ permutationSize :: Permutation -> Int -permutationSize (Permutation ar) = snd $ bounds ar +permutationSize (Permutation ar) = _bound ar +-- | Returns the image @sigma(k)@ of @k@ under the permutation @sigma@. +-- +-- Note: we don't check the bounds! It may even crash if you index out of bounds! +lookupPermutation :: Permutation -> Int -> Int +lookupPermutation (Permutation ar) idx = ar .! idx + +-- infix 8 !!! + +-- | Infix version of 'lookupPermutation' +(!!!) :: Permutation -> Int -> Int +(!!!) (Permutation ar) idx = ar .! idx + instance HasWidth Permutation where width = permutationSize -- | Checks whether the permutation is the identity permutation isIdentityPermutation :: Permutation -> Bool -isIdentityPermutation (Permutation ar) = (elems ar == [1..n]) where - (1,n) = bounds ar +isIdentityPermutation (Permutation ar) = (_elems ar == [1..n]) where + n = _bound ar -- | Given a permutation of @n@ and a permutation of @m@, we return -- a permutation of @n+m@ resulting by putting them next to each other. @@ -255,18 +321,18 @@ -- | 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) = genericTwoLineNotation $ zip [1..] (elems arr) +twoLineNotation (Permutation arr) = genericTwoLineNotation $ zip [1..] (_elems arr) -- | 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)@. +-- as the bottom row of @twoLineNotation (inversePermutation perm)@. -- inverseTwoLineNotation :: Permutation -> ASCII inverseTwoLineNotation (Permutation arr) = - genericTwoLineNotation $ sortBy (comparing snd) $ zip [1..] (elems arr) + genericTwoLineNotation $ sortBy (comparing snd) $ zip [1..] (_elems arr) -- | Two-line notation for any set of numbers genericTwoLineNotation :: [(Int,Int)] -> ASCII @@ -301,7 +367,7 @@ numberOfCycles = numberOfCycles . permutationToDisjointCycles disjointCyclesToPermutation :: Int -> DisjointCycles -> Permutation -disjointCyclesToPermutation n (DisjointCycles cycles) = Permutation perm where +disjointCyclesToPermutation n (DisjointCycles cycles) = Permutation $ fromUArray perm where pairs :: [Int] -> [(Int,Int)] pairs xs@(x:_) = worker (xs++[x]) where @@ -329,7 +395,7 @@ permutationToDisjointCycles :: Permutation -> DisjointCycles permutationToDisjointCycles (Permutation perm) = res where - (1,n) = bounds perm + n = _bound perm -- we don't want trivial cycles f :: [Int] -> Bool @@ -357,7 +423,7 @@ worker :: STUArray s Int Bool -> Int -> Int -> [Int] -> ST s [Int] worker tag k l cyc = do writeArray tag l True - let m = perm ! l + let m = perm .! l if m == k then return cyc else worker tag k m (m:cyc) @@ -365,7 +431,7 @@ isEvenPermutation :: Permutation -> Bool isEvenPermutation (Permutation perm) = res where - (1,n) = bounds perm + n = _bound perm res = runST $ do tag <- newArray (1,n) False cycles <- unfoldM (step tag) 1 @@ -387,7 +453,7 @@ worker :: STUArray s Int Bool -> Int -> Int -> Int -> ST s Int worker tag k l cyclen = do writeArray tag l True - let m = perm ! l + let m = perm .! l if m == k then return cyclen else worker tag k m (1+cyclen) @@ -426,8 +492,8 @@ -- 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 ] + n = _bound arr + list = [ (i,j) | i<-[1..n-1], j<-[i+1..n], arr.!i > arr.!j ] -- | Returns the number of inversions: -- @@ -442,8 +508,8 @@ -- This should be @O(n*log(n))@ -- numberOfInversionsMerge :: Permutation -> Int -numberOfInversionsMerge (Permutation arr) = fst (sortCnt n $ elems arr) where - (_,n) = bounds arr +numberOfInversionsMerge (Permutation arr) = fst (sortCnt n $ _elems arr) where + n = _bound arr -- | First argument is length of the list. -- Returns also the inversion count. @@ -471,8 +537,8 @@ -- 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 ] + n = _bound 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: -- @@ -491,7 +557,7 @@ -- bubbleSort :: Permutation -> [Int] bubbleSort perm@(Permutation tgt) = runST action where - (_,n) = bounds tgt + n = _bound tgt action :: forall s. ST s [Int] action = do @@ -502,7 +568,7 @@ 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 + 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@ @@ -534,11 +600,11 @@ -- | The permutation @[n,n-1,n-2,...,2,1]@. Note that it is the inverse of itself. reversePermutation :: Int -> Permutation -reversePermutation n = Permutation $ listArray (1,n) [n,n-1..1] +reversePermutation n = Permutation $ fromPermListN n [n,n-1..1] -- | Checks whether the permutation is the reverse permutation @[n,n-1,n-2,...,2,1]. isReversePermutation :: Permutation -> Bool -isReversePermutation (Permutation arr) = elems arr == [n,n-1..1] where (1,n) = bounds arr +isReversePermutation (Permutation arr) = _elems arr == [n,n-1..1] where n = _bound arr -- | A transposition (swapping two elements). -- @@ -547,7 +613,7 @@ transposition :: Int -> (Int,Int) -> Permutation transposition n (i,j) = if i>=1 && j>=1 && i<=n && j<=n - then Permutation $ listArray (1,n) [ f k | k<-[1..n] ] + then Permutation $ fromPermListN n [ f k | k<-[1..n] ] else error "transposition: index out of range" where f k | k == i = j @@ -559,7 +625,7 @@ -- > transpositions n list == multiplyMany [ transposition n pair | pair <- list ] -- transpositions :: Int -> [(Int,Int)] -> Permutation -transpositions n list = Permutation (runSTUArray action) where +transpositions n list = Permutation (fromUArray $ runSTUArray action) where action :: ST s (STUArray s Int Int) action = do @@ -584,7 +650,7 @@ -- > adjacentTranspositions n list == multiplyMany [ adjacentTransposition n idx | idx <- list ] -- adjacentTranspositions :: Int -> [Int] -> Permutation -adjacentTranspositions n list = Permutation (runSTUArray action) where +adjacentTranspositions n list = Permutation (fromUArray $ runSTUArray action) where action :: ST s (STUArray s Int Int) action = do @@ -610,7 +676,7 @@ -- > ( 2 3 4 5 1 ) -- cycleLeft :: Int -> Permutation -cycleLeft n = Permutation $ listArray (1,n) $ [2..n] ++ [1] +cycleLeft n = Permutation $ fromPermListN n ([2..n] ++ [1]) -- | The permutation which cycles a list right by one step: -- @@ -622,52 +688,52 @@ -- > ( 5 1 2 3 4 ) -- cycleRight :: Int -> Permutation -cycleRight n = Permutation $ listArray (1,n) $ n : [1..n-1] +cycleRight n = Permutation $ fromPermListN n (n : [1..n-1]) -------------------------------------------------------------------------------- -- * Permutation groups --- | Multiplies two permutations together: @p `multiply` q@ +-- | Multiplies two permutations together: @p `multiplyPermutation` q@ -- means the permutation when we first apply @p@, and then @q@ -- (that is, the natural action is the /right/ action) -- --- See also 'permute' for our conventions. +-- See also 'permuteArray' for our conventions. -- -multiply :: Permutation -> Permutation -> Permutation -multiply pi1@(Permutation perm1) pi2@(Permutation perm2) = +multiplyPermutation :: Permutation -> Permutation -> Permutation +multiplyPermutation pi1@(Permutation perm1) pi2@(Permutation perm2) = if (n==m) - then Permutation result - else error "multiply: permutations of different sets" + then Permutation $ fromUArray result + else error "multiplyPermutation: permutations of different sets" where - (_,n) = bounds perm1 - (_,m) = bounds perm2 - result = permute pi2 perm1 + n = _bound perm1 + m = _bound perm2 + result = permuteArray pi2 (toUArray perm1) -infixr 7 `multiply` +infixr 7 `multiplyPermutation` -- | The inverse permutation. -inverse :: Permutation -> Permutation -inverse (Permutation perm1) = Permutation result +inversePermutation :: Permutation -> Permutation +inversePermutation (Permutation perm1) = Permutation $ fromUArray result where - result = array (1,n) $ map swap $ assocs perm1 - (_,n) = bounds perm1 + result = array (1,n) $ map swap $ _assocs perm1 + n = _bound perm1 -- | The identity (or trivial) permutation. -identity :: Int -> Permutation -identity n = Permutation $ listArray (1,n) [1..n] +identityPermutation :: Int -> Permutation +identityPermutation n = Permutation $ fromPermListN n [1..n] -- | Multiply together a /non-empty/ list of permutations (the reason for requiring the list to -- be non-empty is that we don\'t know the size of the result). See also 'multiplyMany''. -multiplyMany :: [Permutation] -> Permutation -multiplyMany [] = error "multiplyMany: empty list, we don't know size of the result" -multiplyMany ps = foldl1' multiply ps +productOfPermutations :: [Permutation] -> Permutation +productOfPermutations [] = error "productOfPermutations: empty list, we don't know size of the result" +productOfPermutations ps = foldl1' multiplyPermutation ps -- | Multiply together a (possibly empty) list of permutations, all of which has size @n@ -multiplyMany' :: Int -> [Permutation] -> Permutation -multiplyMany' n [] = identity n -multiplyMany' n ps@(p:_) = if n == permutationSize p - then foldl1' multiply ps - else error "multiplyMany': incompatible permutation size(s)" +productOfPermutations' :: Int -> [Permutation] -> Permutation +productOfPermutations' n [] = identityPermutation n +productOfPermutations' n ps@(p:_) = if n == permutationSize p + then foldl1' multiplyPermutation ps + else error "productOfPermutations': incompatible permutation size(s)" -------------------------------------------------------------------------------- -- * Action of the permutation group @@ -682,75 +748,75 @@ -- We adopt the convention that permutations act /on the right/ -- (as in Knuth): -- --- > permute pi2 (permute pi1 set) == permute (pi1 `multiply` pi2) set +-- > permuteArray pi2 (permuteArray pi1 set) == permuteArray (pi1 `multiplyPermutation` pi2) set -- --- Synonym to 'permuteRight' +-- Synonym to 'permuteArrayRight' -- -{-# SPECIALIZE permute :: Permutation -> Array Int b -> Array Int b #-} -{-# SPECIALIZE permute :: Permutation -> UArray Int Int -> UArray Int Int #-} -permute :: IArray arr b => Permutation -> arr Int b -> arr Int b -permute = permuteRight +{-# SPECIALIZE permuteArray :: Permutation -> Array Int b -> Array Int b #-} +{-# SPECIALIZE permuteArray :: Permutation -> UArray Int Int -> UArray Int Int #-} +permuteArray :: IArray arr b => Permutation -> arr Int b -> arr Int b +permuteArray = permuteArrayRight -- | Right action on lists. Synonym to 'permuteListRight' -- permuteList :: Permutation -> [a] -> [a] -permuteList = permuteRightList +permuteList = permuteListRight -- | The right (standard) action of permutations on sets. -- --- > permuteRight pi2 (permuteRight pi1 set) == permuteRight (pi1 `multiply` pi2) set +-- > permuteArrayRight pi2 (permuteArrayRight pi1 set) == permuteArrayRight (pi1 `multiplyPermutation` pi2) set -- -- The second argument should be an array with bounds @(1,n)@. -- The function checks the array bounds. -- -{-# SPECIALIZE permuteRight :: Permutation -> Array Int b -> Array Int b #-} -{-# SPECIALIZE permuteRight :: Permutation -> UArray Int Int -> UArray Int Int #-} -permuteRight :: IArray arr b => Permutation -> arr Int b -> arr Int b -permuteRight pi@(Permutation perm) ar = +{-# SPECIALIZE permuteArrayRight :: Permutation -> Array Int b -> Array Int b #-} +{-# SPECIALIZE permuteArrayRight :: Permutation -> UArray Int Int -> UArray Int Int #-} +permuteArrayRight :: IArray arr b => Permutation -> arr Int b -> arr Int b +permuteArrayRight pi@(Permutation perm) ar = if (a==1) && (b==n) - then listArray (1,n) [ ar!(perm!i) | i <- [1..n] ] - else error "permuteRight: array bounds do not match" + then listArray (1,n) [ ar!(perm.!i) | i <- [1..n] ] + else error "permuteArrayRight: array bounds do not match" where - (_,n) = bounds perm + n = _bound perm (a,b) = bounds ar -- | The right (standard) action on a list. The list should be of length @n@. -- --- > fromPermutation perm == permuteRightList perm [1..n] +-- > fromPermutation perm == permuteListRight perm [1..n] -- -permuteRightList :: forall a . Permutation -> [a] -> [a] -permuteRightList perm xs = elems $ permuteRight perm $ arr where +permuteListRight :: forall a . Permutation -> [a] -> [a] +permuteListRight perm xs = elems $ permuteArrayRight perm $ arr where arr = listArray (1,n) xs :: Array Int a n = permutationSize perm -- | The left (opposite) action of the permutation group. -- --- > permuteLeft pi2 (permuteLeft pi1 set) == permuteLeft (pi2 `multiply` pi1) set +-- > permuteArrayLeft pi2 (permuteArrayLeft pi1 set) == permuteArrayLeft (pi2 `multiplyPermutation` pi1) set -- --- It is related to 'permuteLeft' via: +-- It is related to 'permuteLeftArray' via: -- --- > permuteLeft pi arr == permuteRight (inverse pi) arr --- > permuteRight pi arr == permuteLeft (inverse pi) arr +-- > permuteArrayLeft pi arr == permuteArrayRight (inversePermutation pi) arr +-- > permuteArrayRight pi arr == permuteArrayLeft (inversePermutation pi) arr -- -{-# SPECIALIZE permuteLeft :: Permutation -> Array Int b -> Array Int b #-} -{-# SPECIALIZE permuteLeft :: Permutation -> UArray Int Int -> UArray Int Int #-} -permuteLeft :: IArray arr b => Permutation -> arr Int b -> arr Int b -permuteLeft pi@(Permutation perm) ar = +{-# SPECIALIZE permuteArrayLeft :: Permutation -> Array Int b -> Array Int b #-} +{-# SPECIALIZE permuteArrayLeft :: Permutation -> UArray Int Int -> UArray Int Int #-} +permuteArrayLeft :: IArray arr b => Permutation -> arr Int b -> arr Int b +permuteArrayLeft pi@(Permutation perm) ar = -- permuteRight (inverse pi) ar if (a==1) && (b==n) - then array (1,n) [ ( perm!i , ar!i ) | i <- [1..n] ] - else error "permuteLeft: array bounds do not match" + then array (1,n) [ ( perm.!i , ar!i ) | i <- [1..n] ] + else error "permuteArrayLeft: array bounds do not match" where - (_,n) = bounds perm + n = _bound perm (a,b) = bounds ar -- | The left (opposite) action on a list. The list should be of length @n@. -- --- > permuteLeftList perm set == permuteList (inverse perm) set --- > fromPermutation (inverse perm) == permuteLeftList perm [1..n] +-- > permuteListLeft perm set == permuteList (inversePermutation perm) set +-- > fromPermutation (inversePermutation perm) == permuteListLeft perm [1..n] -- -permuteLeftList :: forall a. Permutation -> [a] -> [a] -permuteLeftList perm xs = elems $ permuteLeft perm $ arr where +permuteListLeft :: forall a. Permutation -> [a] -> [a] +permuteListLeft perm xs = elems $ permuteArrayLeft perm $ arr where arr = listArray (1,n) xs :: Array Int a n = permutationSize perm @@ -843,7 +909,7 @@ forM_ [1..n] $ \i -> writeArray ar i i rnd' <- worker n (if isSattolo then n-1 else n) rnd ar perm <- Data.Array.Unsafe.unsafeFreeze ar - return (Permutation perm, rnd') + return (Permutation (fromUArray perm), rnd') worker :: RandomGen g => Int -> Int -> g -> STUArray s Int Int -> ST s g worker n m rnd ar = if n==1 diff --git a/Math/Combinat/RootSystems.hs b/Math/Combinat/RootSystems.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/RootSystems.hs @@ -0,0 +1,319 @@ + +-- | Naive (very inefficient) algorithm to generate the irreducible (Dynkin) root systems +-- +-- Based on + +{-# LANGUAGE BangPatterns, FlexibleInstances, TypeSynonymInstances, FlexibleContexts #-} +module Math.Combinat.RootSystems where + +-------------------------------------------------------------------------------- + +import Control.Monad + +import Data.Array + +import Data.Set (Set) +import qualified Data.Set as Set + +import Data.List +import Data.Ord + +import Math.Combinat.Numbers.Primes +import Math.Combinat.Sets + +-------------------------------------------------------------------------------- +-- * Half-integers + +-- | The type of half-integers (internally represented by their double) +-- +-- TODO: refactor this into its own module +newtype HalfInt + = HalfInt Int + deriving (Eq,Ord) + +half :: HalfInt +half = HalfInt 1 + +divByTwo :: Int -> HalfInt +divByTwo n = HalfInt n + +mulByTwo :: HalfInt -> Int +mulByTwo (HalfInt n) = n + +scaleBy :: Int -> HalfInt -> HalfInt +scaleBy k (HalfInt n) = HalfInt (k*n) + +instance Show HalfInt where + show (HalfInt n) = case divMod n 2 of + (k,0) -> show k + (_,1) -> show n ++ "/2" + +instance Num HalfInt where + fromInteger = HalfInt . (*2) . fromInteger + a + b = divByTwo $ mulByTwo a + mulByTwo b + a - b = divByTwo $ mulByTwo a - mulByTwo b + a * b = case divMod (mulByTwo a * mulByTwo b) 4 of + (k,0) -> HalfInt (2*k) + (k,2) -> HalfInt (2*k+1) + _ -> error "the result of multiplication is not a half-integer" + negate = divByTwo . negate . mulByTwo + signum = divByTwo . signum . mulByTwo + abs = divByTwo . abs . mulByTwo + +-------------------------------------------------------------------------------- +-- * Vectors of half-integers + +type HalfVec = [HalfInt] + +instance Num HalfVec where + fromInteger = error "HalfVec/fromInteger" + (+) = safeZip (+) + (-) = safeZip (-) + (*) = safeZip (*) + negate = map negate + abs = map abs + signum = map signum + +scaleVec :: Int -> HalfVec -> HalfVec +scaleVec k = map (scaleBy k) + +negateVec :: HalfVec -> HalfVec +negateVec = map negate + +-- dotProd :: HalfVec -> HalfVec +-- dotProd xs ys = foldl' (+) 0 $ safeZip (*) xs ys + +safeZip :: (a -> b -> c) -> [a] -> [b] -> [c] +safeZip f = go where + go (x:xs) (y:ys) = f x y : go xs ys + go [] [] = [] + go _ _ = error "safeZip: the lists do not have equal length" + +-------------------------------------------------------------------------------- +-- * Dynkin diagrams + +data Dynkin + = A !Int + | B !Int + | C !Int + | D !Int + | E6 | E7 | E8 + | F4 + | G2 + deriving (Eq,Show) + +-------------------------------------------------------------------------------- +-- * The roots of root systems + +-- | The ambient dimension of (our representation of the) system (length of the vector) +ambientDim :: Dynkin -> Int +ambientDim d = case d of + A n -> n+1 -- it's an n dimensional subspace of (n+1) dimensions + B n -> n + C n -> n + D n -> n + E6 -> 6 + E7 -> 8 -- sublattice of E8 ? + E8 -> 8 + F4 -> 4 + G2 -> 3 -- it's a 2 dimensional subspace of 3 dimensions + +simpleRootsOf :: Dynkin -> [HalfVec] +simpleRootsOf d = + + case d of + + A n -> [ e i - e (i+1) | i <- [1..n] ] + + B n -> [ e i - e (i+1) | i <- [1..n-1] ] ++ [e n] + + C n -> [ e i - e (i+1) | i <- [1..n-1] ] ++ [scaleVec 2 (e n)] + + D n -> [ e i - e (i+1) | i <- [1..n-1] ] ++ [e (n-1) + e n] + + E6 -> simpleRootsE6_123 + E7 -> simpleRootsE7_12 + E8 -> simpleRootsE8_even + + F4 -> [ [ 1,-1, 0, 0] + , [ 0, 1,-1, 0] + , [ 0, 0, 1, 0] + , [-h,-h,-h,-h] + ] + + G2 -> [ [ 1,-1, 0] + , [-1, 2,-1] + ] + + where + h = half + n = ambientDim d + + e :: Int -> HalfVec + e i = replicate (i-1) 0 ++ [1] ++ replicate (n-i) 0 + +positiveRootsOf :: Dynkin -> Set HalfVec +positiveRootsOf = positiveRoots . simpleRootsOf + +negativeRootsOf :: Dynkin -> Set HalfVec +negativeRootsOf = Set.map negate . positiveRootsOf + +allRootsOf :: Dynkin -> Set HalfVec +allRootsOf dynkin = Set.unions [ pos , neg ] where + simple = simpleRootsOf dynkin + pos = positiveRoots simple + neg = Set.map negate pos + +-------------------------------------------------------------------------------- +-- * Positive roots + +-- | Finds a vector, which is hopefully not orthognal to any root +-- (generated by the given simple roots), and has positive dot product with each of them. +findPositiveHyperplane :: [HalfVec] -> [Double] +findPositiveHyperplane vs = w where + n = length (head vs) + w0 = map (fromIntegral . mulByTwo) (foldl1 (+) vs) :: [Double] + w = zipWith (+) w0 perturb + perturb = map small $ map fromIntegral $ take n primes + small :: Double -> Double + small x = x / (10**10) + +positiveRoots :: [HalfVec] -> Set HalfVec +positiveRoots simples = Set.fromList pos where + roots = mirrorClosure simples + w = findPositiveHyperplane simples + pos = [ r | r <- Set.toList roots , dot4 r > 0 ] where + + dot4 :: HalfVec -> Double + dot4 a = foldl' (+) 0 $ safeZip (*) w $ map (fromIntegral . mulByTwo) a + +basisOfPositives :: Set HalfVec -> [HalfVec] +basisOfPositives set = Set.toList (Set.difference set set2) where + set2 = Set.fromList [ a + b | [a,b] <- choose 2 (Set.toList set) ] + + +-------------------------------------------------------------------------------- +-- * Operations on half-integer vectors + +-- | bracket b a = (a,b)/(a,a) +bracket :: HalfVec -> HalfVec -> HalfInt +bracket b a = + case divMod (2*a_dot_b) (a_dot_a) of + (n,0) -> divByTwo n + _ -> error "bracket: result is not a half-integer" + where + a_dot_b = foldl' (+) 0 $ safeZip (*) (map mulByTwo a) (map mulByTwo b) + a_dot_a = foldl' (+) 0 $ safeZip (*) (map mulByTwo a) (map mulByTwo a) + +-- | mirror b a = b - 2*(a,b)/(a,a) * a +mirror :: HalfVec -> HalfVec -> HalfVec +mirror b a = b - scaleVec (mulByTwo $ bracket b a) a + +-- | Cartan matrix of a list of (simple) roots +cartanMatrix :: [HalfVec] -> Array (Int,Int) Int +cartanMatrix list = array ((1,1),(n,n)) [ ((i,j), f i j) | i<-[1..n] , j<-[1..n] ] where + n = length list + arr = listArray (1,n) list + f !i !j = mulByTwo $ bracket (arr!j) (arr!i) + +printMatrix :: Show a => Array (Int,Int) a -> IO () +printMatrix arr = do + let ((1,1),(n,m)) = bounds arr + arr' = fmap show arr + let ks = [ 1 + maximum [ length (arr'!(i,j)) | i<-[1..n] ] | j<-[1..m] ] + forM_ [1..n] $ \i -> do + putStrLn $ flip concatMap [1..m] $ \j -> extendTo (ks!!(j-1)) $ arr' ! (i,j) + where + extendTo n s = replicate (n-length s) ' ' ++ s + +-------------------------------------------------------------------------------- +-- * Mirroring + +-- | We mirror stuff until there is no more things happening +-- (very naive algorithm, but seems to work) +mirrorClosure :: [HalfVec] -> Set HalfVec +mirrorClosure = go . Set.fromList where + + go set + | n' > n = go set' + | n'' > n = go set'' + | otherwise = set + where + n = Set.size set + n' = Set.size set' + n'' = Set.size set'' + set' = mirrorStep set + set'' = Set.union set (Set.map negateVec set) + +mirrorStep :: Set HalfVec -> Set HalfVec +mirrorStep old = Set.union old new where + new = Set.fromList [ mirror b a | [a,b] <- choose 2 $ Set.toList old ] + +-------------------------------------------------------------------------------- +-- * E6, E7 and E8 + +-- | This is a basis of E6 as the subset of the even E8 root system +-- where the first three coordinates agree (they are consolidated +-- into the first coordinate here) +simpleRootsE6_123:: [HalfVec] +simpleRootsE6_123 = roots where + h = half + roots = + [ [-h,-h,-h,-h,-h,-h,-h,-h] + , [ h, h, h, h, h, h,-h,-h] + , [ 0, 0, 0, 0,-1, 0, 1, 0] + , [ 0, 0, 0, 0, 0, 0,-1, 1] + , [-h,-h,-h, h, h, h, h,-h] + , [ 0, 0, 0,-1, 1, 0, 0, 0] + ] + +-- | This is a basis of E8 as the subset of the even E8 root system +-- where the first two coordinates agree (they are consolidated +-- into the first coordinate here) +simpleRootsE7_12:: [HalfVec] +simpleRootsE7_12 = roots where + h = half + roots = + [ [-h,-h,-h,-h,-h,-h,-h,-h] + , [ h, h, h, h, h, h,-h,-h] + , [ h, h,-h,-h,-h,-h, h, h] + , [-h,-h, h, h,-h, h, h,-h] + , [ 0, 0, 0,-1, 1, 0, 0, 0] + , [ 0, 0,-1, 1, 0, 0, 0, 0] + , [ 0, 0, 0, 0, 0, 0,-1, 1] + ] + +-- | This is a basis of E7 as the subset of the even E8 root system +-- for which the sum of coordinates sum to zero +simpleRootsE7_diag :: [HalfVec] +simpleRootsE7_diag = roots where + roots = [ e i - e (i+1) | i <-[2..7] ] ++ [[h,h,h,h,-h,-h,-h,-h]] + h = half + n = 8 + + e :: Int -> HalfVec + e i = replicate (i-1) 0 ++ [1] ++ replicate (n-i) 0 + +simpleRootsE8_even :: [HalfVec] +simpleRootsE8_even = roots where + roots = [v1,v2,v3,v4,v5,v7,v8,v6] + + [v1,v2,v3,v4,v5,v6,v7,v8] = roots0 + roots0 = [ e i - e (i+1) | i <-[1..6] ] ++ [ e 6 + e 7 , replicate 8 (-h) ] + + h = half + n = 8 + + e :: Int -> HalfVec + e i = replicate (i-1) 0 ++ [1] ++ replicate (n-i) 0 + +simpleRootsE8_odd :: [HalfVec] +simpleRootsE8_odd = roots where + roots = [ e i - e (i+1) | i <-[1..7] ] ++ [[-h,-h,-h,-h,-h , h,h,h]] + h = half + n = 8 + + e :: Int -> HalfVec + e i = replicate (i-1) 0 ++ [1] ++ replicate (n-i) 0 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Tableaux/LittlewoodRichardson.hs b/Math/Combinat/Tableaux/LittlewoodRichardson.hs --- a/Math/Combinat/Tableaux/LittlewoodRichardson.hs +++ b/Math/Combinat/Tableaux/LittlewoodRichardson.hs @@ -212,8 +212,8 @@ lrScalar lambdaMu alphaBeta = _lrScalar (fromSkewPartition lambdaMu) (fromSkewPartition alphaBeta) _lrScalar :: (Partition,Partition) -> (Partition,Partition) -> Int -_lrScalar (plam @(Partition lam ) , pmu @(Partition mu0) ) - (palpha@(Partition alpha) , pbeta@(Partition beta)) = +_lrScalar ( plam@( Partition lam ) , pmu@( Partition mu0 ) ) + ( palpha@(Partition alpha) , pbeta@(Partition beta) ) = if not (pmu `isSubPartitionOf` plam ) || not (pbeta `isSubPartitionOf` palpha) || (sum' lam + sum' beta) /= (sum' alpha + sum' mu0) -- equivalent to (lambda-mu) /= (alpha-beta) diff --git a/cbits/c_compact_partition.c b/cbits/c_compact_partition.c deleted file mode 100644 --- a/cbits/c_compact_partition.c +++ /dev/null @@ -1,24 +0,0 @@ - -#include - -// ----------------------------------------------------------------------------- - -uint64_t c_dual_nibble(uint64_t word) -{ - uint64_t n = (word & 15); // length - uint64_t h = (word >> 60); // height - uint64_t dual = h; // length of dual = height of original - uint64_t w = word - n; // zero out the low nibble - - uint64_t o = 60; - for(uint64_t i=0; i> (64-4*(n-i))) - - (w >> (60-4*(n-i))) ) & 15 ; // diff - for(uint64_t j=0;j= 4 && < 5, array >= 0.5, containers, random, transformers + Build-Depends: base >= 4 && < 5, + array >= 0.5, containers, random, transformers, + compact-word-vectors >= 0.2.0.2 Exposed-Modules: Math.Combinat Math.Combinat.Classes @@ -72,6 +74,7 @@ Math.Combinat.Trees.Nary Math.Combinat.Trees.Graphviz Math.Combinat.LatticePaths + Math.Combinat.RootSystems Math.Combinat.ASCII Math.Combinat.Helper Math.Combinat.TypeLevel @@ -85,8 +88,6 @@ Hs-Source-Dirs: . - C-Sources: cbits/c_compact_partition.c - ghc-options: -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports @@ -106,12 +107,14 @@ Tests.SkewTableaux Tests.Thompson Tests.Partitions.Integer - Tests.Partitions.Compact + Tests.Partitions.Compact Tests.Partitions.Skew Tests.Partitions.Ribbon - + Tests.Numbers.Primes + Tests.Numbers.Sequences + build-depends: base >= 4 && < 5, array >= 0.5, containers >= 0.5, random, transformers, - combinat, + combinat, compact-word-vectors >= 0.2.0.2, test-framework, test-framework-quickcheck2, QuickCheck >= 2, tasty, tasty-quickcheck, tasty-hunit diff --git a/test/TestSuite.hs b/test/TestSuite.hs --- a/test/TestSuite.hs +++ b/test/TestSuite.hs @@ -22,6 +22,17 @@ main :: IO () main = defaultMain tests +{- + +----- missing (because tasty, not test-framework): ----- + +Partitions.Compact +Numbers.Primes +Numbers.Sequences + +-} + + tests :: [Test] tests = [ testgroup_Permutations diff --git a/test/Tests/Braid.hs b/test/Tests/Braid.hs --- a/test/Tests/Braid.hs +++ b/test/Tests/Braid.hs @@ -226,7 +226,7 @@ prop_braidPerm_comp :: TwoBraids -> Bool prop_braidPerm_comp (TwoBraids b1 b2) = (p == q) where p = braidPermutation (compose b1 b2) - q = braidPermutation b1 `P.multiply` braidPermutation b2 + q = braidPermutation b1 `P.multiplyPermutation` braidPermutation b2 prop_link_positive :: PositiveBraid -> Bool prop_link_positive (PositiveB braid) = all (>=0) $ elems $ linkingMatrix braid diff --git a/test/Tests/Numbers/Primes.hs b/test/Tests/Numbers/Primes.hs new file mode 100644 --- /dev/null +++ b/test/Tests/Numbers/Primes.hs @@ -0,0 +1,106 @@ + + +-- | Tests for number theory +-- + +{-# LANGUAGE CPP, GeneralizedNewtypeDeriving, DataKinds, KindSignatures #-} +module Tests.Numbers.Primes where + +-------------------------------------------------------------------------------- + +{- +import Test.Framework +import Test.Framework.Providers.QuickCheck2 +import Test.QuickCheck +-} + +import Test.Tasty +import Test.Tasty.HUnit as U +import Test.Tasty.QuickCheck as Q + +import System.Random + +import Data.List +import Data.Ratio + +import GHC.TypeLits +import Data.Proxy + +import Math.Combinat.Sign +import Math.Combinat.Numbers +import Math.Combinat.Numbers.Primes +import Math.Combinat.Helper + +-------------------------------------------------------------------------------- + +prop_primes_sum_100 = sum (take 100 primes) @=? 24133 +prop_primes_sum_1000 = sum (take 1000 primes) @=? 3682913 +prop_primes_sum_10000 = sum (take 10000 primes) @=? 496165411 +prop_primes_sum_100000 = sum (take 100000 primes) @=? 62260698721 + +prop_divisorsigma1_sum_1000 = sum [ divisorSum n | n <- [1..1000] ] @=? 823081 +prop_divisorsigma1_sum_1000_b = sum [ divisorSum' 1 n | n <- [1..1000] ] @=? 823081 +prop_divisorsigma2_sum_1000 = sum [ divisorSum' 2 n | n <- [1..1000] ] @=? 401382971 +prop_divisorsigma3_sum_1000 = sum [ divisorSum' 3 n | n <- [1..1000] ] @=? 271161435595 + +prop_divisors_def n = sort (divisors n) == [ d | d<-[1..n] , d `divides` n ] + +prop_moebius_inversion n = sum [ moebiusMu d | d <- divisors n ] == (if n==1 then 1 else 0) + +prop_totient_divisorsum n = n == sum [ eulerTotient d | d <- divisors n ] + +prop_totient_mobius_inv n = eulerTotient n == sum [ moebiusMu d * div n d | d <- divisors n ] + +prop_Liouville_squaredivs n = liouvilleLambda n == rhs where + rhs = sum [ moebiusMu (div n d2) + | d <- divisors n , let d2 = d*d , d2 <= n , d2 `divides` n + ] + +prop_Liouville_sum n = sum [ liouvilleLambda d | d <- divisors n ] == (if isSquare n then 1 else 0) + +-------------------------------------------------------------------------------- + +prop_product_of_factors n = productOfFactors (factorize n) == n +prop_factorize_vs_naive n = factorize n == factorizeNaive n + +-------------------------------------------------------------------------------- + +main = defaultMain tests + +tests :: TestTree +tests = testGroup "Tests" [{-properties,-} unitTests] + +unitTests :: TestTree +unitTests = testGroup "Primes module" + [ unitTests1 + , unitTests2 + ] + +unitTests1 :: TestTree +unitTests1 = testGroup "Elementary number theory unit tests " + [ testCase "sum first 100 primes" $ prop_primes_sum_100 + , testCase "sum first 1000 primes" $ prop_primes_sum_1000 + , testCase "sum first 10000 primes" $ prop_primes_sum_10000 + , testCase "sum first 100000 primes" $ prop_primes_sum_100000 + , testCase "divisor set" $ allTrue [ prop_divisors_def n | n<-[1..1000] ] + , testCase "sum 1000 divisor sigma_1" $ prop_divisorsigma1_sum_1000 + , testCase "sum 1000 divisor sigma_1 /b" $ prop_divisorsigma1_sum_1000_b + , testCase "sum 1000 divisor sigma_2" $ prop_divisorsigma2_sum_1000 + , testCase "sum 1000 divisor sigma_3" $ prop_divisorsigma3_sum_1000 + , testCase "moebius inversion" $ allTrue [ prop_moebius_inversion n | n<-[1..1000] ] + , testCase "totient divisor sum" $ allTrue [ prop_totient_divisorsum n | n<-[1..1000] ] + , testCase "totient moebius inversion" $ allTrue [ prop_totient_mobius_inv n | n<-[1..1000] ] + , testCase "Liouville square divs sumn" $ allTrue [ prop_Liouville_squaredivs n | n<-[1..1000] ] + , testCase "Liouville divisor sum" $ allTrue [ prop_Liouville_sum n | n<-[1..1000] ] + ] + +unitTests2 :: TestTree +unitTests2 = testGroup "Integer factorization" + [ testCase "productOfFactors . factorize = id" $ allTrue [ prop_product_of_factors n | n<-[1..1000] ] + , testCase "factorize vs. factorizeNaive" $ allTrue [ prop_factorize_vs_naive n | n<-[1..1000] ] + ] + +allTrue :: [Bool] -> Assertion +allTrue bools = (and bools @=? True) + +-------------------------------------------------------------------------------- diff --git a/test/Tests/Numbers/Sequences.hs b/test/Tests/Numbers/Sequences.hs new file mode 100644 --- /dev/null +++ b/test/Tests/Numbers/Sequences.hs @@ -0,0 +1,66 @@ + + +-- | Tests for integer sequences +-- + +{-# LANGUAGE CPP, GeneralizedNewtypeDeriving, DataKinds, KindSignatures #-} +module Tests.Numbers.Sequences where + +-------------------------------------------------------------------------------- + +{- +import Test.Framework +import Test.Framework.Providers.QuickCheck2 +import Test.QuickCheck +-} + +import Test.Tasty +import Test.Tasty.HUnit as U +import Test.Tasty.QuickCheck as Q + +import System.Random + +import Data.List +import Data.Ratio + +import GHC.TypeLits +import Data.Proxy + +import Math.Combinat.Sign +import Math.Combinat.Numbers.Sequences +import Math.Combinat.Numbers.Primes +import Math.Combinat.Helper + +-------------------------------------------------------------------------------- + +prop_factorial_naive_vs_split n = factorialNaive n == factorialSplit n +prop_factorial_split_vs_swing n = factorialSplit n == factorialSwing n + +prop_fac_exponents_vs_naive n = factorialPrimeExponents n == factorialPrimeExponentsNaive n +prop_factorial_vs_fac_expos n = productOfFactors (factorialPrimeExponents n) == factorial n + +prop_double_factorial_naive_vs_split n = doubleFactorialNaive n == doubleFactorialSplit n + +prop_binomial_naive_vs_split n k = binomialNaive n k == binomialSplit n k + +-------------------------------------------------------------------------------- + +main = defaultMain tests + +tests :: TestTree +tests = testGroup "Tests" [{-properties,-} unitTests] + +unitTests :: TestTree +unitTests = testGroup "Numbers.Sequences module" + [ testCase "naive vs. split factorial" $ allTrue [ prop_factorial_naive_vs_split n | n<-[1..1000] ] + , testCase "split vs. swing factorial" $ allTrue [ prop_factorial_split_vs_swing n | n<-[1..1000] ] + , testCase "naive vs. fast factorial expos" $ allTrue [ prop_fac_exponents_vs_naive n | n<-[1..1000] ] + , testCase "factorial vs. fac exponents" $ allTrue [ prop_factorial_vs_fac_expos n | n<-[1..1000] ] + , testCase "naive vs. split double factorial" $ allTrue [ prop_double_factorial_naive_vs_split n | n<-[1..1000] ] + , testCase "naive vs. split binomial" $ allTrue [ prop_binomial_naive_vs_split n k | n<-[1..100] , k<-[-3..n+3] ] + ] + +allTrue :: [Bool] -> Assertion +allTrue bools = (and bools @=? True) + +-------------------------------------------------------------------------------- diff --git a/test/Tests/Partitions/Compact.hs b/test/Tests/Partitions/Compact.hs --- a/test/Tests/Partitions/Compact.hs +++ b/test/Tests/Partitions/Compact.hs @@ -3,16 +3,18 @@ -------------------------------------------------------------------------------- -import Data.List hiding ( uncons ) +import Data.List hiding ( uncons , singleton ) import Data.Ord import Test.Tasty import Test.Tasty.HUnit as U import Test.Tasty.QuickCheck as Q +import qualified Data.Vector.Compact.WordVec as Vec + import Math.Combinat.Partitions.Integer.Compact -import Math.Combinat.Partitions.Integer as P --- import qualified Math.Combinat.Parititions.Integer.IntList as P +-- import qualified Math.Combinat.Partitions.Integer as P +import qualified Math.Combinat.Partitions.Integer.IntList as P -------------------------------------------------------------------------------- @@ -23,6 +25,11 @@ -------------------------------------------------------------------------------- +cmp :: Partition -> Partition -> Ordering +cmp (Partition x) (Partition y) = Vec.cmpExtZero x y + +-------------------------------------------------------------------------------- + unitTests :: TestTree unitTests = testGroup "Unit tests" [ testCase "toList . fromList == id /1" $ allTrue [ xs == toList (fromDescList xs) | xs <- _testPartitions ] @@ -80,7 +87,7 @@ , testCase "dominates/15" $ allTrue [ P._dominates xs1 xs2 == dominates p1 p2 | xs1 <- _parts 15 , xs2 <- _parts 15 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ] , testCase "dominates/16" $ allTrue [ P._dominates xs1 xs2 == dominates p1 p2 | xs1 <- _parts 16 , xs2 <- _parts 16 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ] , testCase "dominates/17" $ allTrue [ P._dominates xs1 xs2 == dominates p1 p2 | xs1 <- _parts 17 , xs2 <- _parts 17 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ] - , testCase "pieriRuleSingleBox" $ allTrue [ pieriRuleSingleBox p =%%= map fromDescList (P._pieriRule xs 1) | xs <- _testPartitions , let p = fromDescList xs ] +-- , testCase "pieriRuleSingleBox" $ allTrue [ pieriRuleSingleBox p =%%= map fromDescList (P._pieriRule xs 1) | xs <- _testPartitions , let p = fromDescList xs ] , testCase "pieriRule" $ allTrue [ pieriRule p k =%%= map fromDescList (P._pieriRule xs k) | xs <- every10th _testPartitionsSmall , let p = fromDescList xs , k <- [1..2] ] ] @@ -122,7 +129,7 @@ , prop "fromExpo . toExpo == id" $ \p -> fromExponentialForm (toExponentialForm p) == p , prop "isSubPartitionOf" $ \p q -> isSubPartitionOf p q == P._isSubPartitionOf (toList p) (toList q) , prop "dominates" $ \p q -> dominates p q == P._dominates (toList p) (toList q) - , prop "pieriRuleSingleBox" $ \p -> map toList (pieriRuleSingleBox p) =%%= P._pieriRule (toList p) 1 +-- , prop "pieriRuleSingleBox" $ \p -> map toList (pieriRuleSingleBox p) =%%= P._pieriRule (toList p) 1 , localOption (QuickCheckTests 100) $ prop "pieriRule /1" $ \p (PieriK k) -> map toList (pieriRule p k) =%%= P._pieriRule (toList p) k diff --git a/test/Tests/Partitions/Integer.hs b/test/Tests/Partitions/Integer.hs --- a/test/Tests/Partitions/Integer.hs +++ b/test/Tests/Partitions/Integer.hs @@ -48,13 +48,20 @@ n <- choose (0, fromInteger (natVal (Proxy :: Proxy n))) myMkGen' Part (randomPartition n) +newtype Expo = Expo [Int] deriving (Eq,Show) + +instance Arbitrary Expo where + arbitrary = do + n <- choose (0, 10) + es <- replicateM n $ choose (0,4) + return $ Expo es + -------------------------------------------------------------------------------- -- * Types and instances newtype PartitionWeight = PartitionWeight Int deriving (Eq,Show) data PartitionWeightPair = PartitionWeightPair Int Int deriving (Eq,Show) data PartitionIntPair = PartitionIntPair Partition Int deriving (Eq,Show) - maxPartitionSize :: Int maxPartitionSize = 44 @@ -79,6 +86,16 @@ return (PartitionIntPair part k) -------------------------------------------------------------------------------- + +-- {- CONJUGATE LEXICOGRAPHIC ordering is a refinement of dominance partial ordering -} +-- let test n = [ ConjLex p >= ConjLex q | p <- partitions n , q <-partitions n , p `dominates` q ] +-- and (test 20) + +-- {- LEXICOGRAPHIC ordering is a refinement of dominance partial ordering -} +-- let test n = [ p >= q | p <- partitions n , q <-partitions n , p `dominates` q ] +-- and (test 20) + +-------------------------------------------------------------------------------- -- * test group testgroup_IntegerPartitions :: Test @@ -94,6 +111,11 @@ , testProperty "dominating partitions" prop_dominating_list , testProperty "counting partitions" prop_countParts , testProperty "union/sum duality" prop_union_sum_duality + -- + , testProperty "to/from expo vector" prop_to_from_expo_vector + , testProperty "to/from expo form" prop_to_from_expo_form + , testProperty "from/to expo vector" prop_from_to_expo_vector + , testProperty "from/to expo form" prop_from_to_expo_form ] -------------------------------------------------------------------------------- @@ -133,3 +155,12 @@ -------------------------------------------------------------------------------- +prop_to_from_expo_vector p = fromExponentVector (toExponentVector p) == p +prop_to_from_expo_form p = fromExponentialForm (toExponentialForm p) == p + +prop_from_to_expo_vector (Expo es) = toExponentVector (fromExponentVector es) == dropTailingZeros es + +prop_from_to_expo_form p = let ef = toExponentialForm p + in toExponentialForm (fromExponentialForm ef) == ef + +-------------------------------------------------------------------------------- diff --git a/test/Tests/Permutations.hs b/test/Tests/Permutations.hs --- a/test/Tests/Permutations.hs +++ b/test/Tests/Permutations.hs @@ -2,7 +2,7 @@ -- | Tests for permutations. -- -{-# LANGUAGE CPP, ScopedTypeVariables, GeneralizedNewtypeDeriving, FlexibleContexts #-} +{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, GeneralizedNewtypeDeriving, FlexibleContexts #-} module Tests.Permutations where -------------------------------------------------------------------------------- @@ -46,7 +46,7 @@ n = permutationSize perm permInternalSet :: Permutation -> Array Int Elem -permInternalSet perm@(Permutation arr) = listArray (1,n) [ Elem (arr!i) | i<-[1..n] ] where +permInternalSet perm@(Permutation arr) = listArray (1,n) [ Elem (perm !!! i) | i<-[1..n] ] where n = permutationSize perm sameSize :: Permutation -> Permutation -> Bool @@ -134,6 +134,7 @@ , testProperty "sign of permutation is multiplicative" prop_mulSign , testProperty "inverse is compatible with multiplication" prop_invMul + , testProperty "sign of permutation is parity of inversions" prop_sign_inversions , testProperty "parity of cyclic permutaiton" prop_cyclSign , testProperty "random permutation is valid" prop_permIsPerm @@ -178,38 +179,40 @@ prop_randCyclic cycl = ( isCyclicPermutation (fromCyclic cycl) ) -prop_inverse perm = ( perm == inverse (inverse perm) ) +prop_inverse perm = ( perm == inversePermutation (inversePermutation perm) ) prop_mulPerm (SameSize perm1 perm2) = - ( permute perm2 (permute perm1 set) == permute (perm1 `multiply` perm2) set ) + ( permuteArray perm2 (permuteArray perm1 set) == permuteArray (perm1 `multiplyPermutation` perm2) set ) where set = naturalSet perm1 prop_mulPermRight (SameSize perm1 perm2) = - ( permuteRight perm2 (permuteRight perm1 set) == permuteRight (perm1 `multiply` perm2) set ) + ( permuteArrayRight perm2 (permuteArrayRight perm1 set) == permuteArrayRight (perm1 `multiplyPermutation` perm2) set ) where set = naturalSet perm1 prop_mulPermLeft (SameSize perm1 perm2) = - ( permuteLeft perm2 (permuteLeft perm1 set) == permuteLeft (perm2 `multiply` perm1) set ) + ( permuteArrayLeft perm2 (permuteArrayLeft perm1 set) == permuteArrayLeft (perm2 `multiplyPermutation` 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_perm perm = permuteArray perm (naturalSet perm) == permInternalSet perm +prop_permLeft perm = permuteArrayLeft perm (permInternalSet perm) == naturalSet perm +prop_permRight perm = permuteArrayRight perm (naturalSet perm) == permInternalSet perm +prop_permLeftRight perm = permuteArrayLeft (inversePermutation perm) (naturalSet perm) == permuteArrayRight (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) ) + ( sgn perm1 * sgn perm2 == sgn (perm1 `multiplyPermutation` perm2) ) where sgn = signValue . signOfPermutation :: Permutation -> Int +prop_sign_inversions perm = signOfPermutation perm == paritySign (numberOfInversions perm) + prop_invMul (SameSize perm1 perm2) = - ( inverse perm2 `multiply` inverse perm1 == inverse (perm1 `multiply` perm2) ) + ( inversePermutation perm2 `multiplyPermutation` inversePermutation perm1 == inversePermutation (perm1 `multiplyPermutation` perm2) ) prop_cyclSign cycl = ( isEvenPermutation perm == odd n ) where perm = fromCyclic cycl @@ -221,17 +224,17 @@ isEvenAlternative p = even $ sum $ map (\x->x-1) $ map length $ fromDisjointCycles $ permutationToDisjointCycles p -prop_bubbleSort perm = multiplyMany' n (map (adjacentTransposition n) $ bubbleSort perm) == perm where +prop_bubbleSort perm = productOfPermutations' n (map (adjacentTransposition n) $ bubbleSort perm) == perm where n = permutationSize perm -prop_bubbleSort2 perm = multiplyMany' n (map (transposition n) $ bubbleSort2 perm) == perm where +prop_bubbleSort2 perm = productOfPermutations' n (map (transposition n) $ bubbleSort2 perm) == perm where n = permutationSize perm prop_bubble_inversions perm = length (bubbleSort perm) == numberOfInversions perm prop_number_inversions perm = length (inversions perm) == numberOfInversions perm -prop_ninversions_inverse perm = numberOfInversions perm == numberOfInversions (inverse perm) +prop_ninversions_inverse perm = numberOfInversions perm == numberOfInversions (inversePermutation perm) prop_merge_inversions perm = (numberOfInversionsMerge perm == numberOfInversionsNaive perm) diff --git a/test/Tests/Series.hs b/test/Tests/Series.hs --- a/test/Tests/Series.hs +++ b/test/Tests/Series.hs @@ -159,7 +159,7 @@ newtype NZRat = NZRat { fromNZRat :: Rational } deriving (Eq,Ord,Show,Num,Fractional) -- type parameter is for controlling the size (length), because some tests are too slow -newtype Ser (n :: Nat) = Ser { fromSer' :: [Integer] } deriving (Eq,Ord,Show) +newtype Ser (n :: Nat) = Ser { fromSer' :: [Integer] } deriving (Eq,Ord,Show) newtype SerR (n :: Nat) = SerR { fromSerR' :: [Rational] } deriving (Eq,Ord,Show) newtype Exp = Exp { fromExp :: Int } deriving (Eq,Ord,Show,Num,Random) @@ -371,7 +371,7 @@ prop_diff_integrate ser = (xs =!= differentiateSeries (integrateSeries xs)) where xs = fromSerR ser -prop_integrate_diff ser = (0 : tail xs =!= integrateSeries (differentiateSeries xs)) where +prop_integrate_diff ser = (null xs) || (0 : tail xs =!= integrateSeries (differentiateSeries xs)) where xs = fromSerR ser prop_cos_vs_cos2 = (cosSeries =!= (cosSeries2 :: [Rational]))