packages feed

sym (empty) → 0.1

raw patch · 16 files changed

+2030/−0 lines, 16 filesdep +basedep +randomdep +vectorsetup-changed

Dependencies added: base, random, vector

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2012, Anders Claesson++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Anders Claesson nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Math/Sym.hs view
@@ -0,0 +1,342 @@+-- |+-- Module      : Math.Sym+-- Copyright   : (c) Anders Claesson 2012+-- License     : BSD-style+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>+-- +-- Provides an efficient definition of standard permutations,+-- 'StPerm', together with an abstract class, 'Perm', whose+-- functionality is largely inherited from 'StPerm' using a group+-- action and the standardization map.++module Math.Sym+    (+    -- * Standard permutations+      StPerm+    , toVector        -- :: StPerm -> Vector Int+    , fromVector      -- :: Vector Int -> StPerm+    , toList          -- :: StPerm -> [Int]+    , fromList        -- :: [Int] -> StPerm+    , (/-/)           -- :: StPerm -> StPerm -> StPerm+    , unrankStPerm    -- :: Int -> Integer -> StPerm+    , sym             -- :: Int -> [StPerm]++    -- * The permutation typeclass+    , Perm (..)++    -- * Generalize+    , generalize      -- :: Perm a => (StPerm -> StPerm) -> a -> a++    -- * Generating permutations+    , unrankPerm      -- :: Perm a => a -> Integer -> a+    , randomPerm      -- :: Perm a => a -> IO a+    , perms           -- :: Perm a => a -> [a]++    -- * Sorting operators+    , stackSort       -- :: Perm a => a -> a+    , bubbleSort      -- :: Perm a => a -> a++    -- * Permutation patterns+    , copies          -- :: Perm a => StPerm -> a -> [Set]+    , avoids          -- :: Perm a => [StPerm] -> a -> Bool+    , avoiders        -- :: Perm a => [StPerm] -> [a] -> [a]+    , av              -- :: [StPerm] -> Int -> [StPerm]++    -- * Subsets+    , Set+    , subsets         -- :: Int -> Int -> [Set]+    ) where++import Control.Monad (liftM)+import Data.Ord (comparing)+import Data.Monoid (Monoid(..))+import Data.Bits (Bits, bitSize, testBit, popCount, shiftL)+import Data.List (sort, sortBy)+import Data.Vector.Storable (Vector)+import qualified Data.Vector.Storable as SV (Vector, toList, fromList, fromListN, empty, map, (++))+import qualified Math.Sym.Internal as I+import Foreign.C.Types (CUInt(..))+++-- Standard permutations+-- ---------------------++-- | By a /standard permutation/ we shall mean a permutations of+-- @[0..k-1]@.+newtype StPerm = StPerm { perm0 :: I.Perm0 } deriving (Eq, Ord)++instance Show StPerm where+    show = show . toVector++instance Monoid StPerm where+    mempty = fromVector SV.empty+    mappend u v = fromVector $ (SV.++) u' v'+        where+          u' = toVector u+          v' = SV.map ( + size u) $ toVector v++-- | Convert a standard permutation to a vector.+toVector :: StPerm -> Vector Int+toVector = perm0++-- | Convert a vector to a standard permutation. The vector should a+-- permutation of the elements @[0..k-1]@ for some positive @k@. No+-- checks for this are done.+fromVector :: Vector Int -> StPerm+fromVector = StPerm++-- | Convert a standard permutation to a list.+toList :: StPerm -> [Int]+toList = SV.toList . toVector++-- | Convert a list to a standard permutation. The list should a+-- permutation of the elements @[0..k-1]@ for some positive @k@. No+-- checks for this are done.+fromList :: [Int] -> StPerm+fromList = fromVector . SV.fromList++infixl 6 /-/++-- | The /skew sum/ of two permutations. (A definition of the+-- /direct sum/ is provided by the Monoid instance.)+(/-/) :: StPerm -> StPerm -> StPerm+u /-/ v = fromVector $ (SV.++) u' v'+    where+      u' = SV.map ( + size v) $ toVector u+      v' = toVector v++-- | @unrankStPerm n rank@ is the @rank@-th (Myrvold & Ruskey)+-- permutation of @[0..n-1]@. E.g.,+-- +-- > unrankStPerm 16 19028390 == fromList [6,15,4,11,7,8,9,2,5,0,10,3,12,13,14,1]+-- +unrankStPerm :: Int -> Integer -> StPerm+unrankStPerm n = fromVector . I.unrankPerm n++-- | The list of standard permutations of the given size (the symmetric group). E.g.,+-- +-- > sym 2 == [fromList [0,1], fromList [1,0]]+-- +sym :: Int -> [StPerm]+sym n = map (unrankStPerm n) [0 .. product [1 .. toInteger n] - 1]+++-- The permutation typeclass+-- -------------------------++-- | The class of permutations. Minimal complete definition: 'st' and+-- 'act'. The default implementations of 'size' and 'idperm' can be+-- somewhat slow, so you may want to consider implementing them as+-- well.+class Perm a where++    -- | The standardization map. If there is an underlying linear+    -- order on @a@ then @st@ is determined by the unique order+    -- preserving map from @[0..]@ to that order. In any case, the+    -- standardization map should be equivariant with respect to the+    -- group action defined below; i.e., it should hold that+    -- +    -- > st (u `act` v) == u `act` st v+    -- +    st :: a -> StPerm++    -- | A (left) /group action/ of 'StPerm' on @a@. As for any group+    -- action it should hold that+    -- +    -- > (u `act` v) `act` w == u `act` (v `act` w)   &&   idperm u `act` v == v+    -- +    act :: StPerm -> a -> a++    -- | The size of a permutation. The default implementation derived from+    -- +    -- > size == size . st+    -- +    -- This is not a circular definition as 'size' on 'StPerm' is+    -- implemented independently. If the implementation of 'st' is+    -- slow, then it can be worth while to override the standard+    -- definiton; any implementation should, however, satisfy the+    -- identity above.+    {-# INLINE size #-}+    size :: a -> Int+    size = size . st++    -- | The identity permutation on the same underlying set as the+    -- given permutation. It should hold that+    -- +    -- > st (idperm u) == idperm (st u)+    -- +    -- Group theoretically, it should also hold that @u . inverse u ==+    -- idperm u@. In terms of the group action this means+    -- +    -- > idperm u == inverse (st u) `act` u+    -- +    -- and this is the default implementation.+    {-# INLINE idperm #-}+    idperm :: a -> a+    idperm u = inverse (st u) `act` u++    -- | The group theoretical inverse. It should hold that+    -- +    -- > inverse u == inverse (st u) `act` idperm u+    -- +    -- and this is the default implementation.+    {-# INLINE inverse #-}+    inverse :: a -> a+    inverse u = inverse (st u) `act` idperm u++    -- | Predicate determining if two permutations are+    -- order-isomorphic. The default implementation uses+    -- +    -- > u `ordiso` v  ==  u == st v+    -- +    -- Equivalently, one could use+    -- +    -- > u `ordiso` v  ==  inverse u `act` v == idperm v+    -- +    {-# INLINE ordiso #-}+    ordiso :: StPerm -> a -> Bool+    ordiso u v = u == st v++instance Perm StPerm where+    st         = id+    act u v    = fromVector $ I.act (toVector u) (toVector v)+    size       = I.size . toVector+    idperm     = fromVector . I.idperm . size+    inverse    = fromVector . I.inverse . toVector+    ordiso     = (==)++-- Auxiliary function: @w = act' u v@ iff @w[u[i]] = v[i]@.+-- Caveat: @act'@ is not a proper group action.+act' :: Ord a => [a] -> [b] -> [b]+act' u = map snd . sortBy (comparing fst) . zip u++instance (Enum a, Ord a) => Perm [a] where+    st         = fromVector . I.st . I.fromList . map fromEnum+    act u      = act' $ toList (inverse u)+    inverse v  = act' v (idperm v)+    size       = length+    idperm     = sort+++-- Generalize+-- ----------++-- | Generalize a function on 'StPerm' to a function on any permutations:+-- +-- > generalize f v = f (st v) `act` idperm v+-- +-- Note that this will only work as intended if @f@ is size preserving.+generalize :: Perm a => (StPerm -> StPerm) -> a -> a+generalize f v = f (st v) `act` idperm v+++-- Generating permutations+-- -----------------------++-- | @unrankPerm u rank@ is the @rank@-th (Myrvold & Ruskey)+-- permutation of @u@. E.g.,+-- +-- > unrankPerm ['1'..'9'] 88888 == "561297843"+-- +unrankPerm :: Perm a => a -> Integer -> a+unrankPerm u = (`act` u) . fromVector . I.unrankPerm (size u)++-- | @randomPerm u@ is a random permutation of @u@.+randomPerm :: Perm a => a -> IO a+randomPerm u = ((`act` u) . fromVector . I.fromLehmercode) `liftM` I.randomLehmercode (size u)++-- | All permutations of a given permutation. E.g.,+-- +-- > perms "123" == ["123","213","321","132","231","312"]+-- +perms :: Perm a => a -> [a]+perms u = map (`act` u) $ sym (size u)+++-- Sorting operators+-- -----------------++-- | One pass of stack-sort.+stackSort :: Perm a => a -> a+stackSort = generalize (fromVector . I.stackSort . toVector)++-- | One pass of bubble-sort.+bubbleSort :: Perm a => a -> a+bubbleSort = generalize (fromVector . I.bubbleSort . toVector)+++-- Permutation patterns+-- --------------------++-- | @copies p w@ is the list of (indices of) copies of the pattern+-- @p@ in the permutation @w@. E.g.,+-- +-- > copies (st "21") "2431" == [fromList [1,2],fromList [0,3],fromList [1,3],fromList [2,3]]+-- +copies :: Perm a => StPerm -> a -> [Set]+copies p w = I.copies subsets (toVector p) (toVector $ st w)++-- | @avoids ps w@ is a predicate determining if @w@ avoids the patterns @ps@.+avoids :: Perm a => [StPerm] -> a -> Bool+avoids ps w = all null [ copies p w | p <- ps ]++-- | @avoiders ps v@ is the list of permutations of @v@ avoiding the+-- patterns @ps@. This is equivalent to the definition+-- +-- > avoiders ps = filter (avoids ps)+-- +-- but is usually much faster.+avoiders :: Perm a => [StPerm] -> [a] -> [a]+avoiders ps = I.avoiders subsets (toVector . st) (map toVector ps)++-- | @av ps n@ is the list of permutations of @[0..n-1]@ avoiding the+-- patterns @ps@. E.g.,+-- +-- > map (length . av [st "132", st "321"]) [1..8] == [1,2,4,7,11,16,22,29]+-- +av :: [StPerm] -> Int -> [StPerm]+av ps = avoiders ps . sym+++-- Subsets+-- -------++-- | A set is represented by an increasing vector of non-negative+-- integers.+type Set = SV.Vector Int++-- A sub-class of 'Bits' used internally. Minimal complete definiton: 'next'.+class (Bits a, Integral a) => Bitmask a where+    -- | Lexicographically, the next bitmask with the same Hamming weight.+    next :: a -> a++    -- | @ones k m@ is the set of indices whose bits are set in+    -- @m@. Default implementation:+    -- +    -- > ones m = fromListN (popCount m) $ filter (testBit m) [0..]+    -- +    ones :: a -> Set+    ones m = SV.fromListN (popCount m) $ filter (testBit m) [0..]++instance Bitmask CUInt where+    next = I.nextCUInt+    ones = I.onesCUInt++instance Bitmask Integer where+    next = I.nextIntegral++-- @bitmasks n k@ is the list of bitmasks with Hamming weight @k@ and+-- size less than @2^n@.+bitmasks :: Bitmask a => Int -> Int -> [a]+bitmasks n k = take binomial (iterate next ((1 `shiftL` k) - 1))+    where+      n' = toInteger n+      k' = toInteger k+      binomial = fromIntegral $ product [n', n'-1 .. n'-k'+1] `div` product [1..k']++-- | @subsets n k@ is the list of subsets of @[0..n-1]@ with @k@+-- elements.+subsets :: Int -> Int -> [Set]+subsets n k = if n <= bitSize (0 :: CUInt)+              then map ones (bitmasks n k :: [CUInt])+              else map ones (bitmasks n k :: [Integer])
+ Math/Sym/D8.hs view
@@ -0,0 +1,118 @@+-- |+-- Module      : Math.Sym.D8+-- Copyright   : (c) Anders Claesson 2012+-- License     : BSD-style+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>+-- +-- The dihedral group of order 8 acting on permutations.+-- +-- To avoid name clashes this module is best imported @qualified@;+-- e.g.+-- +-- > import qualified Math.Sym.D8 as D8+-- ++module Math.Sym.D8+    (+    -- * The group elements+      r0            -- :: Perm a => a -> a+    , r1            -- :: Perm a => a -> a+    , r2            -- :: Perm a => a -> a+    , r3            -- :: Perm a => a -> a+    , s0            -- :: Perm a => a -> a+    , s1            -- :: Perm a => a -> a+    , s2            -- :: Perm a => a -> a+    , s3            -- :: Perm a => a -> a++    -- * D8, the klein four-group, and orbits+    , d8            -- :: Perm a => [a -> a]+    , klein4        -- :: Perm a => [a -> a]+    , orbit         -- :: Ord a => Perm a => [a -> a] -> a -> [a]++    -- * Aliases+    , id            -- :: Perm a => a -> a+    , rotate        -- :: Perm a => a -> a+    , complement    -- :: Perm a => a -> a+    , reverse       -- :: Perm a => a -> a+    , inverse       -- :: Perm a => a -> a+    ) where++import Prelude hiding (reverse, id)+import Data.List (group, sort)+import Math.Sym (Perm (size), fromVector, act)+import qualified Math.Sym (inverse)+import Math.Sym.Internal (revIdperm)++-- The group elements+-- ------------------++r0, r1, r2, r3, s0, s1, s2, s3 :: Perm a => a -> a++-- | Ration by 0 degrees, i.e. the identity map.+r0 w = w++-- | Ration by 90 degrees clockwise.+r1 = s2 . s1++-- | Ration by 2*90 = 180 degrees clockwise.+r2 = r1 . r1++-- | Ration by 3*90 = 270 degrees clockwise.+r3 = r2 . r1++-- | Reflection through a horizontal axis (also called 'complement').+s0 = r1 . s2++-- | Reflection through a vertical axis (also called 'reverse').+s1 w = (fromVector . revIdperm . size) w `act` w++-- | Reflection through the main diagonal (also called 'inverse').+s2 = Math.Sym.inverse++-- | Reflection through the anti-diagonal.+s3 = s1 . r1++-- D8, the klein four-group, and orbits+-- ------------------------------------++d8, klein4 :: Perm a => [a -> a]++-- | The dihedral group of order 8 (the symmetries of a square); that is,+-- +-- > d8 = [r0, r1, r2, r3, s0, s1, s2, s3]+-- +d8 = [r0, r1, r2, r3, s0, s1, s2, s3]++-- | The Klein four-group (the symmetries of a non-equilateral+-- rectangle); that is,+-- +-- > klein4 = [r0, r2, s0, s1]+-- +klein4 = [r0, r2, s0, s1]++-- | @orbit fs x@ is the orbit of @x@ under the functions in @fs@. E.g.,+-- +-- > orbit klein4 "2314" == ["1423","2314","3241","4132"]+-- +orbit :: Ord a => Perm a => [a -> a] -> a -> [a]+orbit fs x = map head . group $ sort [ f x | f <- fs ]++-- Aliases+-- -------++id, rotate, complement, reverse, inverse :: Perm a => a -> a++-- | @id = r0@+id = r0++-- | @rotate = r1@+rotate = r1++-- | @complement = s0@+complement = s0++-- | @reverse = s1@+reverse = s1++-- | @inverse = s2@+inverse = s2
+ Math/Sym/Internal.hs view
@@ -0,0 +1,481 @@+{-# LANGUAGE ForeignFunctionInterface #-}++-- |+-- Module      : Math.Sym.Internal+-- Copyright   : (c) Anders Claesson 2012+-- License     : BSD-style+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>+-- +-- An internal module used by the sym package.+-- +-- A Lehmercode is a vector of integers @w@ such @w!i <= length w - 1 - i@+-- for each @i@ in @[0..length w - 1]@; such a vector encodes a permutation.+-- This module implements /O(n)/ algorithms for unranking Lehmercodes and+-- permutations; the algorithms are due to W. Myrvold and F. Ruskey+-- [Ranking and Unranking Permutations in Linear Time, Information Processing+-- Letters, 79 (2001) 281-284].+-- +-- In addition, this module implements sorting operators, the+-- symmetries in D8 acting on permutations, as well as most of the+-- common permutation statistics.++module Math.Sym.Internal+    (+      Lehmercode+    , Perm0++    -- * Lehmercodes+    , unrankLehmercode+    , fromLehmercode+    , randomLehmercode+    , lehmercodes++    -- * Permutations+    , size+    , toList+    , fromList+    , act+    , unrankPerm+    , randomPerm+    , sym+    , idperm+    , revIdperm+    , sti+    , st+    , ordiso+    , copies+    , avoiders++    -- * Permutation symmetries+    , reverse+    , complement+    , inverse+    , rotate++    -- * Permutation statistics+    , asc     -- ascents+    , des     -- descents+    , exc     -- excedances+    , fp      -- fixed points+    , inv     -- inversions+    , maj     -- the major index+    , peak    -- peaks+    , vall    -- valleys+    , dasc    -- double ascents+    , ddes    -- double descents+    , lmin    -- left-to-right minima+    , lmax    -- left-to-right maxima+    , rmin    -- right-to-left minima+    , rmax    -- right-to-left maxima+    , head    -- the first element+    , last    -- the last element+    , lir     -- left-most increasing run+    , ldr     -- left-most decreasing run+    , rir     -- right-most increasing run+    , rdr     -- right-most decreasing run+    , comp    -- components+    , ep      -- rank a la Elizalde & Pak++    -- * Sorting operators+    , stackSort+    , bubbleSort++    -- * Bitmasks+    , onesCUInt+    , nextCUInt+    , nextIntegral+    ) where++import Prelude hiding (reverse, head, last)+import qualified Prelude (head)+import System.Random (getStdRandom, randomR)+import Control.Monad (forM_, liftM)+import Control.Monad.ST (runST)+import Data.List (group)+import Data.Bits (Bits, shiftR, (.|.), (.&.), popCount)+import qualified Data.Vector.Storable as SV+    ( Vector, toList, fromList, length, (!), thaw+    , unsafeFreeze, unsafeWith, enumFromN, enumFromStepN+    , head, last, filter, maximum, minimum, null, reverse, map+    )+import qualified Data.Vector.Storable.Mutable as MV+    ( unsafeNew, unsafeWrite, unsafeWith, swap, replicate+    )+import Foreign (Ptr, castPtr)+import System.IO.Unsafe (unsafePerformIO)+import Foreign.C.Types (CLong(..), CInt(..), CUInt(..))+import Foreign.Marshal.Utils (toBool)++-- | A Lehmercode is a vector of integers @w@ such @w!i <= length w - 1 - i@+-- for each @i@ in @[0..length w - 1]@.+type Lehmercode = SV.Vector Int++-- | By convention, a member of @Perm0@ is a permutation of some+-- finite subset of @[0..]@.+type Perm0 = SV.Vector Int+++-- Lehmercodes+-- -----------++-- | @unrankLehmercode n rank@ is the @rank@-th Lehmercode of length @n@.+unrankLehmercode :: Int -> Integer -> Lehmercode+unrankLehmercode n rank = runST $ do+  v <- MV.unsafeNew n+  iter v n rank (toInteger n)+  SV.unsafeFreeze v+    where+      {-# INLINE iter #-}+      iter _ 0 _ _ = return ()+      iter v i r m = do+        let (r',j) = quotRem r m+        MV.unsafeWrite v (n-i) (fromIntegral j)+        iter v (i-1) r' (m-1)++-- | Build a permutation from its Lehmercode.+fromLehmercode :: Lehmercode -> Perm0+fromLehmercode code = runST $ do+  let n = SV.length code+  v <- MV.unsafeNew n+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v i i+  forM_ [0..n-1] $ \i -> MV.swap v i (i + (SV.!) code i)+  SV.unsafeFreeze v++-- | A random Lehmercode of the given length.+randomLehmercode :: Int -> IO Lehmercode+randomLehmercode n = unrankLehmercode n `liftM` getStdRandom (randomR (0, factorial n - 1))++-- | The list of Lehmercodes of a given length.+lehmercodes :: Int -> [Lehmercode]+lehmercodes n = map (unrankLehmercode n) [0 .. factorial n - 1]+++-- Permutations+-- ------------++-- | The size of a permutation; the number of elements.+size :: Perm0 -> Int+size = SV.length++-- | The list of images of a permutation.+toList :: Perm0 -> [Int]+toList = SV.toList++-- | Make a permutation from a list of images.+fromList :: [Int] -> Perm0+fromList = SV.fromList++-- | @act u v@ is the permutation /w/ defined by /w(u(i)) = v(i)/.+act :: Perm0 -> Perm0 -> Perm0+act u v = runST $ do+  let n = SV.length u+  w <- MV.unsafeNew n+  forM_ [0..n-1] $ \i -> MV.unsafeWrite w i ((SV.!) v ((SV.!) u i))+  SV.unsafeFreeze w++factorial :: Integral a => a -> Integer+factorial = product . enumFromTo 1 . toInteger ++-- | @unrankPerm n rank@ is the @rank@-th (Myrvold & Ruskey) permutation of length @n@.+unrankPerm :: Int -> Integer -> Perm0+unrankPerm n = fromLehmercode . unrankLehmercode n++-- | A random permutation of the given length.+randomPerm :: Int -> IO Perm0+randomPerm n = fromLehmercode `liftM` randomLehmercode n++-- | @sym n@ is the list of permutations of @[0..n-1]@ (the symmetric group).+sym :: Int -> [Perm0]+sym n = map (unrankPerm n) [0 .. factorial n - 1]++-- | The identity permutation of the given length.+idperm :: Int -> Perm0+idperm = SV.enumFromN 0++-- | The reverse of the identity permutation.+revIdperm :: Int -> Perm0+revIdperm n = SV.enumFromStepN (n-1) (-1) n++-- | @sti w@ is the inverse of the standardization of @w@ (a+-- permutation on @[0..length w-1]@). E.g., @sti \<4,9,2\> ==+-- \<2,0,1\>@.+sti :: Perm0 -> Perm0+sti w = runST $ do+  let a = if SV.null w then 0 else SV.minimum w+  let b = if SV.null w then 0 else SV.maximum w+  let n = SV.length w+  v <- MV.replicate (1 + b - a) (-1)+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v ((SV.!) w i - a) i+  SV.filter (>=0) `liftM` SV.unsafeFreeze v++-- | The standardization map.+st :: Perm0 -> Perm0+st = inverse . sti++foreign import ccall unsafe "ordiso.h ordiso" c_ordiso+    :: Ptr CLong -> Ptr CLong -> Ptr CLong -> CLong -> CInt++-- | @ordiso u v m@ determines whether the subword in @v@ specified by+-- @m@ is order isomorphic to @u@.+ordiso :: Perm0 -> Perm0 -> SV.Vector Int -> Bool+ordiso u v m =+    let k = fromIntegral (SV.length u)+    in  unsafePerformIO $+        SV.unsafeWith u $ \u' ->+        SV.unsafeWith v $ \v' ->+        SV.unsafeWith m $ \m' ->+        return . toBool $ c_ordiso (castPtr u') (castPtr v') (castPtr m') k++-- | @copies subsets p w@ is the list of bitmasks that represent copies of @p@ in @w@.+copies :: (Int -> Int -> [SV.Vector Int]) -> Perm0 -> Perm0 -> [SV.Vector Int]+copies subsets p w = filter (ordiso p w) $ subsets n k+    where+      n = SV.length w+      k = SV.length p++avoiders1 :: (Int -> Int -> [SV.Vector Int]) -> (a -> Perm0) -> Perm0 -> [a] -> [a]+avoiders1 subsets f p ws =+    let ws0 = map f ws+        ws2 = zip ws0 ws+    in case group (map SV.length ws0) of+         []  -> []+         [_] -> let k = SV.length p+                    n = SV.length (Prelude.head ws0)+                in  [ v | (v0,v) <- ws2,  not $ any (ordiso p v0) (subsets n k) ]+         _   ->     [ v | (v0,v) <- ws2, null $ copies subsets p v0 ] ++-- | @avoiders subsets st ps ws@ is the list of permutations in @ws@+-- avoiding the patterns in @ps@.+avoiders :: (Int -> Int -> [SV.Vector Int]) -> (a -> Perm0) -> [Perm0] -> [a] -> [a]+avoiders _       _   []   ws = ws+avoiders subsets f (p:ps) ws = avoiders subsets f ps $ avoiders1 subsets f p ws+++-- Permutation symmetries+-- ----------------------++-- | @reverse \<a_1,...,a_n\> == \<a_n,,...,a_1\>@. E.g., @reverse+-- \<9,3,7,2\> == \<2,7,3,9\>@.+reverse :: Perm0 -> Perm0+reverse = SV.reverse++-- | @complement \<a_1,...,a_n\> == \<b_1,,...,b_n\>@, where @b_i = n - a_i - 1@.+-- E.g., @complement \<3,4,0,1,2\> == \<1,0,4,3,2\>@.+complement :: Perm0 -> Perm0+complement w = SV.map (\x -> SV.length w - x - 1) w++-- | @inverse w@ is the group theoretical inverse of @w@. E.g.,+-- @inverse \<1,2,0\> == \<2,0,1\>@.+inverse :: Perm0 -> Perm0+inverse w = runST $ do+  let n = SV.length w+  v <- MV.unsafeNew n+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v ((SV.!) w i) i+  SV.unsafeFreeze v++-- | The clockwise rotatation through 90 degrees. E.g.,+-- @rotate \<1,0,2\> == \<1,2,0\>@.+rotate :: Perm0 -> Perm0+rotate w = runST $ do+  let n = SV.length w+  v <- MV.unsafeNew n+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v ((SV.!) w (n-1-i)) i+  SV.unsafeFreeze v+++-- Permutation statistics+-- ----------------------++foreign import ccall unsafe "stat.h asc" c_asc+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h des" c_des+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h exc" c_exc+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h fp" c_fp+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h inv" c_inv+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h maj" c_maj+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h peak" c_peak+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h vall" c_vall+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h dasc" c_dasc+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h ddes" c_ddes+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h lmin" c_lmin+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h lmax" c_lmax+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h lir" c_lir+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h ldr" c_ldr+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h comp" c_comp+    :: Ptr CLong -> CLong -> CLong++foreign import ccall unsafe "stat.h ep" c_ep+    :: Ptr CLong -> CLong -> CLong++-- Marshal a permutation statistic defined in C to on in Haskell.+stat :: (Ptr CLong -> CLong -> CLong) -> Perm0 -> Int+stat f w = unsafePerformIO $+           SV.unsafeWith w $ \ptr ->+               return . fromIntegral $ f (castPtr ptr) (fromIntegral (SV.length w))++-- | First (left-most) value of a permutation.+head :: Perm0 -> Int+head = SV.head++-- | Last (right-most) value of a permutation.+last :: Perm0 -> Int+last = SV.last++-- | The number of right-to-left minima.+rmin :: Perm0 -> Int+rmin = lmin . SV.reverse++-- | The number of right-to-left maxima.+rmax :: Perm0 -> Int+rmax = lmax . SV.reverse++-- | The right-most increasing run.+rir :: Perm0 -> Int+rir = ldr . SV.reverse++-- | The right-most decreasing run.+rdr :: Perm0 -> Int+rdr = lir . SV.reverse++-- | The number of ascents.+asc :: Perm0 -> Int+asc = stat c_asc++-- | The number of descents.+des :: Perm0 -> Int+des = stat c_des++-- | The number of inversions.+inv :: Perm0 -> Int+inv = stat c_inv++-- | The major index.+maj :: Perm0 -> Int+maj = stat c_maj++-- | The number of peaks.+peak :: Perm0 -> Int+peak = stat c_peak++-- | The number of valleys.+vall :: Perm0 -> Int+vall = stat c_vall++-- | The number of double ascents.+dasc :: Perm0 -> Int+dasc = stat c_dasc++-- | The number of double descents.+ddes :: Perm0 -> Int+ddes = stat c_ddes++-- | The number of left-to-right minima.+lmin :: Perm0 -> Int+lmin = stat c_lmin++-- | The number of left-to-right maxima.+lmax :: Perm0 -> Int+lmax = stat c_lmax++-- | The left-most increasing run.+lir :: Perm0 -> Int+lir = stat c_lir++-- | The left-most decreasing run.+ldr :: Perm0 -> Int+ldr = stat c_ldr++-- | The number of excedances.+exc :: Perm0 -> Int+exc = stat c_exc++-- | The number of fixed points.+fp :: Perm0 -> Int+fp = stat c_fp++-- | The number of components.+comp :: Perm0 -> Int+comp = stat c_comp++-- | Rank as defined by Elizalde & Pak.+ep :: Perm0 -> Int+ep = stat c_ep+++-- Sorting operators+-- -----------------++foreign import ccall unsafe "sortop.h stacksort" c_stacksort+    :: Ptr CLong -> CLong -> IO ()++foreign import ccall unsafe "sortop.h bubblesort" c_bubblesort+    :: Ptr CLong -> CLong -> IO ()++-- Marshal a sorting operator defined in C to on in Haskell.+sortop :: (Ptr CLong -> CLong -> IO ()) -> Perm0 -> Perm0+sortop f w = unsafePerformIO $ do+               v <- SV.thaw w+               MV.unsafeWith v $ \ptr -> do+                 f (castPtr ptr) (fromIntegral (SV.length w))+                 SV.unsafeFreeze v++-- | One pass of stack-sort.+stackSort :: Perm0 -> Perm0+stackSort = sortop c_stacksort++-- | One pass of bubble-sort.+bubbleSort :: Perm0 -> Perm0+bubbleSort = sortop c_bubblesort+++-- Bitmasks+-- --------++foreign import ccall unsafe "bit.h next" c_next :: CUInt -> CUInt++-- | Lexicographically, the next 'CUInt' with the same Hamming weight.+nextCUInt :: CUInt -> CUInt+nextCUInt = c_next++foreign import ccall unsafe "bit.h ones" c_ones :: Ptr CUInt -> CUInt -> IO ()++-- | @onesCUInt k m@ gives the @k@ smallest indices whose bits are set in @m@.+onesCUInt :: CUInt -> SV.Vector Int+onesCUInt m = SV.map fromIntegral . unsafePerformIO $ do+                v <- MV.unsafeNew (popCount m)+                MV.unsafeWith v $ \ptr -> do+                  c_ones ptr m+                  SV.unsafeFreeze v++-- | Lexicographically, the next integral number with the same Hamming weight.+nextIntegral :: (Integral a, Bits a) => a -> a+nextIntegral a =+    let b = (a .|. (a - 1)) + 1+    in  b .|. ((((b .&. (-b)) `div` (a .&. (-a))) `shiftR` 1) - 1)
+ Math/Sym/Stat.hs view
@@ -0,0 +1,154 @@+-- |+-- Module      : Math.Sym.Stat+-- Copyright   : (c) Anders Claesson 2012+-- License     : BSD-style+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>+-- +-- Common permutation statistics. Please contact the maintainer if you+-- feel that there is a statistic that is missing; even better, send a+-- patch or make a pull request.+-- +-- To avoid name clashes this module is best imported @qualified@;+-- e.g.+-- +-- > import qualified Math.Sym.Stat as S+-- +-- For any permutation statistic @f@, below, it holds that @f w == f+-- (st w)@, and therefore the explanations below will be done on+-- standard permutations for convenience.++module Math.Sym.Stat +    (+      asc     -- ascents+    , des     -- descents+    , exc     -- excedances+    , fp      -- fixed points+    , inv     -- inversions+    , maj     -- the major index+    , peak    -- peaks+    , vall    -- valleys+    , dasc    -- double ascents+    , ddes    -- double descents+    , lmin    -- left-to-right minima+    , lmax    -- left-to-right maxima+    , rmin    -- right-to-left minima+    , rmax    -- right-to-left maxima+    , head    -- the first element+    , last    -- the last element+    , lir     -- left-most increasing run+    , ldr     -- left-most decreasing run+    , rir     -- right-most increasing run+    , rdr     -- right-most decreasing run+    , comp    -- components+    , ep      -- rank a la Elizalde & Pak+    ) where++import Prelude hiding (head, last)+import Math.Sym (Perm, toVector, st)+import Math.Sym.Internal (Perm0)+import qualified Math.Sym.Internal as I +    ( asc, des, exc, fp, inv, maj, peak, vall, dasc, ddes, lmin, lmax, rmin, rmax+    , head, last, lir, ldr, rir, rdr, comp, ep+    )++generalize :: Perm a => (Perm0 -> Int) -> a -> Int+generalize f = f . toVector . st++-- | The number of ascents. An /ascent/ in @w@ is an index @i@ such+-- that @w[i] \< w[i+1]@.+asc :: Perm a => a -> Int+asc = generalize I.asc++-- | The number of descents. A /descent/ in @w@ is an index @i@ such+-- that @w[i] > w[i+1]@.+des :: Perm a => a -> Int+des = generalize I.des++-- | The number of /excedances/: positions @i@ such that @w[i] > i@;+exc :: Perm a => a -> Int+exc = generalize I.exc++-- | The number of /fixed points/: positions @i@ such that @w[i] == i@;+fp :: Perm a => a -> Int+fp = generalize I.fp++-- | The number of /inversions/: pairs @\(i,j\)@ such that @i \< j@ and @w[i] > w[j]@+inv :: Perm a => a -> Int+inv = generalize I.inv++-- | /The major index/ is the sum of descents.+maj :: Perm a => a -> Int+maj = generalize I.maj++-- | The number of /peaks/: positions @i@ such that @w[i-1] \< w[i]@ and @w[i] \> w[i+1]@.+peak :: Perm a => a -> Int+peak = generalize I.peak++-- | The number of /valleys/: positions @i@ such that @w[i-1] \> w[i]@ and @w[i] \< w[i+1]@.+vall :: Perm a => a -> Int+vall = generalize I.vall++-- | The number of /double ascents/: positions @i@ such that @w[i-1] \<  w[i] \< w[i+1]@.+dasc :: Perm a => a -> Int+dasc = generalize I.dasc++-- | The number of /double descents/: positions @i@ such that @w[i-1] \>  w[i] \> w[i+1]@.+ddes :: Perm a => a -> Int+ddes = generalize I.ddes++-- | The number of /left-to-right minima/: positions @i@ such that @w[i] \< w[j]@ for all @j \< i@.+lmin :: Perm a => a -> Int+lmin = generalize I.lmin++-- | The number of /left-to-right maxima/: positions @i@ such that @w[i] \> w[j]@ for all @j \< i@.+lmax :: Perm a => a -> Int+lmax = generalize I.lmax++-- | The number of /right-to-left minima/: positions @i@ such that @w[i] \< w[j]@ for all @j \> i@.+rmin :: Perm a => a -> Int+rmin = generalize I.rmin++-- | The number of /right-to-left maxima/: positions @i@ such that @w[i] \> w[j]@ for all @j \> i@.+rmax :: Perm a => a -> Int+rmax = generalize I.rmax++-- | The first (left-most) element in the standardization. E.g., @head \"231\" = head (fromList [1,2,0]) = 1@.+head :: Perm a => a -> Int+head = generalize I.head++-- | The last (right-most) element in the standardization. E.g., @last \"231\" = last (fromList [1,2,0]) = 0@.+last :: Perm a => a -> Int+last = generalize I.last++-- | Length of the left-most increasing run: largest @i@ such that+-- @w[0] \< w[1] \< ... \< w[i-1]@.+lir :: Perm a => a -> Int+lir = generalize I.lir++-- | Length of the left-most decreasing run: largest @i@ such that+-- @w[0] \> w[1] \> ... \> w[i-1]@.+ldr :: Perm a => a -> Int+ldr = generalize I.ldr++-- | Length of the right-most increasing run: largest @i@ such that+-- @w[n-i] \< ... \< w[n-2] \< w[n-1]@.+rir :: Perm a => a -> Int+rir = generalize I.rir++-- | Length of the right-most decreasing run: largest @i@ such that+-- @w[n-i] \> ... \> w[n-2] \> w[n-1]@.+rdr :: Perm a => a -> Int+rdr = generalize I.rdr++-- | The number of components. E.g., @[2,0,3,1,4,6,7,5]@ has three+-- components: @[2,0,3,1]@, @[4]@ and @[6,7,5]@.+comp :: Perm a => a -> Int+comp = generalize I.comp++-- | The rank as defined by Elizalde and Pak [Bijections for+-- refined restricted permutations, /J. Comb. Theory, Ser. A/, 2004]:+-- +-- > maximum [ k | k <- [0..n-1], w[i] >= k for all i < k ]+-- +ep :: Perm a => a -> Int+ep = generalize I.ep
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ cbits/bit.c view
@@ -0,0 +1,19 @@+ #include <strings.h>++/* Lexicographically, the next bitmask with the same Hamming weight */+unsigned int+next(const unsigned int v)+{+	unsigned int t = v | (v - 1);+	return ((t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1)));+}++/* Positions of bits set */+void+ones(unsigned int *u, const unsigned int a)+{+	unsigned int b;++	for (b = a; b; b &= b-1)+		*u++ = ffs(b) - 1;+}
+ cbits/ordiso.c view
@@ -0,0 +1,27 @@+#include <stdlib.h>++/*+ * Determines whether the subword in v specified by m is order+ * isomorphic to u; len is the length of u.+ */+int+ordiso(const long *u, const long *v, const long *m, long len)+{+	register int i;+	long *w = malloc(len*sizeof(*w));+	long *w0 = w;++        /* Let w = v.m.u^{-1} */+	for (i = 0; i < len; i++, u++, m++)+		w[(int)*u] = v[(int)*m];++        /* Return 1 if w is increasing, 0 otherwise */+        for (; len > 1; len--, w++) {+                if (*w > *(w+1)) {+			free(w0);+			return 0;+		}+        }+	free(w0);+	return 1;+}
+ cbits/sortop.c view
@@ -0,0 +1,35 @@++/* One pass of stack-sort implemented a la Petter Br\"and\'en [Actions+ * on permutations and unimodality of descent polynomials, European+ * J. Combin. 29 (2008)]+ */+void+stacksort(long *w, long len) {+        int i = 0;+        int j = 0;+        int y;+        while (i < len) {+                j = i;+                y = w[j];+                while (y > w[j+1] && j+1 < len) {+                        w[j] = w[j+1];+                        j++;+                }+                w[j] = y;+                if (j == i)+			i++;+        }+}++/* On pass of bubble-sort */+void+bubblesort(long *w, long len) {+        int tmp;+	for (; len > 1; len--, w++) {+		if (*w > *(w+1)) {+			tmp    = *w;+			*w     = *(w+1);+			*(w+1) = tmp;+		}+	}+}
+ cbits/stat.c view
@@ -0,0 +1,244 @@+#include <string.h>++/* The number of ascents */+long+asc(const long *w, long len)+{+	long acc = 0;++	for (; len > 1; len--, w++) {+		if (*w < *(w+1))+			acc++;+	}+	return acc;+}+++/* The number of descents */+long+des(const long *w, long len)+{+	long acc = 0;++	for (; len > 1; len--, w++) {+		if (*w > *(w+1))+			acc++;+	}+	return acc;+}+++/* The number of excedances */+long+exc(const long *w, long len)+{+	long i, acc = 0;++	for (i = 0; i < len; i++, w++) {+		if (*w > i)+			acc++;+	}+	return acc;+}+++/* The number of fixed points */+long+fp(const long *w, long len)+{+	long i, acc = 0;++	for (i = 0; i < len; i++, w++) {+		if (*w == i)+			acc++;+	}+	return acc;+}+++/* The number of inversions */+long+inv(const long *w, long len)+{+	long i, j;+	long acc = 0;+	long *v;++	for (i = 0; i < len; i++, w++) {+		for (j = i+1, v = (long*)w+1; j < len; j++, v++) {+			if (*w > *v)+				acc++;+		}+	}+	return acc;+}+++/* The major index */+long+maj(const long *w, long len)+{+	long i, sum = 0;++	for (i = 1; i < len; i++, w++) {+		if (*w > *(w+1))+			sum += i;+	}+	return sum;+}+++/* The number of peaks */+long+peak(const long *w, long len)+{+	long acc = 0;++	for (; len > 2; len--, w++) {+		if (*w < *(w+1) && *(w+1) > *(w+2))+			acc++;+	}+	return acc;+}+++/* The number of valleys */+long+vall(const long *w, long len)+{+	long acc = 0;++	for (; len > 2; len--, w++) {+		if (*w > *(w+1) && *(w+1) < *(w+2))+			acc++;+	}+	return acc;+}+++/* The number of double ascents */+long+dasc(const long *w, long len)+{+	long acc = 0;++	for (; len > 2; len--, w++) {+		if (*w < *(w+1) && *(w+1) < *(w+2))+			acc++;+	}+	return acc;+}+++/* The number of double descents */+long+ddes(const long *w, long len)+{+	long acc = 0;++	for (; len > 2; len--, w++) {+		if (*w > *(w+1) && *(w+1) > *(w+2))+			acc++;+	}+	return acc;+}+++/* The number of left-to-right minima */+long+lmin(const long *w, long len)+{+	long m = *w + 1;+	long acc = 0;++	for (; len > 0; len--, w++) {+		if (*w < m) {+			m = *w;+			acc++;+		}+	}+	return acc;+}+++/* The number of left-to-right maxima */+long+lmax(const long *w, long len)+{+	long m = *w - 1;+	long acc = 0;++	for (; len > 0; len--, w++) {+		if (*w > m) {+			m = *w;+			acc++;+		}+	}+	return acc;+}++/* The left-most increasing run */+long+lir(const long *w, long len)+{+	long acc;++	if (len == 0)+		return 0;++	for (acc = 1; len > 1 && *w < *(w+1); len--, w++)+		acc++;++	return acc;+}++/* The left-most decreasing run */+long+ldr(const long *w, long len)+{+	long acc;++	if (len == 0)+		return 0;++	for (acc = 1; len > 1 && *w > *(w+1); len--, w++)+		acc++;++	return acc;+}++/* The number of components; O(len) */+long+comp(const long *w, long len)+{+	long i;+	long m = *w - 1;+	long acc = 0;++	for (i = 0; i < len; i++, w++) {+		if (*w > m)+			m = *w;+		if (m == i)+			acc++;+	}+	return acc;+}++/* rank as defined by Elizalde & Pak */+long+ep(const long *w, long len)+{+	long i;+	long m = *w;++	if (len == 0)+		return 0;++	for (i = 0; i < len; i++, w++) {+		if (*w <= m)+			m = *w;+		if (m <= i)+			return i;+	}+	return len;+}
+ include/bit.h view
@@ -0,0 +1,2 @@+unsigned int next(unsigned int);+void ones(unsigned int *, const unsigned int);
+ include/ordiso.h view
@@ -0,0 +1,1 @@+int ordiso(const long *, const long *, const long *, long);
+ include/sortop.h view
@@ -0,0 +1,1 @@+void stacksort(long *, long);
+ include/stat.h view
@@ -0,0 +1,16 @@+long asc  (const long *, long); /* ascents */+long des  (const long *, long); /* descents */+long exc  (const long *, long); /* excedances */+long fp   (const long *, long); /* fixed points */+long inv  (const long *, long); /* inversions */+long maj  (const long *, long); /* major index */+long peak (const long *, long); /* peaks */+long vall (const long *, long); /* valleys */+long dasc (const long *, long); /* double ascents */+long ddes (const long *, long); /* double descents */+long lmin (const long *, long); /* left-to-right minima */+long lmax (const long *, long); /* left-to-right maxima */+long lir  (const long *, long); /* left-most increasing run */+long ldr  (const long *, long); /* left-most decreasing run */+long comp (const long *, long); /* components */+long ep   (const long *, long); /* rank a la Elizalde & Pak */
+ sym.cabal view
@@ -0,0 +1,50 @@+Name:                sym+Version:             0.1+Synopsis:            Permutations, patterns, and statistics+Description:         +  Definitions for permutations with an emphasis on permutation+  patterns and statistics.+  .+  ["Math.Sym"] Provides an efficient definition of standard+  permutations (@StPerm@) together with a typeclass (@Perm@) whose+  functionality is largely inherited from @StPerm@ using a group+  action and the standardization map.+  .+  ["Math.Sym.D8"] The dihedral group of order 8 acting on permutations.+  .+  ["Math.Sym.Stat"] Common permutation statistics, such as @des@,+  @inv@, @exc@, @maj@, @fp@, @comp@, @lmin@, @lmax@, ...++Homepage:            http://github.com/akc/sym++License:             BSD3+License-file:        LICENSE+Author:              Anders Claesson+Maintainer:          anders.claesson@gmail.com+Category:            Math+Build-type:          Simple++Extra-source-files:  tests/Properties.hs++Cabal-version:       >=1.6++source-repository head+  type:                git+  location:            git://github.com/akc/sym.git++Library+  Exposed-modules:     Math.Sym+                       Math.Sym.D8+                       Math.Sym.Stat+                       Math.Sym.Internal++  Build-depends:       base >= 3 && < 5, random, vector+  +  ghc-prof-options:    -auto-all -caf-all+  ghc-options:         -Wall -O2+  cc-options:          -Wall++  c-sources:           cbits/stat.c, cbits/sortop.c, cbits/ordiso.c, cbits/bit.c+  include-dirs:        include+  includes:            stat.h, sortop.h, ordiso.h, bit.h+  install-includes:    stat.h, sortop.h, ordiso.h, bit.h
+ tests/Properties.hs view
@@ -0,0 +1,508 @@+-- |+-- Copyright   : (c) Anders Claesson 2012+-- License     : BSD-style+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>++import Data.List+import Data.Monoid+import Control.Monad+import qualified Math.Sym as Sym+import qualified Math.Sym.D8 as D8+import qualified Math.Sym.Stat as S+import qualified Math.Sym.Internal as I+import qualified Data.Vector.Storable as SV+import Test.QuickCheck++check :: Testable prop => prop -> IO ()+check = quickCheck++---------------------------------------------------------------------------------+-- Generators+---------------------------------------------------------------------------------++rank :: Int -> Gen Integer+rank n = choose (0, product [1..fromIntegral n] - 1)++lenRank :: Gen (Int, Integer)+lenRank = sized $ \m -> do+            n <- choose (0, m)+            r <- rank n+            return (n, r)++lenRank2 :: Gen (Int, Integer, Integer)+lenRank2 = do (n, r1) <- lenRank+              r2 <- rank n+              return (n, r1, r2)++moreThan :: Int -> Gen Int+moreThan x = (\d -> x + abs d) `liftM` choose (1, 100)++vecFrom :: Int -> Int -> Gen [Int]+vecFrom 0 _ = return []+vecFrom n x = moreThan x >>= liftM (x:) . vecFrom (n-1)++incVec :: Int -> Gen [Int]+incVec n = arbitrary >>= vecFrom n++-- The sub-permutation determined by a set of indices.+subperm :: Sym.Set -> Sym.StPerm -> Sym.StPerm+subperm m w = Sym.fromVector . I.st $ SV.map ((SV.!) (Sym.toVector w)) m++subperms :: Int -> Sym.StPerm -> [Sym.StPerm]+subperms k w = [ subperm m w | m <- Sym.subsets (Sym.size w) k ]++instance Arbitrary Sym.StPerm where+    arbitrary = uncurry Sym.unrankStPerm `liftM` lenRank+    shrink w = nub $ [0 .. Sym.size w - 1] >>= \k -> subperms k w++perm2 :: Gen (Sym.StPerm, [Int])+perm2 = do u <- arbitrary+           v <- incVec (Sym.size u)+           return (u, v)++perm3 :: Gen (Sym.StPerm, Sym.StPerm, [Int])+perm3 = do (n,r1,r2) <- lenRank2+           let u = Sym.unrankStPerm n r1+           let v = Sym.unrankStPerm n r2+           w <- incVec n+           return (u, v, w)++perm :: Gen [Int]+perm = liftM (uncurry Sym.act) perm2++newtype Symmetry = Symmetry (Sym.StPerm -> Sym.StPerm, String)++d8Symmetries :: [Symmetry]+d8Symmetries = [ Symmetry (D8.r0, "r0")+               , Symmetry (D8.r1, "r1")+               , Symmetry (D8.r2, "r2")+               , Symmetry (D8.r3, "r3")+               , Symmetry (D8.s0, "s0")+               , Symmetry (D8.s1, "s1")+               , Symmetry (D8.s2, "s2")+               , Symmetry (D8.s3, "s3")+               ]++instance Show Symmetry where+    show (Symmetry (_,s)) = s++instance Arbitrary Symmetry where+    arbitrary = liftM (d8Symmetries !!) $ choose (0, length d8Symmetries - 1)+++---------------------------------------------------------------------------------+-- Properties for Math.Sym+---------------------------------------------------------------------------------++prop_monoid_mempty1 w = mempty <> w == (w :: Sym.StPerm)+prop_monoid_mempty2 w = w <> mempty == (w :: Sym.StPerm)+prop_monoid_associative u v w = u <> (v <> w) == (u <> v) <> (w :: Sym.StPerm)++newtype S = S {unS :: Sym.StPerm} deriving (Eq, Show)++instance Arbitrary S where+    arbitrary = liftM S arbitrary++prop_monoid_mempty1_S w = mempty <> w == (w :: S)+prop_monoid_mempty2_S w = w <> mempty == (w :: S)+prop_monoid_associative_S u v w = u <> (v <> w) == (u <> v) <> (w :: S)++instance Monoid S where+    mempty = S $ Sym.fromVector SV.empty+    mappend u v = S $ (Sym./-/) (unS u) (unS v)++prop_unrankStPerm_distinct =+    forAll lenRank $ \(n, r) ->+        let w = Sym.toList (Sym.unrankStPerm n r) in nub w == w++prop_unrankStPerm_injective =+    forAll lenRank2 $ \(n, r1, r2) ->+        (Sym.unrankStPerm n r1 :: Sym.StPerm) /= Sym.unrankStPerm n r2 || r1 == r2++prop_sym = and [ sort (Sym.sym n) == sort (sym' n) | n<-[0..6] ]+    where+      sym' n = map Sym.fromList $ Data.List.permutations [0..fromIntegral n - 1]++prop_perm =+    and [ sort (Sym.perms [1..n]) == sort (permutations [1..n]) | n<-[0..6] ]++prop_st =+    forAll perm2 $ \(u,v) -> Sym.st (u `Sym.act` v) == u `Sym.act` Sym.st v++prop_act_def =+    forAll perm2 $ \(u,v) -> u `Sym.act` v == map (v!!) (Sym.toList u)++prop_act_id =+    forAll perm2 $ \(u,v) -> Sym.idperm u `Sym.act` v == v++prop_act_associative =+    forAll perm3 $ \(u,v,w) -> (u `Sym.act` v) `Sym.act` w == u `Sym.act` (v `Sym.act` w)++prop_size =+    forAll perm $ \v -> Sym.size v == Sym.size (Sym.st v)++prop_idperm =+    forAll perm2 $ \(u,v) -> Sym.idperm u == Sym.inverse (Sym.st u) `Sym.act` u++prop_inverse =+    forAll perm $ \v -> Sym.inverse v == Sym.inverse (Sym.st v) `Sym.act` Sym.idperm v++prop_ordiso1 =+    forAll perm2 $ \(u,v) -> u `Sym.ordiso` v  ==  (u == Sym.st v)++prop_ordiso2 =+    forAll perm2 $ \(u,v) -> u `Sym.ordiso` v  ==  (Sym.inverse u `Sym.act` v == Sym.idperm v)++prop_unrankPerm =+    forAll perm $ \w ->+    forAll (choose (0, product [1..fromIntegral (length w) - 1])) $ \r ->+        Sym.st (Sym.unrankPerm (sort w) r) == Sym.unrankStPerm (length w) r++prop_stackSort = forAll perm $ \v -> Sym.stackSort v == stack v++prop_stackSort_231 =+    forAll perm $ \v -> (Sym.stackSort v == Sym.idperm v) == (Sym.avoids [Sym.st "231"] v)++prop_bubbleSort = forAll perm $ \v -> Sym.bubbleSort v == bubble v++prop_bubbleSort_231_321 =+    forAll perm $ \v -> (Sym.bubbleSort v == Sym.idperm v) == (Sym.avoids [Sym.st "231", Sym.st "321"] v)++prop_subperm_copies p =+    forAll (resize 21 perm) $ \w -> and [ subperm m (Sym.st w) == p | m <- Sym.copies p w ]++prop_copies =+    forAll (resize  6 arbitrary) $ \p ->+    forAll (resize 12 perm)      $ \w ->+        sort (Sym.copies p w) == sort (map I.fromList $ copies (Sym.toList p) w)++prop_copies_self =+    forAll perm $ \v -> Sym.copies (Sym.st v) v == [SV.fromList [0 .. length v - 1]]++prop_copies_d8 (Symmetry (f,_)) =+    forAll (resize  6 arbitrary) $ \p ->+    forAll (resize 20 perm)      $ \w ->+        let p' = f p+            w' = Sym.generalize f w+        in length (Sym.copies p w) == length (Sym.copies p' w')++prop_avoiders_avoid =+    forAll (resize 20 arbitrary) $ \ws ->+    forAll (resize  6 arbitrary) $ \ps ->+        all (Sym.avoids ps) $ Sym.avoiders ps (ws :: [Sym.StPerm])++prop_avoiders_idempotent =+    forAll (resize 18 arbitrary) $ \vs ->+    forAll (resize  5 arbitrary) $ \ps ->+        let ws = Sym.avoiders ps (vs :: [Sym.StPerm]) in ws == Sym.avoiders ps ws++prop_avoiders_d8 (Symmetry (f,_)) =+    forAll (choose (0, 5))      $ \n ->+    forAll (resize 5 arbitrary) $ \p ->+        let ws = Sym.sym n in sort (map f $ Sym.avoiders [p] ws) == sort (Sym.avoiders [f p] ws)++prop_avoiders_d8' (Symmetry (f,_)) =+    forAll (choose (0, 5))      $ \n ->+    forAll (resize 5 arbitrary) $ \ps ->+        let ws = Sym.sym n in sort (map f $ Sym.avoiders ps ws) == sort (Sym.avoiders (map f ps) (map f ws))++prop_avoiders_d8'' (Symmetry (f,_)) =+    forAll (resize 18 arbitrary) $ \ws ->+    forAll (resize  5 arbitrary) $ \ps ->+        sort (map f $ Sym.avoiders ps ws) == sort (Sym.avoiders (map f ps) (map f ws :: [Sym.StPerm]))++prop_av_cardinality =+    forAll (resize 3 arbitrary) $ \p ->+        let spec = [ length $ Sym.av [p] n | n<-[0..6] ]+        in case Sym.size p of+             0 -> spec == [0,0,0,0,0,0,0]+             1 -> spec == [1,0,0,0,0,0,0]+             2 -> spec == [1,1,1,1,1,1,1]+             3 -> spec == [1,1,2,5,14,42,132]+             _ -> True++binomial n k = fromIntegral $ product [n', n'-1 .. n'-k'+1] `div` product [1..k']+    where+      n' = toInteger n+      k' = toInteger k++kSubsequences :: Int -> [a] -> [[a]]+kSubsequences 0 _      = [[]]+kSubsequences _ []     = []+kSubsequences k (x:xs) = map (x:) (kSubsequences (k-1) xs) ++ kSubsequences k xs++copies :: [Int] -> [Int] -> [[Int]]+copies p w = [ is | js <- u, let (is, q) = unzip (f js (zip [0..] w)), st q == p ]+    where+      k = length p+      n = length w+      u = kSubsequences k [0..n-1]+      f s@(j:t) ((i,x):v) = if i == j then (i,x) : f t v else f s v+      f _       _         = []++prop_subsets1 =+    forAll (choose (0,14)) $ \n ->+    forAll (choose (0,14)) $ \k ->+        sort (kSubsequences k [0..n-1]) == sort (map SV.toList $ Sym.subsets n k)++prop_subsets2 =+    forAll (choose (0,35)) $ \n ->+    forAll (choose (0,3))  $ \k ->+        sort (kSubsequences k [0..n-1]) == sort (map SV.toList $ Sym.subsets n k)++prop_subsets_singleton =+    forAll (choose (0,500)) $ \n ->+        let [v] = Sym.subsets n n in SV.toList v == [0..n-1]++prop_subsets_cardinality1 =+    forAll (choose (0,20)) $ \n ->+    forAll (choose (0,20)) $ \k ->+        length (Sym.subsets n k) == binomial n k++prop_subsets_cardinality2 =+    forAll (choose (0,20)) $ \n ->+    forAll (choose (0,20)) $ \k ->+        let cs = map (SV.length) (Sym.subsets n k) in ((k > n) && null cs) || ([k] == nub cs)++testsPerm =+    [ ("monoid/mempty/1",                check prop_monoid_mempty1)+    , ("monoid/mempty/2",                check prop_monoid_mempty2)+    , ("monoid/mempty/associative",      check prop_monoid_associative)+    , ("monoid/mempty/1/skew",           check prop_monoid_mempty1_S)+    , ("monoid/mempty/2/skew",           check prop_monoid_mempty2_S)+    , ("monoid/mempty/associative/skew", check prop_monoid_associative_S)+    , ("unrankStPerm/distinct",          check prop_unrankStPerm_distinct)+    , ("unrankStPerm/injective",         check prop_unrankStPerm_injective)+    , ("sym",                            check prop_sym)+    , ("perm",                           check prop_perm)+    , ("st",                             check prop_st)+    , ("act/def",                        check prop_act_def)+    , ("act/id",                         check prop_act_id)+    , ("act/associative",                check prop_act_associative)+    , ("size",                           check prop_size)+    , ("idperm",                         check prop_idperm)+    , ("inverse",                        check prop_inverse)+    , ("ordiso/1",                       check prop_ordiso1)+    , ("ordiso/2",                       check prop_ordiso2)+    , ("unrankPerm",                     check prop_unrankPerm)+    , ("stackSort",                      check prop_stackSort)+    , ("stackSort/231",                  check prop_stackSort_231)+    , ("bubbleSort",                     check prop_bubbleSort)+    , ("bubbleSort/231&321",             check prop_bubbleSort_231_321)+    , ("subperm/copies",                 check prop_subperm_copies)+    , ("copies",                         check prop_copies)+    , ("copies/self",                    check prop_copies_self)+    , ("copies/D8",                      check prop_copies_d8)+    , ("avoiders/avoid",                 check prop_avoiders_avoid)+    , ("avoiders/idempotent",            check prop_avoiders_idempotent)+    , ("avoiders/D8/0",                  check prop_avoiders_d8)+    , ("avoiders/D8/1",                  check prop_avoiders_d8')+    , ("avoiders/D8/2",                  check prop_avoiders_d8'')+    , ("av/cardinality",                 check prop_av_cardinality)+    , ("subsets/1",                      check prop_subsets1)+    , ("subsets/2",                      check prop_subsets2)+    , ("subsets/singleton",              check prop_subsets_singleton)+    , ("subsets/cardinality/1",          check prop_subsets_cardinality1)+    , ("subsets/cardinality/2",          check prop_subsets_cardinality2)+    ]++---------------------------------------------------------------------------------+-- Properties for Math.Sym.D8+---------------------------------------------------------------------------------++prop_D8_orbit fs w = all (`elem` orbD8) $ D8.orbit (map fn fs) w+    where+      orbD8 = D8.orbit D8.d8 w+      fn (Symmetry (f,_)) = f++prop_D8_reverse w    = I.reverse    (Sym.toVector w) == Sym.toVector (D8.reverse w)+prop_D8_complement w = I.complement (Sym.toVector w) == Sym.toVector (D8.complement w)+prop_D8_inverse w    = I.inverse    (Sym.toVector w) == Sym.toVector (D8.inverse w)+prop_D8_rotate w     = I.rotate     (Sym.toVector w) == Sym.toVector (D8.rotate w)++testsD8 =+    [ ("D8/orbit",       check prop_D8_orbit)+    , ("D8/reverse",     check prop_D8_reverse)+    , ("D8/complement",  check prop_D8_complement)+    , ("D8/inverse",     check prop_D8_inverse)+    , ("D8/rotate",      check prop_D8_rotate)+    ]++---------------------------------------------------------------------------------+-- Properties for Math.Sym.Stat+---------------------------------------------------------------------------------++-- the group theoretical inverse of w+inverse :: (Ord a) => [a] -> [Int]+inverse w = map snd . sort $ zip w [0..]++-- the standardization of w+st :: (Ord a) => [a] -> [Int]+st = inverse . inverse++ascents, descents :: (Ord a) => [a] -> [(a, a)]+ascents  w = filter (uncurry (<)) $ zip w (tail w)+descents w = filter (uncurry (>)) $ zip w (tail w)++peaks          w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x < y, y > z ]+valleys        w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x > y, y < z ]+doubleAscents  w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x < y, y < z ]+doubleDescents w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x > y, y > z ]++inversions :: (Ord a) => [a] -> [(a, a)]+inversions w = init (tails w) >>= \(x:xs) -> [ (x,y) | y<-xs, x > y ]++records :: (a -> a -> Bool) -> [a] -> [a]+records f []     = []+records f (x:xs) = records' f [x] xs where+    records' f recs       []     = recs+    records' f recs@(r:_) (x:xs) = records' f (if f r x then x:recs else recs) xs++lMinima, lMaxima, rMinima, rMaxima :: (Ord a) => [a] -> [a]++lMinima = reverse . records (>)+lMaxima = reverse . records (<)+rMinima = records (>) . reverse+rMaxima = records (<) . reverse++excedances  xs = map fst . filter (\(i,a)->i <  fromIntegral a) $ zip [0..] xs+fixedpoints xs = map fst . filter (\(i,a)->i == fromIntegral a) $ zip [0..] xs++exc, fp :: [Int] -> Int+exc = length . excedances . st+fp  = length . fixedpoints . st++runs :: Ord a => (a -> a -> Bool) -> [a] -> [a] -> [[a]]+runs _ [] [] = []+runs _ rs [] = [rs]+runs f [] (x:xs) = runs f [x] xs+runs f u@(r:_) v@(x:xs) | f r x = runs f (x:u) xs+                        | otherwise = u : runs f [x] xs++decruns :: Ord a => [a] -> [[a]]+decruns = runs (>) []++incruns :: Ord a => [a] -> [[a]]+incruns = runs (<) []++ldr, rdr, lir, rir :: (Ord a) => [a] -> Int++ldr [] = 0+ldr xs = length . head $ decruns xs++rdr [] = 0+rdr xs = length . last $ decruns xs++lir [] = 0+lir xs = length . head $ incruns xs++rir [] = 0+rir xs = length . last $ incruns xs++-- The stack-sort operator+stack [] = []+stack xs = stack left ++ stack right ++ [n]+    where+      (left, n:right) = span ( < maximum xs) xs++-- The bubble-sort operator; i.e. one pass of the classical bubble+-- sort algorithm+bubble :: Ord a => [a] -> [a]+bubble = bub []+    where+      bub xs []       = reverse xs+      bub [] (y:ys)   = bub [y] ys+      bub (x:xs) (y:ys)+          | x < y     = bub (y:x:xs) ys+          | otherwise = bub (x:y:xs) ys++-- Like Data.List.intersect, but by assuming that the lists are sorted+-- uses a faster algorithm+cap :: Ord a => [a] -> [a] -> [a]+cap [] ys = []+cap xs [] = []+cap xs@(x:xt) ys@(y:yt) = case compare x y of+                            EQ -> x : cap xt yt+                            LT -> cap xt ys+                            GT -> cap xs yt++-- The number of components in a permutation+comp w = length $ lMaxima w `cap` rMinima (bubble w)++-- rank a la Elizalde+ep = fst . last . filter (\(k,ys) -> all (k<=) ys) . zip [0..] . inits . st++des, asc, inv, lmin, lmax, rmin, rmax, peak, vall :: [Int] -> Int+dasc, ddes, maj, comp, ep :: [Int] -> Int++maj w = sum [ i | (i,x,y) <- zip3 [1..] w (tail w), x > y ]+des  = length . descents+asc  = length . ascents+inv  = length . inversions+lmin = length . lMinima+lmax = length . lMaxima+rmin = length . rMinima+rmax = length . rMaxima+peak = length . peaks+vall = length . valleys+dasc = length . doubleAscents+ddes = length . doubleDescents++prop_asc  = forAll perm $ \w -> asc  w == S.asc  w+prop_des  = forAll perm $ \w -> des  w == S.des  w+prop_exc  = forAll perm $ \w -> exc  w == S.exc  w+prop_fp   = forAll perm $ \w -> fp   w == S.fp   w+prop_inv  = forAll perm $ \w -> inv  w == S.inv  w+prop_maj  = forAll perm $ \w -> maj  w == S.maj  w+prop_lmin = forAll perm $ \w -> lmin w == S.lmin w+prop_lmax = forAll perm $ \w -> lmax w == S.lmax w+prop_rmin = forAll perm $ \w -> rmin w == S.rmin w+prop_rmax = forAll perm $ \w -> rmax w == S.rmax w+prop_head = forAll perm $ \w -> not (null w) ==> head (st w) == S.head w+prop_last = forAll perm $ \w -> not (null w) ==> last (st w) == S.last w+prop_peak = forAll perm $ \w -> peak w == S.peak w+prop_vall = forAll perm $ \w -> vall w == S.vall w+prop_dasc = forAll perm $ \w -> dasc w == S.dasc w+prop_ddes = forAll perm $ \w -> ddes w == S.ddes w+prop_ep   = forAll perm $ \w -> ep  w == S.ep  w+prop_lir  = forAll perm $ \w -> lir  w == S.lir  w+prop_ldr  = forAll perm $ \w -> ldr  w == S.ldr  w+prop_rir  = forAll perm $ \w -> rir  w == S.rir  w+prop_rdr  = forAll perm $ \w -> rdr  w == S.rdr  w+prop_comp = forAll perm $ \w -> comp w == S.comp w++prop_inv_21 = forAll perm $ \w -> S.inv w == length (Sym.copies (Sym.st "21") w)++testsStat =+    [ ("asc",    check prop_asc)+    , ("des",    check prop_des)+    , ("exc",    check prop_exc)+    , ("fp",     check prop_fp)+    , ("inv",    check prop_inv)+    , ("maj",    check prop_maj)+    , ("lmin",   check prop_lmin)+    , ("lmax",   check prop_lmax)+    , ("rmin",   check prop_rmin)+    , ("rmax",   check prop_rmax)+    , ("head",   check prop_head)+    , ("last",   check prop_last)+    , ("peak",   check prop_peak)+    , ("vall",   check prop_vall)+    , ("dasc",   check prop_dasc)+    , ("ddes",   check prop_ddes)+    , ("ep",     check prop_ep)+    , ("lir",    check prop_lir)+    , ("ldr",    check prop_ldr)+    , ("rir",    check prop_rir)+    , ("rdr",    check prop_rdr)+    , ("comp",   check prop_comp)+    , ("inv/21", check prop_inv_21)+    ]++---------------------------------------------------------------------------------+-- Main+---------------------------------------------------------------------------------++tests = testsPerm ++ testsD8 ++ testsStat++main = mapM_ (\(name, t) -> putStr (name ++ ":\t") >> t) tests