combinat 0.2.9.0 → 0.2.10.0

raw patch · 22 files changed

+1230/−767 lines, 22 filesdep +compact-word-vectorsdep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: compact-word-vectors

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• Math/Combinat/Groups/Braid.hs +9/−9
• Math/Combinat/Groups/Braid/NF.hs +9/−9
• Math/Combinat/Helper.hs +41/−1
• Math/Combinat/Numbers.hs +6/−3
• Math/Combinat/Numbers/Primes.hs +108/−6
• Math/Combinat/Numbers/Sequences.hs +118/−9
• Math/Combinat/Partitions/Integer.hs +74/−1
• Math/Combinat/Partitions/Integer/Compact.hs +99/−563
• Math/Combinat/Partitions/Integer/Naive.hs +12/−0
• Math/Combinat/Permutations.hs +172/−106
• Math/Combinat/RootSystems.hs +319/−0
• Math/Combinat/Tableaux/LittlewoodRichardson.hs +2/−2
• cbits/c_compact_partition.c +0/−24
• combinat.cabal +13/−10
• test/TestSuite.hs +11/−0
• test/Tests/Braid.hs +1/−1
• test/Tests/Numbers/Primes.hs +106/−0
• test/Tests/Numbers/Sequences.hs +66/−0
• test/Tests/Partitions/Compact.hs +12/−5
• test/Tests/Partitions/Integer.hs +32/−1
• test/Tests/Permutations.hs +18/−15
• test/Tests/Series.hs +2/−2

@@ -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 

@@ -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 && k<n  = (P.reversePermutation n `P.multiply` P.adjacentTransposition n k)+  | k>0 && k<n  = (P.reversePermutation n `P.multiplyPermutation` P.adjacentTransposition n k)   | otherwise   = error "posLPerm: index out of range"  -- | The permutation @posR k :: Braid n@ is realizing posRPerm :: Int -> Int -> Permutation posRPerm n k -  | k>0 && k<n  = (P.adjacentTransposition n k `P.multiply` P.reversePermutation n )+  | k>0 && k<n  = (P.adjacentTransposition n k `P.multiplyPermutation` P.reversePermutation n )   | otherwise   = error "posRPerm: index out of range"  --------------------------------------------------------------------------------@@ -333,15 +333,15 @@ -- > 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

@@ -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 

@@ -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+

@@ -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)]

@@ -4,39 +4,148 @@ -- See the \"On-Line Encyclopedia of Integer Sequences\", -- <https://oeis.org> . +{-# 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)

@@ -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<j +    then 0 : go (i+1) ef+    else e : go (i+1) rest+-}++-------------------------------------------------------------------------------- -- * Union and sum  -- | This is simply the union of parts. For example @@ -232,6 +275,14 @@ -------------------------------------------------------------------------------- -- * Dominating \/ dominated partitions +-- | Dominance partial ordering as a partial ordering.+dominanceCompare :: Partition -> 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

@@ -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 <https://en.wikipedia.org/wiki/Pieri's_formula>+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 <https://en.wikipedia.org/wiki/Pieri's_formula> 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  = [] 

@@ -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 

@@ -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 

@@ -0,0 +1,319 @@+
+-- | Naive (very inefficient) algorithm to generate the irreducible (Dynkin) root systems
+--
+-- Based on <https://en.wikipedia.org/wiki/Root_system>
+
+{-# 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 
+
+--------------------------------------------------------------------------------

@@ -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)

@@ -1,24 +0,0 @@--#include <stdint.h>--// -------------------------------------------------------------------------------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<n; i++)-  { uint64_t k = ( (w >> (64-4*(n-i))) -                 - (w >> (60-4*(n-i))) ) & 15 ;   // diff-    for(uint64_t j=0;j<k;j++) { dual |= (n-i) << o ; o -= 4 ; } -  }--  return dual;-}--// ------------------------------------------------------------------------------

@@ -1,5 +1,5 @@ Name:                combinat-Version:             0.2.9.0+Version:             0.2.10.0 Synopsis:            Generate and manipulate various combinatorial objects. Description:         A collection of functions to generate, count, manipulate                      and visualize all kinds of combinatorial objects like @@ -8,12 +8,12 @@ License:             BSD3 License-file:        LICENSE Author:              Balazs Komuves-Copyright:           (c) 2008-2018 Balazs Komuves+Copyright:           (c) 2008-2021 Balazs Komuves Maintainer:          bkomuves (plus) hackage (at) gmail (dot) com Homepage:            http://moire.be/haskell/ Stability:           Experimental Category:            Math-Tested-With:         GHC == 8.0.2+Tested-With:         GHC == 8.6.5 Cabal-Version:       1.24 Build-Type:          Simple @@ -24,13 +24,15 @@  source-repository head   type:                darcs -  location:            http://moire.be/haskell/projects/combinat/+  location:            https://hub.darcs.net/bkomuves/combinat  --------------------------------------------------------------------------------  Library -  Build-Depends:       base >= 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

@@ -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

@@ -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

@@ -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)++--------------------------------------------------------------------------------

@@ -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)++--------------------------------------------------------------------------------

@@ -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 

@@ -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++--------------------------------------------------------------------------------

@@ -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) 

@@ -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]))