diff --git a/Math/Combinat.hs b/Math/Combinat.hs --- a/Math/Combinat.hs +++ b/Math/Combinat.hs @@ -1,7 +1,8 @@ --- | A collection of functions to generate and manipulate --- combinatorial objects like partitions, compositions, --- permutations, Young tableaux, various trees, etc etc. +-- | A collection of functions to generate, manipulate, +-- visualize and count combinatorial objects like partitions, +-- compositions, permutations, braids, Young tableaux, +-- lattice paths, various tree structures, etc etc. -- -- -- See also the @combinat-diagrams@ library for generating @@ -13,18 +14,22 @@ -- -- (1) generate most of the standard structures; -- --- (2) while being efficient; +-- (2) manipulate these structures; -- --- (3) to be able to enumerate the structures --- with constant memory usage; +-- (3) visualize these structures; -- --- (4) and to be able to randomly sample from them. +-- (4) the generation should be efficient; -- --- (5) finally, be a repository of algorithms +-- (5) to be able to enumerate the structures +-- with constant memory usage; -- +-- (6) to be able to randomly sample from them; +-- +-- (7) finally, to be a repository of algorithms. -- +-- -- The short-term goal is simply to generate --- many interesting structures. +-- and manipulate many interesting structures. -- -- -- Naming conventions (subject to change): @@ -40,7 +45,7 @@ -- * \"count\" prefix: counting functions. -- -- --- This module re-exports the most common modules. +-- This module re-exports the most commonly used modules. -- module Math.Combinat diff --git a/Math/Combinat/ASCII.hs b/Math/Combinat/ASCII.hs --- a/Math/Combinat/ASCII.hs +++ b/Math/Combinat/ASCII.hs @@ -11,12 +11,13 @@ -------------------------------------------------------------------------------- -import Data.List +import Data.Char ( isSpace ) +import Data.List ( transpose , intercalate ) import Math.Combinat.Helper -------------------------------------------------------------------------------- --- * The basic type +-- * The basic ASCII type -- | The type of a (rectangular) ASCII figure. -- Internally it is a list of lines of the same length plus the size. @@ -76,45 +77,7 @@ deriving (Eq,Show) data Alignment = Align HAlign VAlign - --------------------------------------------------------------------------------- --- * Extension --- | Extends an ASCII figure with spaces horizontally to the given width -hExtendTo :: HAlign -> Int -> ASCII -> ASCII -hExtendTo halign n0 rect@(ASCII (x,y) ls) = hExtendWith halign (max n0 x - x) rect - --- | Extends an ASCII figure with spaces vertically to the given height -vExtendTo :: VAlign -> Int -> ASCII -> ASCII -vExtendTo valign n0 rect@(ASCII (x,y) ls) = vExtendWith valign (max n0 y - y) rect - --- | Extend horizontally with the given number of spaces -hExtendWith :: HAlign -> Int -> ASCII -> ASCII -hExtendWith alignment d (ASCII (x,y) ls) = ASCII (x+d,y) (map f ls) where - f l = case alignment of - HLeft -> l ++ replicate d ' ' - HRight -> replicate d ' ' ++ l - HCenter -> replicate a ' ' ++ l ++ replicate (d-a) ' ' - a = div d 2 - --- | Extend vertically with the given number of empty lines -vExtendWith :: VAlign -> Int -> ASCII -> ASCII -vExtendWith valign d (ASCII (x,y) ls) = ASCII (x,y+d) (f ls) where - f ls = case valign of - VTop -> ls ++ replicate d emptyline - VBottom -> replicate d emptyline ++ ls - VCenter -> replicate a emptyline ++ ls ++ replicate (d-a) emptyline - a = div d 2 - emptyline = replicate x ' ' - --- | Horizontal indentation -hIndent :: Int -> ASCII -> ASCII -hIndent d = hExtendWith HRight d - --- | Vertical indentation -vIndent :: Int -> ASCII -> ASCII -vIndent d = vExtendWith VBottom d - -------------------------------------------------------------------------------- -- * Separators @@ -155,29 +118,35 @@ VSepEmpty -> [] VSepSpaces k -> replicate k ' ' VSepString s -> s - + -------------------------------------------------------------------------------- --- * Padding +-- * Concatenation --- | Horizontally pads with the given number of spaces, on both sides -hPad :: Int -> ASCII -> ASCII -hPad k (ASCII (x,y) ls) = ASCII (x+2*k,y) (map f ls) where - f l = pad ++ l ++ pad - pad = replicate k ' ' +-- | Horizontal append, centrally aligned, no separation. +(|||) :: ASCII -> ASCII -> ASCII +(|||) p q = hCatWith VCenter HSepEmpty [p,q] --- | Vertically pads with the given number of empty lines, on both sides -vPad :: Int -> ASCII -> ASCII -vPad k (ASCII (x,y) ls) = ASCII (x,y+2*k) (pad ++ ls ++ pad) where - pad = replicate k (replicate x ' ') +-- | Vertical append, centrally aligned, no separation. +(===) :: ASCII -> ASCII -> ASCII +(===) p q = vCatWith HCenter VSepEmpty [p,q] --- | Pads by single empty lines vertically and two spaces horizontally -pad :: ASCII -> ASCII -pad = vPad 1 . hPad 2 +-- | Horizontal concatenation, top-aligned, no separation +hCatTop :: [ASCII] -> ASCII +hCatTop = hCatWith VTop HSepEmpty --------------------------------------------------------------------------------- --- * Concatenation +-- | Horizontal concatenation, bottom-aligned, no separation +hCatBot :: [ASCII] -> ASCII +hCatBot = hCatWith VBottom HSepEmpty --- | Horizontal concatenation +-- | Vertical concatenation, left-aligned, no separation +vCatLeft :: [ASCII] -> ASCII +vCatLeft = vCatWith HLeft VSepEmpty + +-- | Vertical concatenation, right-aligned, no separation +vCatRight :: [ASCII] -> ASCII +vCatRight = vCatWith HRight VSepEmpty + +-- | General horizontal concatenation hCatWith :: VAlign -> HSep -> [ASCII] -> ASCII hCatWith valign hsep rects = ASCII (x',maxy) final where n = length rects @@ -189,7 +158,7 @@ x' = sum' xsz + (n-1)*sepx final = map (intercalate sep) $ transpose (map asciiLines rects1) --- | Vertical concatenation +-- | General vertical concatenation vCatWith :: HAlign -> VSep -> [ASCII] -> ASCII vCatWith halign vsep rects = ASCII (maxx,y') final where n = length rects @@ -202,8 +171,163 @@ final = intercalate fullsep $ map asciiLines rects1 -------------------------------------------------------------------------------- +-- * Padding + +-- | Horizontally pads with the given number of spaces, on both sides +hPad :: Int -> ASCII -> ASCII +hPad k (ASCII (x,y) ls) = ASCII (x+2*k,y) (map f ls) where + f l = pad ++ l ++ pad + pad = replicate k ' ' + +-- | Vertically pads with the given number of empty lines, on both sides +vPad :: Int -> ASCII -> ASCII +vPad k (ASCII (x,y) ls) = ASCII (x,y+2*k) (pad ++ ls ++ pad) where + pad = replicate k (replicate x ' ') + +-- | Pads by single empty lines vertically and two spaces horizontally +pad :: ASCII -> ASCII +pad = vPad 1 . hPad 2 + +-------------------------------------------------------------------------------- +-- * Extension + +-- | Extends an ASCII figure with spaces horizontally to the given width. +-- Note: the alignment is the alignment of the original picture in the new bigger picture! +hExtendTo :: HAlign -> Int -> ASCII -> ASCII +hExtendTo halign n0 rect@(ASCII (x,y) ls) = hExtendWith halign (max n0 x - x) rect + +-- | Extends an ASCII figure with spaces vertically to the given height. +-- Note: the alignment is the alignment of the original picture in the new bigger picture! +vExtendTo :: VAlign -> Int -> ASCII -> ASCII +vExtendTo valign n0 rect@(ASCII (x,y) ls) = vExtendWith valign (max n0 y - y) rect + +-- | Extend horizontally with the given number of spaces. +hExtendWith :: HAlign -> Int -> ASCII -> ASCII +hExtendWith alignment d (ASCII (x,y) ls) = ASCII (x+d,y) (map f ls) where + f l = case alignment of + HLeft -> l ++ replicate d ' ' + HRight -> replicate d ' ' ++ l + HCenter -> replicate a ' ' ++ l ++ replicate (d-a) ' ' + a = div d 2 + +-- | Extend vertically with the given number of empty lines. +vExtendWith :: VAlign -> Int -> ASCII -> ASCII +vExtendWith valign d (ASCII (x,y) ls) = ASCII (x,y+d) (f ls) where + f ls = case valign of + VTop -> ls ++ replicate d emptyline + VBottom -> replicate d emptyline ++ ls + VCenter -> replicate a emptyline ++ ls ++ replicate (d-a) emptyline + a = div d 2 + emptyline = replicate x ' ' + +-- | Horizontal indentation +hIndent :: Int -> ASCII -> ASCII +hIndent d = hExtendWith HRight d + +-- | Vertical indentation +vIndent :: Int -> ASCII -> ASCII +vIndent d = vExtendWith VBottom d + +-------------------------------------------------------------------------------- +-- * Cutting + +-- | Cuts the given number of columns from the picture. +-- The alignment is the alignment of the /picture/, not the cuts. +-- +-- This should be the (left) inverse of 'hExtendWith'. +hCut :: HAlign -> Int -> ASCII -> ASCII +hCut halign k (ASCII (x,y) ls) = ASCII (x',y) (map f ls) where + x' = max 0 (x-k) + f = case halign of + HLeft -> reverse . drop k . reverse + HCenter -> reverse . drop (k-a) . reverse . drop a + HRight -> drop k + a = div k 2 + +-- | Cuts the given number of rows from the picture. +-- The alignment is the alignment of the /picture/, not the cuts. +-- +-- This should be the (left) inverse of 'vExtendWith'. +vCut :: VAlign -> Int -> ASCII -> ASCII +vCut valign k (ASCII (x,y) ls) = ASCII (x,y') (g ls) where + y' = max 0 (y-k) + g = case valign of + VTop -> reverse . drop k . reverse + VCenter -> reverse . drop (k-a) . reverse . drop a + VBottom -> drop k + a = div k 2 + +-------------------------------------------------------------------------------- +-- * Pasting + +-- | Pastes the first ASCII graphics onto the second, keeping the second one's dimension +-- (that is, overlapping parts of the first one are ignored). +-- The offset is relative to the top-left corner of the second picture. +-- Spaces at treated as transparent. +-- +-- Example: +-- +-- > tabulate (HCenter,VCenter) (HSepSpaces 2, VSepSpaces 1) +-- > [ [ caption (show (x,y)) $ +-- > pasteOnto (x,y) (filledBox '@' (4,3)) (asciiBox (7,5)) +-- > | x <- [-4..7] ] +-- > | y <- [-3..5] ] +-- +pasteOnto :: (Int,Int) -> ASCII -> ASCII -> ASCII +pasteOnto = pasteOnto' isSpace + +-- | Pastes the first ASCII graphics onto the second, keeping the second one's dimension. +-- The first argument specifies the transparency condition (on the first picture). +-- The offset is relative to the top-left corner of the second picture. +-- +pasteOnto' + :: (Char -> Bool) -- ^ transparency condition + -> (Int,Int) -- ^ offset relative to the top-left corner of the second picture + -> ASCII -- ^ picture to paste + -> ASCII -- ^ picture to paste onto + -> ASCII +pasteOnto' transparent (xpos,ypos) small big = new where + new = ASCII (xbig,ybig) lines' + (xbig,ybig) = asciiSize big + bigLines = asciiLines big + small' = (if (ypos>=0) then vExtendWith VBottom ypos else vCut VBottom (-ypos)) + $ (if (xpos>=0) then hExtendWith HRight xpos else hCut HRight (-xpos)) + $ small + smallLines = asciiLines small' + lines' = zipWith f bigLines (smallLines ++ repeat "") + f bl sl = zipWith g bl (sl ++ repeat ' ') + g b s = if transparent s then b else s + +-- | A version of 'pasteOnto' where we can specify the corner of the second picture +-- to which the offset is relative: +-- +-- > pasteOntoRel (HLeft,VTop) == pasteOnto +-- +pasteOntoRel :: (HAlign,VAlign) -> (Int,Int) -> ASCII -> ASCII -> ASCII +pasteOntoRel = pasteOntoRel' isSpace + +pasteOntoRel' :: (Char -> Bool) -> (HAlign,VAlign) -> (Int,Int) -> ASCII -> ASCII -> ASCII +pasteOntoRel' transparent (halign,valign) (xpos,ypos) small big = new where + new = pasteOnto' transparent (xpos',ypos') small big + (xsize,ysize) = asciiSize big + xpos' = case halign of + HLeft -> xpos + HCenter -> xpos + div xsize 2 + HRight -> xpos + xsize + ypos' = case valign of + VTop -> ypos + VCenter -> ypos + div ysize 2 + VBottom -> ypos + ysize + +-------------------------------------------------------------------------------- -- * Tabulate +-- | Tabulates the given matrix of pictures. Example: +-- +-- > tabulate (HCenter, VCenter) (HSepSpaces 2, VSepSpaces 1) +-- > [ [ asciiFromLines [ "x=" ++ show x , "y=" ++ show y ] | x<-[7..13] ] +-- > | y<-[98..102] ] +-- tabulate :: (HAlign,VAlign) -> (HSep,VSep) -> [[ASCII]] -> ASCII tabulate (halign,valign) (hsep,vsep) rects0 = final where n = length rects0 @@ -226,7 +350,7 @@ -- autoTabulate :: MatrixOrder -- ^ whether to use row-major or column-major ordering of the elements - -> Either Int Int -- ^ @(Right x)@ creates x columns, while @(Left y)$ creates y rows + -> Either Int Int -- ^ @(Right x)@ creates x columns, while @(Left y)@ creates y rows -> [ASCII] -- ^ list of ASCII rectangles -> ASCII autoTabulate mtxorder ei list = final where @@ -276,21 +400,35 @@ capt = asciiFromString str -------------------------------------------------------------------------------- --- * Testing \/ miscellanea +-- * Ready-made boxes --- | An ASCII box of the given size +-- | An ASCII border box of the given size asciiBox :: (Int,Int) -> ASCII asciiBox (x,y) = ASCII (max x 2, max y 2) (h : replicate (y-2) m ++ [h]) where h = "+" ++ replicate (x-2) '-' ++ "+" m = "|" ++ replicate (x-2) ' ' ++ "|" --- | An \"rounded\" ASCII box of the given size +-- | An \"rounded\" ASCII border box of the given size roundedAsciiBox :: (Int,Int) -> ASCII roundedAsciiBox (x,y) = ASCII (max x 2, max y 2) (a : replicate (y-2) m ++ [b]) where a = "/" ++ replicate (x-2) '-' ++ "\\" m = "|" ++ replicate (x-2) ' ' ++ "|" b = "\\" ++ replicate (x-2) '-' ++ "/" +-- | A box simply filled with the given character +filledBox :: Char -> (Int,Int) -> ASCII +filledBox c (x0,y0) = asciiFromLines $ replicate y (replicate x c) where + x = max 0 x0 + y = max 0 y0 + +-- | A box of spaces +transparentBox :: (Int,Int) -> ASCII +transparentBox = filledBox ' ' + +-------------------------------------------------------------------------------- +-- * Testing \/ miscellanea + +-- | An integer asciiNumber :: Int -> ASCII asciiNumber = asciiShow diff --git a/Math/Combinat/Classes.hs b/Math/Combinat/Classes.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Classes.hs @@ -0,0 +1,66 @@ + +-- | Type classes for some common properties shared by different objects + +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} +module Math.Combinat.Classes where + +-------------------------------------------------------------------------------- + +-- | Emptyness +class CanBeEmpty a where + isEmpty :: a -> Bool + empty :: a + +-------------------------------------------------------------------------------- +-- * Partitions + +-- | Number of parts +class HasNumberOfParts a where + numberOfParts :: a -> Int + +-------------------------------------------------------------------------------- + +class HasWidth a where + width :: a -> Int + +class HasHeight a where + height :: a -> Int + +-------------------------------------------------------------------------------- + +-- | Weight (of partitions, tableaux, etc) +class HasWeight a where + weight :: a -> Int + +-------------------------------------------------------------------------------- + +-- | Duality (of partitions, tableaux, etc) +class HasDuality a where + dual :: a -> a + +-------------------------------------------------------------------------------- +-- * Tableau + +-- | Shape (of tableaux, skew tableaux) +class HasShape a s | a -> s where + shape :: a -> s + +-------------------------------------------------------------------------------- +-- * Trees + +-- | Number of nodes (of trees) +class HasNumberOfNodes t where + numberOfNodes :: t -> Int + +-- | Number of leaves (of trees) +class HasNumberOfLeaves t where + numberOfLeaves :: t -> Int + +-------------------------------------------------------------------------------- +-- * Permutations + +-- | Number of cycles (of partitions) +class HasNumberOfCycles p where + numberOfCycles :: p -> Int + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/FreeGroups.hs b/Math/Combinat/FreeGroups.hs deleted file mode 100644 --- a/Math/Combinat/FreeGroups.hs +++ /dev/null @@ -1,376 +0,0 @@ - --- | Words in free groups (and free powers of cyclic groups). --- This module is not re-exported by "Math.Combinat" --- -{-# LANGUAGE CPP, PatternGuards #-} -module Math.Combinat.FreeGroups where - --------------------------------------------------------------------------------- - --- new Base exports "Word" from Data.Word... -#ifdef MIN_VERSION_base -#if MIN_VERSION_base(4,7,1) -import Prelude hiding ( Word ) -#endif -#elif __GLASGOW_HASKELL__ >= 709 -import Prelude hiding ( Word ) -#endif - -import Data.Char ( chr ) -import Data.List ( mapAccumL ) - -import Control.Monad ( liftM ) -import System.Random - -import Math.Combinat.Numbers -import Math.Combinat.Helper - --------------------------------------------------------------------------------- --- * Words - --- | A generator of a (free) group -data Generator a - = Gen a -- @a@ - | Inv a -- @a^(-1)@ - deriving (Eq,Ord,Show,Read) - --- | The index of a generator -unGen :: Generator a -> a -unGen g = case g of - Gen x -> x - Inv x -> x - --- | A /word/, describing (non-uniquely) an element of a group. --- The identity element is represented (among others) by the empty word. -type Word a = [Generator a] - --------------------------------------------------------------------------------- - --- | Generators are shown as small letters: @a@, @b@, @c@, ... --- and their inverses are shown as capital letters, so @A=a^-1@, @B=b^-1@, etc. -showGen :: Generator Int -> Char -showGen (Gen i) = chr (96+i) -showGen (Inv i) = chr (64+i) - -showWord :: Word Int -> String -showWord = map showGen - --------------------------------------------------------------------------------- - -instance Functor Generator where - fmap f g = case g of - Gen x -> Gen (f x) - Inv y -> Inv (f y) - --------------------------------------------------------------------------------- - --- | The inverse of a generator -inverseGen :: Generator a -> Generator a -inverseGen g = case g of - Gen x -> Inv x - Inv x -> Gen x - --- | The inverse of a word -inverseWord :: Word a -> Word a -inverseWord = map inverseGen . reverse - --- | Lists all words of the given length (total number will be @(2g)^n@). --- The numbering of the generators is @[1..g]@. -allWords - :: Int -- ^ @g@ = number of generators - -> Int -- ^ @n@ = length of the word - -> [Word Int] -allWords g = go where - go 0 = [[]] - go n = [ x:xs | xs <- go (n-1) , x <- elems ] - elems = [ Gen a | a<-[1..g] ] - ++ [ Inv a | a<-[1..g] ] - --- | Lists all words of the given length which do not contain inverse generators --- (total number will be @g^n@). --- The numbering of the generators is @[1..g]@. -allWordsNoInv - :: Int -- ^ @g@ = number of generators - -> Int -- ^ @n@ = length of the word - -> [Word Int] -allWordsNoInv g = go where - go 0 = [[]] - go n = [ x:xs | xs <- go (n-1) , x <- elems ] - elems = [ Gen a | a<-[1..g] ] - --------------------------------------------------------------------------------- --- * Random words - --- | A random group generator (or its inverse) between @1@ and @g@ -randomGenerator - :: RandomGen g - => Int -- ^ @g@ = number of generators - -> g -> (Generator Int, g) -randomGenerator d g0 = (gen,g2) where - (b,g1) = random g0 - (k,g2) = randomR (1,d) g1 - gen = if b then Gen k else Inv k - --- | A random group generator (but never its inverse) between @1@ and @g@ -randomGeneratorNoInv - :: RandomGen g - => Int -- ^ @g@ = number of generators - -> g -> (Generator Int, g) -randomGeneratorNoInv d g0 = (Gen k,g1) where - (k,g1) = randomR (1,d) g0 - --- | A random word of length @n@ using @g@ generators (or their inverses) -randomWord - :: RandomGen g - => Int -- ^ @g@ = number of generators - -> Int -- ^ @n@ = length of the word - -> g -> (Word Int, g) -randomWord d n g0 = (word,g1) where - (g1,word) = mapAccumL (\g _ -> swap (randomGenerator d g)) g0 [1..n] - --- | A random word of length @n@ using @g@ generators (but not their inverses) -randomWordNoInv - :: RandomGen g - => Int -- ^ @g@ = number of generators - -> Int -- ^ @n@ = length of the word - -> g -> (Word Int, g) -randomWordNoInv d n g0 = (word,g1) where - (g1,word) = mapAccumL (\g _ -> swap (randomGeneratorNoInv d g)) g0 [1..n] - --------------------------------------------------------------------------------- --- * The free group on @g@ generators - --- | Multiplication of the free group (returns the reduced result). It is true --- for any two words w1 and w2 that --- --- > multiplyFree (reduceWordFree w1) (reduceWord w2) = multiplyFree w1 w2 --- -multiplyFree :: Eq a => Word a -> Word a -> Word a -multiplyFree w1 w2 = reduceWordFree (w1++w2) - --- | Reduces a word in a free group by repeatedly removing @x*x^(-1)@ and --- @x^(-1)*x@ pairs. The set of /reduced words/ forms the free group; the --- multiplication is obtained by concatenation followed by reduction. --- -reduceWordFree :: Eq a => Word a -> Word a -reduceWordFree = loop where - - loop w = case reduceStep w of - Nothing -> w - Just w' -> loop w' - - reduceStep :: Eq a => Word a -> Maybe (Word a) - reduceStep = go False where - go changed w = case w of - (Gen x : Inv y : rest) | x==y -> go True rest - (Inv x : Gen y : rest) | x==y -> go True rest - (this : rest) -> liftM (this:) $ go changed rest - _ -> if changed then Just w else Nothing - --------------------------------------------------------------------------------- - --- | Counts the number of words of length @n@ which reduce to the identity element. --- --- Generating function is @Gf_g(u) = \\frac {2g-1} { g-1 + g \\sqrt{ 1 - (8g-4)u^2 } }@ --- -countIdentityWordsFree - :: Int -- ^ g = number of generators in the free group - -> Int -- ^ n = length of the unreduced word - -> Integer -countIdentityWordsFree g n = countWordReductionsFree g n 0 - --- | Counts the number of words of length @n@ whose reduced form has length @k@ --- (clearly @n@ and @k@ must have the same parity for this to be nonzero): --- --- > countWordReductionsFree g n k == sum [ 1 | w <- allWords g n, k == length (reduceWordFree w) ] --- -countWordReductionsFree - :: Int -- ^ g = number of generators in the free group - -> Int -- ^ n = length of the unreduced word - -> Int -- ^ k = length of the reduced word - -> Integer -countWordReductionsFree gens_ nn_ kk_ - | nn==0 = if k==0 then 1 else 0 - | even nn && kk == 0 = sum [ ( binomial (nn-i) (n -i) * gg^(i ) * (gg-1)^(n -i ) * ( i) ) `div` (nn-i) | i<-[0..n ] ] - | even nn && even kk = sum [ ( binomial (nn-i) (n-k-i) * gg^(i+1) * (gg-1)^(n+k-i-1) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] - | odd nn && odd kk = sum [ ( binomial (nn-i) (n-k-i) * gg^(i+1) * (gg-1)^(n+k-i ) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] - | otherwise = 0 - where - g = fromIntegral gens_ :: Integer - nn = fromIntegral nn_ :: Integer - kk = fromIntegral kk_ :: Integer - - gg = 2*g - n = div nn 2 - k = div kk 2 - --------------------------------------------------------------------------------- --- * Free powers of cyclic groups - --- | Multiplication in free products of Z2's -multiplyZ2 :: Eq a => Word a -> Word a -> Word a -multiplyZ2 w1 w2 = reduceWordZ2 (w1++w2) - --- | Multiplication in free products of Z3's -multiplyZ3 :: Eq a => Word a -> Word a -> Word a -multiplyZ3 w1 w2 = reduceWordZ3 (w1++w2) - --- | Multiplication in free products of Zm's -multiplyZm :: Eq a => Int -> Word a -> Word a -> Word a -multiplyZm k w1 w2 = reduceWordZm k (w1++w2) - --------------------------------------------------------------------------------- - --- | Reduces a word, where each generator @x@ satisfies the additional relation @x^2=1@ --- (that is, free products of Z2's) -reduceWordZ2 :: Eq a => Word a -> Word a -reduceWordZ2 = loop where - loop w = case reduceStep w of - Nothing -> w - Just w' -> loop w' - - reduceStep :: Eq a => Word a -> Maybe (Word a) - reduceStep = go False where - go changed w = case w of - (Gen x : Gen y : rest) | x==y -> go True rest - (Gen x : Inv y : rest) | x==y -> go True rest - (Inv x : Gen y : rest) | x==y -> go True rest - (Inv x : Inv y : rest) | x==y -> go True rest - (this : rest) -> liftM (this:) $ go changed rest - _ -> if changed then Just w else Nothing - --- | Reduces a word, where each generator @x@ satisfies the additional relation @x^3=1@ --- (that is, free products of Z3's) -reduceWordZ3 :: Eq a => Word a -> Word a -reduceWordZ3 = loop where - loop w = case reduceStep w of - Nothing -> w - Just w' -> loop w' - - reduceStep :: Eq a => Word a -> Maybe (Word a) - reduceStep = go False where - go changed w = case w of - (Gen x : Inv y : rest) | x==y -> go True rest - (Inv x : Gen y : rest) | x==y -> go True rest - (Gen x : Gen y : Gen z : rest) | x==y && y==z -> go True rest - (Inv x : Inv y : Inv z : rest) | x==y && y==z -> go True rest - (Gen x : Gen y : rest) | x==y -> go True (Inv x : rest) -- !!! - (Inv x : Inv y : rest) | x==y -> go True (Gen x : rest) - (this : rest) -> liftM (this:) $ go changed rest - _ -> if changed then Just w else Nothing - --- | Reduces a word, where each generator @x@ satisfies the additional relation @x^m=1@ --- (that is, free products of Zm's) -reduceWordZm :: Eq a => Int -> Word a -> Word a -reduceWordZm m = loop where - - loop w = case reduceStep w of - Nothing -> w - Just w' -> loop w' - - halfm = div m 2 -- if we encounter strictly more than m/2 equal elements in a row, we replace them by the inverses - - reduceStep :: Eq a => Word a -> Maybe (Word a) - reduceStep = go False where - go changed w = case w of - (Gen x : Inv y : rest) | x==y -> go True rest - (Inv x : Gen y : rest) | x==y -> go True rest --- something | Just rest <- dropk w -> go True rest - something | Just (k,rest) <- dropIfMoreThanHalf w -> go True (replicate (m-k) (inverseGen (head w)) ++ rest) - (this : rest) -> liftM (this:) $ go changed rest - _ -> if changed then Just w else Nothing - - dropIfMoreThanHalf :: Eq a => Word a -> Maybe (Int, Word a) - dropIfMoreThanHalf w = - let (k,rest) = dropWhileEqual w - in if k > halfm then Just (k,rest) - else Nothing - - dropWhileEqual :: Eq a => Word a -> (Int, Word a) - dropWhileEqual [] = (0,[]) - dropWhileEqual (x0:rest) = go 1 rest where - go k [] = (k,[]) - go k xxs@(x:xs) = if k==m then (m,xxs) - else if x==x0 then go (k+1) xs - else (k,xxs) - -{- - dropm :: Eq a => Word a -> Maybe (Word a) - dropm [] = Nothing - dropm (x:xs) = go (m-1) xs where - go 0 rest = Just rest - go j (y:ys) = if y==x - then go (j-1) ys - else Nothing - go j [] = Nothing --} - --------------------------------------------------------------------------------- - --- | Counts the number of words (without inverse generators) of length @n@ --- which reduce to the identity element, using the relations @x^2=1@. --- --- Generating function is @Gf_g(u) = \\frac {2g-2} { g-2 + g \\sqrt{ 1 - (4g-4)u^2 } }@ --- --- The first few @g@ cases: --- --- > A000984 = [ countIdentityWordsZ2 2 (2*n) | n<-[0..] ] = [1,2,6,20,70,252,924,3432,12870,48620,184756...] --- > A089022 = [ countIdentityWordsZ2 3 (2*n) | n<-[0..] ] = [1,3,15,87,543,3543,23823,163719,1143999,8099511,57959535...] --- > A035610 = [ countIdentityWordsZ2 4 (2*n) | n<-[0..] ] = [1,4,28,232,2092,19864,195352,1970896,20275660,211823800,2240795848...] --- > A130976 = [ countIdentityWordsZ2 5 (2*n) | n<-[0..] ] = [1,5,45,485,5725,71445,925965,12335685,167817405,2321105525,32536755565...] --- -countIdentityWordsZ2 - :: Int -- ^ g = number of generators in the free group - -> Int -- ^ n = length of the unreduced word - -> Integer -countIdentityWordsZ2 g n = countWordReductionsZ2 g n 0 - --- | Counts the number of words (without inverse generators) of length @n@ whose --- reduced form in the product of Z2-s (that is, for each generator @x@ we have @x^2=1@) --- has length @k@ --- (clearly @n@ and @k@ must have the same parity for this to be nonzero): --- --- > countWordReductionsZ2 g n k == sum [ 1 | w <- allWordsNoInv g n, k == length (reduceWordZ2 w) ] --- -countWordReductionsZ2 - :: Int -- ^ g = number of generators in the free group - -> Int -- ^ n = length of the unreduced word - -> Int -- ^ k = length of the reduced word - -> Integer -countWordReductionsZ2 gens_ nn_ kk_ - | nn==0 = if k==0 then 1 else 0 - | even nn && kk == 0 = sum [ ( binomial (nn-i) (n -i) * g^(i ) * (g-1)^(n -i ) * ( i) ) `div` (nn-i) | i<-[0..n ] ] - | even nn && even kk = sum [ ( binomial (nn-i) (n-k-i) * g^(i+1) * (g-1)^(n+k-i-1) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] - | odd nn && odd kk = sum [ ( binomial (nn-i) (n-k-i) * g^(i+1) * (g-1)^(n+k-i ) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] - | otherwise = 0 - where - g = fromIntegral gens_ :: Integer - nn = fromIntegral nn_ :: Integer - kk = fromIntegral kk_ :: Integer - - n = div nn 2 - k = div kk 2 - --- | Counts the number of words (without inverse generators) of length @n@ --- which reduce to the identity element, using the relations @x^3=1@. --- --- > countIdentityWordsZ3NoInv g n == sum [ 1 | w <- allWordsNoInv g n, 0 == length (reduceWordZ2 w) ] --- --- In mathematica, the formula is: @Sum[ g^k * (g-1)^(n-k) * k/n * Binomial[3*n-k-1, n-k] , {k, 1,n} ]@ --- -countIdentityWordsZ3NoInv - :: Int -- ^ g = number of generators in the free group - -> Int -- ^ n = length of the unreduced word - -> Integer -countIdentityWordsZ3NoInv gens_ nn_ - | nn==0 = 1 - | mod nn 3 == 0 = sum [ ( binomial (3*n-i-1) (n-i) * g^i * (g-1)^(n-i) * i ) `div` n | i<-[1..n] ] - | otherwise = 0 - where - g = fromIntegral gens_ :: Integer - nn = fromIntegral nn_ :: Integer - - n = div nn 3 - --------------------------------------------------------------------------------- - diff --git a/Math/Combinat/Groups/Braid.hs b/Math/Combinat/Groups/Braid.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Groups/Braid.hs @@ -0,0 +1,651 @@ + +-- | Braids. See eg. +-- +-- +-- Based on: +-- +-- * Joan S. Birman, Tara E. Brendle: BRAIDS - A SURVEY +-- +-- +-- +-- Note: This module GHC 7.8, since we use type-level naturals +-- to parametrize the 'Braid' type. +-- + + +{-# LANGUAGE + CPP, BangPatterns, + ScopedTypeVariables, ExistentialQuantification, + DataKinds, KindSignatures, Rank2Types, + TypeOperators, TypeFamilies #-} + +module Math.Combinat.Groups.Braid where + +-------------------------------------------------------------------------------- + +import Data.Proxy +import GHC.TypeLits + +import Control.Monad + +import Data.List ( mapAccumL , foldl' ) + +import Data.Array.Unboxed +import Data.Array.ST +import Data.Array.IArray +import Data.Array.MArray +import Data.Array.Unsafe +import Data.Array.Base + +import Control.Monad.ST + +import System.Random + +import Math.Combinat.ASCII +import Math.Combinat.Sign +import Math.Combinat.Helper + +import Math.Combinat.Permutations ( Permutation(..) ) +import qualified Math.Combinat.Permutations as P + +#ifdef QUICKCHECK +import Test.QuickCheck +#endif + +-------------------------------------------------------------------------------- +-- * Artin generators + +-- | A standard Artin generator of a braid: @Sigma i@ represents twisting +-- the neighbour strands @i@ and @(i+1)@, such that @#i@ goes /under/ @#(i+1)@ +-- +-- Note: The strands are numbered @1..n@. +data BrGen + = Sigma !Int + | SigmaInv !Int + deriving (Eq,Ord,Show) + +-- | The strand (more precisely, the first of the two strands) the generator twistes +brGenIdx :: BrGen -> Int +brGenIdx g = case g of + Sigma i -> i + SigmaInv i -> i + +brGenSign :: BrGen -> Sign +brGenSign g = case g of + Sigma _ -> Plus + SigmaInv _ -> Minus + +brGenSignIdx :: BrGen -> (Sign,Int) +brGenSignIdx g = case g of + Sigma i -> (Plus ,i) + SigmaInv i -> (Minus,i) + +-- | The inverse of a braid generator +invBrGen :: BrGen -> BrGen +invBrGen g = case g of + Sigma i -> SigmaInv i + SigmaInv i -> Sigma i + +-------------------------------------------------------------------------------- +-- * The braid type + +-- | The braid group @B_n@ on @n@ strands. +-- The number @n@ is encoded as a type level natural in the type parameter. +-- +-- Braids are represented as words in the standard generators and their +-- inverses. +newtype Braid (n :: Nat) = Braid [BrGen] deriving (Show) + +-- | The number of strands in the braid +numberOfStrands :: KnownNat n => Braid n -> Int +numberOfStrands = fromInteger . natVal . braidProxy where + braidProxy :: Braid n -> Proxy n + braidProxy _ = Proxy + +-- | Sometimes we want to hide the type-level parameter @n@, for example when +-- dynamically creating braids whose size is known only at runtime. +data SomeBraid = forall n. KnownNat n => SomeBraid (Braid n) + +someBraid :: Int -> (forall (n :: Nat). KnownNat n => Braid n) -> SomeBraid +someBraid n polyBraid = + case snat of + SomeNat pxy -> SomeBraid (asProxyTypeOf1 polyBraid pxy) + where + snat = case someNatVal (fromIntegral n :: Integer) of + Just sn -> sn + Nothing -> error "someBraid: input is not a natural number" + +withSomeBraid :: SomeBraid -> (forall n. KnownNat n => Braid n -> a) -> a +withSomeBraid sbraid f = case sbraid of SomeBraid braid -> f braid + +-------------------------------------------------------------------------------- + +braidWord :: Braid n -> [BrGen] +braidWord (Braid gs) = gs + +braidWordLength :: Braid n -> Int +braidWordLength (Braid gs) = length gs + +-- | Embeds a smaller braid group into a bigger braid group +extend :: (n1 <= n2) => Braid n1 -> Braid n2 +extend (Braid gs) = Braid gs + +-- | Apply \"free reduction\" to the word, that is, iteratively remove @sigma_i sigma_i^-1@ pairs. +-- The resulting braid is clearly equivalent to the original. +freeReduceBraidWord :: Braid n -> Braid n +freeReduceBraidWord (Braid orig) = Braid (loop orig) where + + loop w = case reduceStep w of + Nothing -> w + Just w' -> loop w' + + reduceStep :: [BrGen] -> Maybe [BrGen] + reduceStep = go False where + go !changed w = case w of + (Sigma x : SigmaInv y : rest) | x==y -> go True rest + (SigmaInv x : Sigma y : rest) | x==y -> go True rest + (this : rest) -> liftM (this:) $ go changed rest + _ -> if changed then Just w else Nothing + +-------------------------------------------------------------------------------- +-- * Some specific braids + +-- | The braid generator @sigma_i@ as a braid +sigma :: KnownNat n => Int -> Braid (n :: Nat) +sigma k = braid where + braid = if k > 0 && k < numberOfStrands braid + then Braid [Sigma k] + else error "sigma: braid generator index out of range" + +-- | The braid generator @sigma_i^(-1)@ as a braid +sigmaInv :: KnownNat n => Int -> Braid (n :: Nat) +sigmaInv k = braid where + braid = if k > 0 && k < numberOfStrands braid + then Braid [SigmaInv k] + else error "sigma: braid generator index out of range" + +-- | @doubleSigma s t@ (for s Int -> Int -> Braid (n :: Nat) +doubleSigma s t = braid where + n = numberOfStrands braid + braid + | s < 1 || s > n = error "doubleSigma: s index out of range" + | t < 1 || t > n = error "doubleSigma: t index out of range" + | s >= t = error "doubleSigma: s >= t" + | otherwise = Braid $ + [ Sigma i | i<-[t-1,t-2..s] ] ++ [ SigmaInv i | i<-[s+1..t-1] ] + +-- | @positiveWord [2,5,1]@ is shorthand for the word @sigma_2*sigma_5*sigma_1@. +positiveWord :: KnownNat n => [Int] -> Braid (n :: Nat) +positiveWord idxs = braid where + braid = Braid (map gen idxs) + n = numberOfStrands braid + gen i = if i>0 && i Braid n +halfTwist = braid where + braid = Braid $ map Sigma $ _halfTwist n + n = numberOfStrands braid + +-- | The untyped version of 'halfTwist' +_halfTwist :: Int -> [Int] +_halfTwist n = gens where + gens = concat [ sub k | k<-[1..n-1] ] + sub k = [ j | j<-[n-1,n-2..k] ] + +-- | Synonym for 'halfTwist' +theGarsideBraid :: KnownNat n => Braid n +theGarsideBraid = halfTwist + +-- | The inner automorphism defined by @X -> Delta^-1 X Delta@, +-- where @Delta@ is the positive half-twist. +-- +-- This sends each generator @sigma_j@ to @sigma_(n-j)@. +-- +tau :: KnownNat n => Braid n -> Braid n +tau braid@(Braid gens) = Braid (map f gens) where + n = numberOfStrands braid + f (Sigma i) = Sigma (n-i) + f (SigmaInv i) = SigmaInv (n-i) + +-------------------------------------------------------------------------------- +-- * Group operations + +-- | The trivial braid +identity :: Braid n +identity = Braid [] + +-- | The inverse of a braid. Note: we do not perform reduction here, +-- as a word is reduced if and only if its inverse is reduced. +inverse :: Braid n -> Braid n +inverse = Braid . reverse . map invBrGen . braidWord + +-- | Composes two braids, doing free reduction on the result +-- (that is, removing @(sigma_k * sigma_k^-1)@ pairs@) +compose :: Braid n -> Braid n -> Braid n +compose (Braid gs) (Braid hs) = freeReduceBraidWord $ Braid (gs++hs) + +composeMany :: [Braid n] -> Braid n +composeMany = freeReduceBraidWord . Braid . concat . map braidWord + +-- | Composes two braids without doing any reduction. +composeDontReduce :: Braid n -> Braid n -> Braid n +composeDontReduce (Braid gs) (Braid hs) = Braid (gs++hs) + +-------------------------------------------------------------------------------- +-- * Braid permutations + +-- | A braid is pure if its permutation is trivial +isPureBraid :: KnownNat n => Braid n -> Bool +isPureBraid braid = (braidPermutation braid == P.identity n) where + n = numberOfStrands braid + +-- | Returns the left-to-right permutation associated to the braid. +-- We follow the strands /from the left to the right/ (or from the top to the +-- bottom), and return the permutation taking the left side to the right side. +-- +-- This is compatible with /right/ (standard) action of the permutations: +-- @permuteRight (braidPermutationRight b1)@ corresponds to the left-to-right +-- permutation of the strands; also: +-- +-- > (braidPermutation b1) `multiply` (braidPermutation b2) == braidPermutation (b1 `compose` b2) +-- +-- Writing the right numbering of the strands below the left numbering, +-- we got the two-line notation of the permutation. +-- +braidPermutation :: KnownNat n => Braid n -> Permutation +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 + + action :: forall s. ST s (STUArray s Int Int) + action = do + arr <- newArray_ (1,n) + forM_ [1..n] $ \i -> writeArray arr i i + worker arr idxs + return arr + + worker arr = go where + go [] = return arr + go (i:is) = do + a <- readArray arr i + b <- readArray arr (i+1) + writeArray arr i b + writeArray arr (i+1) a + go is + +-------------------------------------------------------------------------------- +-- * Permutation braids + +-- | A positive braid word contains only positive (@Sigma@) generators. +isPositiveBraidWord :: KnownNat n => Braid n -> Bool +isPositiveBraidWord (Braid gs) = all (isPlus . brGenSign) gs + +-- | A /permutation braid/ is a positive braid where any two strands cross +-- at most one, and /positively/. +-- +isPermutationBraid :: KnownNat n => Braid n -> Bool +isPermutationBraid braid = isPositiveBraidWord braid && crosses where + crosses = and [ check i j | i<-[1..n-1], j<-[i+1..n] ] + check i j = zeroOrOne (lkMatrix ! (i,j)) + zeroOrOne a = (a==1 || a==0) + lkMatrix = linkingMatrix braid + n = numberOfStrands braid + +-- | Untyped version of 'isPermutationBraid' for positive words. +_isPermutationBraid :: Int -> [Int] -> Bool +_isPermutationBraid n gens = crosses where + crosses = and [ check i j | i<-[1..n-1], j<-[i+1..n] ] + check i j = zeroOrOne (lkMatrix ! (i,j)) + zeroOrOne a = (a==1 || a==0) + lkMatrix = _linkingMatrix n $ map Sigma gens + +-- | For any permutation this functions returns a /permutation braid/ realizing +-- that permutation. Note that this is not unique, so we make an arbitrary choice +-- (except for the permutation @[n,n-1..1]@ reversing the order, in which case +-- the result must be the half-twist braid). +-- +-- The resulting braid word will have a length at most @choose n 2@ (and will have +-- that length only for the permutation @[n,n-1..1]@) +-- +-- > braidPermutationRight (permutationBraid perm) == perm +-- > isPermutationBraid (permutationBraid perm) == True +-- +permutationBraid :: KnownNat n => Permutation -> Braid n +permutationBraid perm = braid where + n1 = numberOfStrands braid + n2 = P.permutationSize perm + braid = if n1 == n2 + then Braid (map Sigma $ _permutationBraid perm) + else error $ "permutationBraid: incompatible n: " ++ show n1 ++ " vs. " ++ show n2 + +-- | Untyped version of 'permutationBraid' +_permutationBraid :: Permutation -> [Int] +_permutationBraid = concat . _permutationBraid' + +-- | 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 + + action :: forall s. ST s [[Int]] + action = do + + -- cfwd = the current state of strands : cfwd!j = where is strand #j now? + -- cinv = the inverse of that permutation : cinv!i = which strand is on the #i position now? + + cfwd <- newArray_ (1,n) :: ST s (STUArray s Int Int) + cinv <- newArray_ (1,n) :: ST s (STUArray s Int Int) + forM_ [1..n] $ \j -> do + writeArray cfwd j j + writeArray cinv j j + + let doSwap i = do + a <- readArray cinv i + b <- readArray cinv (i+1) + writeArray cinv i b + writeArray cinv (i+1) a + + u <- readArray cfwd a + v <- readArray cfwd b + writeArray cfwd a v + writeArray cfwd b u + + -- at the k-th phase, we move the (inv!k)-th strand, which is the k-th strand /on the RHS/, to correct position. + let worker phase + | phase >= n = return [] + | otherwise = do + let tgt = (arr ! phase) + src <- readArray cfwd tgt + let this = [src-1,src-2..phase] + mapM_ doSwap $ this + rest <- worker (phase+1) + return (this:rest) + + worker 1 + + +-- | We compute the linking numbers between all pairs of strands: +-- +-- > linkingMatrix braid ! (i,j) == strandLinking braid i j +-- +linkingMatrix :: KnownNat n => Braid n -> UArray (Int,Int) Int +linkingMatrix braid@(Braid gens) = _linkingMatrix (numberOfStrands braid) gens where + +-- | Untyped version of 'linkingMatrix' +_linkingMatrix :: Int -> [BrGen] -> UArray (Int,Int) Int +_linkingMatrix n gens = runSTUArray action where + + action :: forall s. ST s (STUArray s (Int,Int) Int) + action = do + perm <- newArray_ (1,n) :: ST s (STUArray s Int Int) + forM_ [1..n] $ \i -> writeArray perm i i + let doSwap :: Int -> ST s () + doSwap i = do + a <- readArray perm i + b <- readArray perm (i+1) + writeArray perm i b + writeArray perm (i+1) a + + mat <- newArray ((1,1),(n,n)) 0 :: ST s (STUArray s (Int,Int) Int) + let doAdd :: Int -> Int -> Int -> ST s () + doAdd i j pm1 = do + x <- readArray mat (i,j) + writeArray mat (i,j) (x+pm1) + writeArray mat (j,i) (x+pm1) + + forM_ gens $ \g -> do + let (sgn,k) = brGenSignIdx g + u <- readArray perm k + v <- readArray perm (k+1) + doAdd u v (signValue sgn) + doSwap k + + return mat + + +-- | The linking number between two strands numbered @i@ and @j@ +-- (numbered such on the /left/ side). +strandLinking :: KnownNat n => Braid n -> Int -> Int -> Int +strandLinking braid@(Braid gens) i0 j0 + | i0 < 1 || i0 > n = error $ "strandLinkingNumber: invalid strand index i: " ++ show i0 + | j0 < 1 || j0 > n = error $ "strandLinkingNumber: invalid strand index j: " ++ show j0 + | i0 == j0 = 0 + | otherwise = go i0 j0 gens + where + n = numberOfStrands braid + + go !i !j [] = 0 + go !i !j (g:gs) + | i == k && j == k+1 = s + go (i+1) (j-1) gs + | j == k && i == k+1 = s + go (i-1) (j+1) gs + | i == k = go (i+1) j gs + | i == k+1 = go (i-1) j gs + | j == k = go i (j+1) gs + | j == k+1 = go i (j-1) gs + | otherwise = go i j gs + where + (sgn,k) = brGenSignIdx g + s = signValue sgn + + +-------------------------------------------------------------------------------- +-- * ASCII diagram + +instance KnownNat n => DrawASCII (Braid n) where + ascii = horizBraidASCII + +-- | Horizontal braid diagram, drawn from left to right, +-- with strands numbered from the bottom to the top +horizBraidASCII :: KnownNat n => Braid n -> ASCII +horizBraidASCII = horizBraidASCII' True + +-- | Horizontal braid diagram, drawn from left to right. +-- The boolean flag indicates whether to flip the strands +-- vertically ('True' means bottom-to-top, 'False' means top-to-bottom) +horizBraidASCII' :: KnownNat n => Bool -> Braid n -> ASCII +horizBraidASCII' flipped braid@(Braid gens) = final where + + n = numberOfStrands braid + + final = vExtendWith VTop 1 $ hCatTop allBlocks + allBlocks = prelude ++ middleBlocks ++ epilogue + prelude = [ numberBlock , spaceBlock , beginEndBlock ] + epilogue = [ beginEndBlock , spaceBlock , numberBlock' ] + middleBlocks = map block gens + + block g = case g of + Sigma i -> block' i $ if flipped then over else under + SigmaInv i -> block' i $ if flipped then under else over + + block' i middle = asciiFromLines $ drop 2 $ concat + $ replicate a horiz ++ [space3, middle] ++ replicate b horiz + where + (a,b) = if flipped then (n-i-1,i-1) else (i-1,n-i-1) + + -- cycleN :: Int -> [a] -> [a] + -- cycleN n = concat . replicate n + + spaceBlock = transparentBox (1,n*3-2) + beginEndBlock = asciiFromLines $ drop 2 $ concat $ replicate n horiz + numberBlock = mkNumbers [1..n] + numberBlock' = mkNumbers $ P.fromPermutation $ braidPermutation braid + + mkNumbers :: [Int] -> ASCII + mkNumbers list = vCatWith HRight (VSepSpaces 2) $ map asciiShow + $ (if flipped then reverse else id) $ list + + under = [ "\\ /" , " / " , "/ \\" ] + over = [ "\\ /" , " \\ " , "/ \\" ] + horiz = [ " " , " " , "___" ] + space3 = [ " " , " " , " " ] + +-------------------------------------------------------------------------------- + +{- this is unusably ugly and vertically loooong + +-- | Vertical braid diagram, drawn from the top to the bottom. +-- Strands are numbered from the left to the right. +-- +-- Writing down the strand numbers from the top and and the bottom +-- gives the two-line notation of the permutation realized by the braid. +-- +verticalBraidASCII :: KnownNat n => Braid n -> ASCII +verticalBraidASCII braid@(Braid gens) = final where + + n = numberOfStrands braid + + final = hExtendWith HLeft 1 $ vCatLeft allBlocks + allBlocks = prelude ++ middleBlocks ++ epilogue + prelude = [ numberBlock , spaceBlock , beginEndBlock ] + epilogue = [ beginEndBlock , spaceBlock , numberBlock' ] + middleBlocks = map block gens + + block g = case g of + Sigma i -> block' i under + SigmaInv i -> block' i over + + block' i middle = asciiFromLines (map f middle) where + f xs = drop 1 $ concat $ h (i-1) ++ [" ",xs] ++ h (n-i-1) + h k = replicate k " |" + + spaceBlock = transparentBox (n*3-2,1) + beginEndBlock = asciiFromLines $ replicate 3 $ drop 1 $ concat (replicate n " |") + numberBlock = mkNumbers [1..n] + numberBlock' = mkNumbers $ P.fromPermutation $ braidPermutation braid + + mkNumbers :: [Int] -> ASCII + mkNumbers list = asciiFromString (drop 1 $ concatMap show3 list) + show3 k = let s = show k + in replicate (3-length s) ' ' ++ s + + under = [ "\\ /" , " / " , "/ \\" ] + over = [ "\\ /" , " \\ " , "/ \\" ] + +-} + +-------------------------------------------------------------------------------- +-- * Random braids + +-- | Random braid word of the given length +randomBraidWord :: (RandomGen g, KnownNat n) => Int -> g -> (Braid n, g) +randomBraidWord len gen = (braid, gen') where + braid = Braid (map sig bjs) + n = numberOfStrands braid + (gen',bjs) = mapAccumL worker gen [1..len] + + worker !g _ = (g'',(b,j)) where + (j, g' ) = randomR (1,n-1) g + (b, g'') = random g' + + sig :: (Bool,Int) -> BrGen + sig (True ,j) = Sigma j + sig (False,j) = SigmaInv j + +-- | Random /positive/ braid word of the given length +randomPositiveBraidWord :: (RandomGen g, KnownNat n) => Int -> g -> (Braid n, g) +randomPositiveBraidWord len gen = (braid, gen') where + braid = Braid (map Sigma js) + n = numberOfStrands braid + (gen',js) = mapAccumL (\(!g) _ -> swap (randomR (1,n-1) g)) gen [1..len] + +-------------------------------------------------------------------------------- + +-- | Given a braid word, we perturb it randomly @m@ times using the braid relations, +-- so that the resulting new braid word is equivalent to the original. +-- +-- Useful for testing. +-- +randomPerturbBraidWord :: forall n g. (RandomGen g, KnownNat n) => Int -> Braid n -> g -> (Braid n, g) +randomPerturbBraidWord m braid@(Braid xs) g = (Braid word' , g') where + + (word',g') = go m (length xs) xs g + + n = numberOfStrands braid + + -- | A random pair cancelling each other + rndE :: g -> ([BrGen],g) + rndE g = (e1,g'') where + (i , g' ) = randomR (1,n-1) g + (b , g'' ) = random g' + e0 = [SigmaInv i, Sigma i] + e1 = if b then reverse e0 else e0 + + brg s i = case s of { Plus -> Sigma i ; Minus -> SigmaInv i } + brginv s i = case s of { Plus -> SigmaInv i ; Minus -> Sigma i } + + go :: Int -> Int -> [BrGen] -> g -> ([BrGen], g) + go !cnt !len !word !g + + | cnt <= 0 = (word, g) + + | len < 2 = let w' = if b1 then (e++word) else (word++e) -- if it is short, we just add a trivial pair somewhere + in continue g4 (len+2) w' + + | abs (i-j) >= 2 = continue g4 len (as ++ v:u:bs) -- they commute, so we just commute them + + | i == j && s/=t = continue g4 (len-2) (as ++ bs ) -- they are inverse of each other, so we kill them + + | abs (i-j) == 1 && s == t = let mid = if b1 + then [ brg s j , brg s i , brg s j , brginv s i ] -- insert pair and + else [ brginv s j , brg s i , brg s j , brg s i ] -- apply ternary relation + in continue g4 (len+2) (as ++ mid ++ bs) + + | otherwise = let mid = if b1 + then (u : e ++ [v]) + else if b2 + then [u,v] ++ e + else e ++ [u,v] + in continue g4 (len+2) (as++(u:e)++[v]++bs) -- otherwise we just insert an trivial pair + + where + + (pos , g1 ) = randomR (0,len-2) g + (b1 :: Bool , g2 ) = random g1 + (b2 :: Bool , g3 ) = random g2 + (e , g4 ) = rndE g3 + (as,u:v:bs) = splitAt pos word + (s,i) = brGenSignIdx u + (t,j) = brGenSignIdx v + + continue g' len' word' = go (cnt-1) len' word' g' + +-------------------------------------------------------------------------------- + +#ifdef QUICKCHECK + +-- | A permutation braid made convenient to use (type-level hackery) +data PermBraid = forall n. KnownNat n => PermBraid Permutation (Braid n) + +mkPermBraid :: Permutation -> PermBraid +mkPermBraid perm = + case snat of + SomeNat pxy -> PermBraid perm (asProxyTypeOf1 (permutationBraid perm) pxy) + where + n = P.permutationSize perm + Just snat = someNatVal (fromIntegral n :: Integer) + +prop_permBraid_perm :: PermBraid -> Bool +prop_permBraid_perm (PermBraid perm braid) = (braidPermutation braid == perm) + +prop_permBraid_valid :: PermBraid -> Bool +prop_permBraid_valid (PermBraid perm braid) = isPermutationBraid braid + +prop_braidPerm_comp :: KnownNat n => Braid n -> Braid n -> Bool +prop_braidPerm_comp b1 b2 = (p == q) where + p = braidPermutation (compose b1 b2) + q = braidPermutation b1 `P.multiply` braidPermutation b2 + + +#endif + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Groups/Braid/NF.hs b/Math/Combinat/Groups/Braid/NF.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Groups/Braid/NF.hs @@ -0,0 +1,544 @@ + +-- | Normal form of braids, take 1. +-- +-- We implement the Adyan-Thurston-ElRifai-Morton solution to the word problem in braid groups. +-- +-- +-- Based on: +-- +-- * [1] Joan S. Birman, Tara E. Brendle: BRAIDS - A SURVEY +-- (chapter 5.1) +-- +-- * [2] Elsayed A. Elrifai, Hugh R. Morton: Algorithms for positive braids +-- + +{-# LANGUAGE + CPP, BangPatterns, + ScopedTypeVariables, ExistentialQuantification, + DataKinds, KindSignatures, Rank2Types #-} + +module Math.Combinat.Groups.Braid.NF + ( BraidNF (..) + , nfReprWord + , braidNormalForm + , braidNormalForm' +#ifdef QUICKCHECK +#endif + ) + where + +-------------------------------------------------------------------------------- + +import Data.Proxy +import GHC.TypeLits + +import Control.Monad + +import Data.List ( mapAccumL , foldl' , (\\) ) + +import Data.Array.Unboxed +import Data.Array.ST +import Data.Array.IArray +import Data.Array.MArray +import Data.Array.Unsafe +import Data.Array.Base + +import Control.Monad.ST + +import Math.Combinat.Helper +import Math.Combinat.Sign + +import Math.Combinat.Permutations ( Permutation(..) , isIdentityPermutation , isReversePermutation ) +import qualified Math.Combinat.Permutations as P + +import Math.Combinat.Groups.Braid + +#ifdef QUICKCHECK +import Test.QuickCheck +#endif + +-------------------------------------------------------------------------------- + +-- | A unique normal form for braids, called the /left-greedy normal form/. +-- It looks like @Delta^i*P@, where @Delta@ is the positive half-twist, @i@ is an integer, +-- and @P@ is a positive word, which can be further decomposed into non-@Delta@ /permutation words/; +-- these words themselves are not unique, but the permutations they realize /are/ unique. +-- +-- This will solve the word problem relatively fast, +-- though it is not the fastest known algorithm. +-- +data BraidNF (n :: Nat) = BraidNF + { _nfDeltaExp :: !Int -- ^ the exponent of @Delta@ + , _nfPerms :: [Permutation] -- ^ the permutations + } + deriving (Eq,Ord,Show) + +-- | A braid word representing the given normal form +nfReprWord :: KnownNat n => BraidNF n -> Braid n +nfReprWord (BraidNF k perms) = freeReduceBraidWord $ composeMany (deltas ++ rest) where + + deltas + | k > 0 = replicate k halfTwist + | k < 0 = replicate (-k) (inverse halfTwist) + | otherwise = [] + + rest = map permutationBraid perms + +-------------------------------------------------------------------------------- + +-- | Computes the normal form of a braid. We apply free reduction first, it should be faster that way. +braidNormalForm :: KnownNat n => Braid n -> BraidNF n +braidNormalForm = braidNormalForm' . freeReduceBraidWord + +-- | This function does not apply free reduction before computing the normal form +braidNormalForm' :: KnownNat n => Braid n -> BraidNF n +braidNormalForm' braid@(Braid gens) = BraidNF (dexp+pexp) perms where + n = numberOfStrands braid + invless = replaceInverses n gens + -- invless = replaceInversesNaive gens + (dexp,posxword) = moveDeltasLeft n invless + factors = leftGreedyFactors n $ expandPosXWord n posxword + (pexp,perms) = normalizePermFactors n $ map (_braidPermutation n) factors + +-------------------------------------------------------------------------------- + +-- | Replaces groups of @sigma_i^-1@ generators by @(Delta^-1 * P)@, +-- where @P@ is a positive word. +-- +-- This should be more clever (resulting in shorter words) than the naive version below +-- +replaceInverses :: Int -> [BrGen] -> [XGen] +replaceInverses n gens = worker gens where + + worker [] = [] + worker xs = replaceNegs neg ++ map (XSigma . brGenIdx) pos ++ worker rest where + (neg,tmp ) = span (isMinus . brGenSign) xs + (pos,rest) = span (isPlus . brGenSign) tmp + + replaceNegs gs = concatMap replaceFac facs where + 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) + + +-- | Replaces @sigma_i^-1@ generators by @(Delta^-1 * L_i)@. +replaceInversesNaive :: [BrGen] -> [XGen] +replaceInversesNaive gens = concatMap f gens where + f (Sigma i) = [ XSigma i ] + f (SigmaInv i) = [ XDelta (-1) , XL i ] + +-------------------------------------------------------------------------------- + +-- | Temporary data structure to be used during the normal form computation +data XGen + = XDelta !Int -- ^ @Delta^k@ + | XSigma !Int -- ^ @Sigma_j@ + | XL !Int -- ^ @L_j = Delta * sigma_j^-1@ + | XTauL !Int -- ^ @tau(L_j)@ + deriving (Eq,Show) + +isXDelta :: XGen -> Bool +isXDelta x = case x of { XDelta {} -> True ; _ -> False } + +-- | We move the all @Delta@'s to the left +moveDeltasLeft :: Int -> [XGen] -> (Int,[XGen]) +moveDeltasLeft n input = (finalExp, finalPosWord) where + + (XDelta finalExp : finalPosWord) = reverse $ worker 0 (reverse input) + + -- we start from the right end, and work towards the left end + worker dexp [] = [ XDelta dexp ] + worker !dexp xs = this' ++ worker dexp' rest where + (delta,notdelta) = span isXDelta xs + (this ,rest ) = span (not . isXDelta) notdelta + dexp' = dexp + sumDeltas delta + this' = if even dexp' + then this + else map xtau this + + sumDeltas :: [XGen] -> Int + sumDeltas xs = foldl' (+) 0 [ k | XDelta k <- xs ] + + -- | The @X -> Delta^-1 * X * Delta@ inner automorphism + xtau :: XGen -> XGen + xtau (XSigma j) = XSigma (n-j) + xtau (XDelta k) = XDelta k + xtau (XL k) = XTauL k + xtau (XTauL k) = XL k + +-------------------------------------------------------------------------------- + +-- | Expands a /positive/ \"X-word\" into a positive braid word +expandPosXWord :: Int -> [XGen] -> [Int] +expandPosXWord n = concatMap f where + + posHalfTwist = _halfTwist n + + jtau :: Int -> Int + jtau j = n-j + + posLTable = listArray (1,n-1) [ _permutationBraid (posLPerm n i) | i<-[1..n-1] ] :: Array Int [Int] + posTauLTable = amap (map jtau) posLTable + + -- posRTable = listArray (1,n-1) [ _permutationBraid (posRPerm n i) | i<-[1..n-1] ] :: Array Int [Int] + + f x = case x of + XSigma i -> [i] + XL i -> posLTable ! i + XTauL i -> posTauLTable ! i + XDelta i + | i > 0 -> concat (replicate i posHalfTwist) + | i < 0 -> error "expandPosXWord: negative delta power" + | otherwise -> [] + + -- word :: Braid n -> [Int] + -- word (Braid gens) = map brGenIdx gens + + +-- | Expands an \"X-word\" into a braid word. Useful for debugging. +expandAnyXWord :: forall n. KnownNat n => [XGen] -> Braid n +expandAnyXWord xgens = braid where + n = numberOfStrands braid + + braid = composeMany (map f xgens) + + posHalfTwist = halfTwist :: Braid n + negHalfTwist = inverse posHalfTwist :: Braid n + + posLTable = listArray (1,n-1) [ permutationBraid (posLPerm n i) | i<-[1..n-1] ] :: Array Int (Braid n) + posTauLTable = amap tau posLTable + + -- posRTable = listArray (1,n-1) [ permutationBraid (posRPerm n i) | i<-[1..n-1] ] :: Array Int (Braid n) + + f :: XGen -> Braid n + f x = case x of + XSigma i -> sigma i + XL i -> posLTable ! i + XTauL i -> posTauLTable ! i + XDelta i + | i > 0 -> composeMany (replicate i posHalfTwist) + | i < 0 -> composeMany (replicate (-i) negHalfTwist) + | otherwise -> identity + +-------------------------------------------------------------------------------- + +-- | @posL k@ (denoted as @L_k@) is a /positive word/ which +-- satisfies @Delta = L_k * sigma_k@, or: +-- +-- > (inverse halfTwist) `compose` (posL k) ~=~ sigmaInv k@ +-- +-- Thus we can replace any word with a positive word plus some @Delta^-1@\'s +-- +posL :: KnownNat n => Int -> Braid n +posL k = braid where + n = numberOfStrands braid + braid = permutationBraid (posLPerm n k) + +-- | @posR k n@ (denoted as @R_k@) is a /permutation braid/ which +-- satisfies @Delta = sigma_k * R_k@ +-- +-- > (posR k) `compose` (inverse halfTwist) ~=~ sigmaInv k@ +-- +-- Thus we can replace any word with a positive word plus some @Delta^-1@'s +-- +posR :: KnownNat n => Int -> Braid n +posR k = braid where + n = numberOfStrands braid + braid = permutationBraid (posRPerm n k) + +-- | The permutation @posL k :: Braid n@ is realizing +posLPerm :: Int -> Int -> Permutation +posLPerm n k + | k>0 && k Int -> Permutation +posRPerm n k + | k>0 && k [[Int]] -> (Int, [[Int]]) +filterDeltaFactors n facs = (exp',facs'') where + + (exp',facs') = go 0 (reverse facs) + + jtau j = n-j + facs'' = reverse facs' + maxlen = div (n*(n-1)) 2 + + go !e [] = (e,[]) + go !e (xs:xxs) + | null xs = go e xxs + | length xs == maxlen = go (e+1) xxs + | otherwise = + if even e + then let (e',yys) = go e xxs in (e' , xs : yys) + else let (e',yys) = go e xxs in (e' , map jtau xs : yys) + +-------------------------------------------------------------------------------- + +-- | The /starting set/ of a positive braid P is the subset of @[1..n-1]@ defined by +-- +-- > S(P) = [ i | P = sigma_i * Q , Q is positive ] = [ i | (sigma_i^-1 * P) is positive ] +-- +-- This function returns the starting set a positive word, assuming it +-- is a /permutation braid/ (see Lemma 2.4 in [2]) +-- +permWordStartingSet :: Int -> [Int] -> [Int] +permWordStartingSet n xs = permWordFinishingSet n (reverse xs) + +-- | The /finishing set/ of a positive braid P is the subset of @[1..n-1]@ defined by +-- +-- > F(P) = [ i | P = Q * sigma_i , Q is positive ] = [ i | (P * sigma_i^-1) is positive ] +-- +-- This function returns the finishing set, assuming the input is a /permutation braid/ +-- +permWordFinishingSet :: Int -> [Int] -> [Int] +permWordFinishingSet n input = runST action where + + action :: forall s. ST s [Int] + action = do + perm <- newArray_ (1,n) :: ST s (STUArray s Int Int) + forM_ [1..n] $ \i -> writeArray perm i i + forM_ input $ \i -> do + a <- readArray perm i + b <- readArray perm (i+1) + writeArray perm i b + writeArray perm (i+1) a + flip filterM [1..n-1] $ \i -> do + a <- readArray perm i + b <- readArray perm (i+1) + return (b permutationStartingSet p == permWordStartingSet n (_permutationBraid p) +-- +permutationStartingSet :: Permutation -> [Int] +permutationStartingSet = permutationFinishingSet . P.inverse + +-- | 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 + +-- | Returns the list of permutations failing Lemma 2.5 in [2] +-- (so an empty list means the implementaton is correct) +fails_lemmma_2_5 :: Int -> [Permutation] +fails_lemmma_2_5 n = [ p | p <- P.permutations n , not (test p) ] where + test p = and [ check i | i<-[1..n-1] ] where + w = _permutationBraid p + s = permWordStartingSet n w + check i = _isPermutationBraid n (i:w) == (not $ elem i s) + +-------------------------------------------------------------------------------- + +-- | Given factors defined as permutation braids, we normalize them +-- to /left-canonical form/ by ensuring that +-- +-- * for each consecutive pair @(P,Q)@ the finishing set F(P) contains the starting set S(Q) +-- +-- * all @Delta@-s (corresponding to the reverse permutation) are moved to the left +-- +-- * all trivial factors are filtered out +-- +-- Unfortunately, it seems that we may need multiple sweeps to do that... +-- +normalizePermFactors :: Int -> [Permutation] -> (Int,[Permutation]) +normalizePermFactors n = go 0 where + go !acc input = + if (exp==0 && input == output) + then (acc,input) + else go (acc+exp) output + where + (exp,output) = normalizePermFactors1 n input + +-- | Does 1 sweep of the above normalization process. +-- Unfortunately, it seems that we may need to do this multiple times... +-- +normalizePermFactors1 :: Int -> [Permutation] -> (Int,[Permutation]) +normalizePermFactors1 n input = (exp, reverse output) where + (exp, output) = worker 0 (reverse input) + + -- Notes: We work in reverse order, from the right to the left. + -- We maintain the number of Delta-s pushed through; the tau involutions + -- are implicit in the parity of this number + -- + worker :: Int -> [Permutation] -> (Int,[Permutation]) + worker = worker' 0 0 + + -- We also maintain additional 0/1 flip flags for the first two permutations + -- this is a little bit of hack but it should work nicely + -- + worker' :: Int -> Int -> Int -> [Permutation] -> (Int,[Permutation]) + worker' !ep !eq !e (!p : rest@(!q : rest')) + + -- check if the very first element is identity or Delta + -- (note: these are tau-invariants) + + | isIdentityPermutation p = worker' eq 0 e rest + | isReversePermutation p = worker' eq 0 (e+1) rest + + -- check if the second element is identity or Delta + -- this is necessary since we "fatten" the second element and it can possibly + -- become Delta after a while (?) + + | isIdentityPermutation q = worker' ep 0 e (p : rest') + | isReversePermutation q = worker' (ep-1) 0 (e+1) (p : rest') + + -- ok so we have something like "... : Q : P" + -- if F(Q) contains S(P) then we can move on; + -- otherwise there is an element j in S(P) \\ F(Q), so we can + -- replace it by "... : Qj : jP" + + | otherwise = + case permutationStartingSet preal \\ permutationFinishingSet qreal of + [] -> 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 + where + preal = oddTau (e+ep) p -- the "real" p + qreal = oddTau (e+eq) q -- the "real" q + + worker' _ _ !e [ ] = (e,[]) + worker' !ep _ !e [p] + | isIdentityPermutation p = (e , []) + | isReversePermutation p = (e+1 , []) + | otherwise = (e , [oddTau (e+ep) p] ) + + oddTau :: Int -> Permutation -> Permutation + oddTau !e p = if even e then p else permTau p + +{- + checkDelta :: Int -> Permutation -> [Permutation] -> (Int,[Permutation]) + checkDelta !e !p !rest + | P.isIdentityPermutation p = worker e rest + | isReversePermutation p = worker (e+1) rest + | otherwise = let (e',rs) = worker e rest in (e', oddTau e p : rs) +-} + +-- | The involution tau on permutation +permTau :: Permutation -> Permutation +permTau (Permutation arr) = Permutation $ listArray (1,n) [ (n+1) - arr!(n-i) | i<-[0..n-1] ] where + (1,n) = bounds arr + +-------------------------------------------------------------------------------- + +-- | Given a /positive/ word, we apply left-greedy factorization of +-- that word into subwords representing /permutation braids/. +-- +-- Example 5.1 from the above handbook: +-- +-- > leftGreedyFactors 7 [1,3,2,2,1,3,3,2,3,2] == [[1,3,2],[2,1,3],[3,2,3],[2]] +-- +leftGreedyFactors :: Int -> [Int] -> [[Int]] +leftGreedyFactors n input = filter (not . null) $ runST (action input) where + + action :: forall s. [Int] -> ST s [[Int]] + action input = do + + perm <- newArray_ (1,n) :: ST s (STUArray s Int Int) + forM_ [1..n] $ \i -> writeArray perm i i + let doSwap :: Int -> ST s () + doSwap i = do + a <- readArray perm i + b <- readArray perm (i+1) + writeArray perm i b + writeArray perm (i+1) a + + mat <- newArray ((1,1),(n,n)) 0 :: ST s (STUArray s (Int,Int) Int) + let clearMat = forM_ [1..n] $ \i -> + forM_ [1..n] $ \j -> writeArray mat (i,j) 0 + + let doAdd1 :: Int -> Int -> ST s Int + doAdd1 i j = do + x <- readArray mat (i,j) + let y = x+1 + writeArray mat (i,j) y + writeArray mat (j,i) y + return y + + let worker :: [Int] -> ST s [[Int]] + worker [] = return [[]] + worker (p:ps) = do + u <- readArray perm p + v <- readArray perm (p+1) + c <- doAdd1 u v + doSwap p + if c<=1 + then do + (f:fs) <- worker ps + return ((p:f):fs) + else do + clearMat + fs <- worker (p:ps) + return ([]:fs) + + worker input + +-------------------------------------------------------------------------------- + +{- + +-- | Finds ternary braid relations, and returns them as a list of indices, decorated +-- with a flag specifying which side of the relation we found, a sign specifying +-- whether it is a relation between positive or negative generators. +-- +findTernaryBraidRelations :: Braid n -> [(Int,Bool,Sign)] +findTernaryBraidRelations (Braid gens) = go 0 gens where + go !k (Sigma a : rest@(Sigma b : Sigma c : _)) + | a==c && b==a+1 = (k,True ,Plus) : go (k+1) rest + | a==c && b==a-1 = (k,False,Plus) : go (k+1) rest + | otherwise = go (k+1) rest + go !k (SigmaInv a : rest@(SigmaInv b : SigmaInv c : _)) + | a==c && b==a+1 = (k,True ,Minus) : go (k+1) rest + | a==c && b==a-1 = (k,False,Minus) : go (k+1) rest + | otherwise = go (k+1) rest + go !k (x:xs) = go (k+1) xs + go _ [] = [] + +-- | Finds subsequences like @(i,i+1,i)@ and @(i+1,i,i+1)@, and returns them +-- and a list of indices, plus a flag specifying which one we found (the first +-- one is 'True', second one is 'False') +-- +_findTernaryBraidRelations :: [Int] -> [(Int,Bool)] +_findTernaryBraidRelations = go 0 where + go !k (a:rest@(b:c:_)) + | a==c && b==a+1 = (k,True ) : go (k+1) rest + | a==c && b==a-1 = (k,False) : go (k+1) rest + | otherwise = go (k+1) rest + go !k (x:xs) = go (k+1) xs + go _ [] = [] + +-} + +-------------------------------------------------------------------------------- + +#ifdef QUICKCHECK + +prop_braidnf_reduce :: KnownNat n => Braid n -> Bool +prop_braidnf_reduce braid = (braidNormalForm' braid == braidNormalForm braid) + +prop_braidnf_reprs :: KnownNat n => Braid n -> Bool +prop_braidnf_reprs braid = (nf == nf') where + nf = braidNormalForm braid + nf' = braidNormalForm braid' + braid' = nfReprWord nf + +#endif + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Groups/Free.hs b/Math/Combinat/Groups/Free.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Groups/Free.hs @@ -0,0 +1,523 @@ + +-- | Words in free groups (and free powers of cyclic groups). +-- +-- This module is not re-exported by "Math.Combinat" +-- +{-# LANGUAGE CPP, BangPatterns, PatternGuards #-} +module Math.Combinat.Groups.Free where + +-------------------------------------------------------------------------------- + +-- new Base exports "Word" from Data.Word... +#ifdef MIN_VERSION_base +#if MIN_VERSION_base(4,7,1) +import Prelude hiding ( Word ) +#endif +#elif __GLASGOW_HASKELL__ >= 709 +import Prelude hiding ( Word ) +#endif + +import Data.Char ( chr ) +import Data.List ( mapAccumL , groupBy ) + +import Control.Monad ( liftM ) +import System.Random + +import Math.Combinat.Numbers +import Math.Combinat.Sign +import Math.Combinat.Helper + +-------------------------------------------------------------------------------- +-- * Words + +-- | A generator of a (free) group, indexed by which \"copy\" of the group we are dealing with. +data Generator idx + = Gen !idx -- @a@ + | Inv !idx -- @a^(-1)@ + deriving (Eq,Ord,Show,Read) + +-- | The index of a generator +genIdx :: Generator idx -> idx +genIdx g = case g of + Gen x -> x + Inv x -> x + +-- | The sign of the (exponent of the) generator (that is, the generator is 'Plus', the inverse is 'Minus') +genSign :: Generator idx -> Sign +genSign g = case g of { Gen _ -> Plus ; Inv _ -> Minus } + +genSignValue :: Generator idx -> Int +genSignValue g = case g of { Gen _ -> (1::Int) ; Inv _ -> (-1::Int) } + +-- | keep the index, but return always the 'Gen' one. +absGen :: Generator idx -> Generator idx +absGen g = case g of + Gen x -> Gen x + Inv x -> Gen x + +-- | A /word/, describing (non-uniquely) an element of a group. +-- The identity element is represented (among others) by the empty word. +type Word idx = [Generator idx] + +-------------------------------------------------------------------------------- + +-- | Generators are shown as small letters: @a@, @b@, @c@, ... +-- and their inverses are shown as capital letters, so @A=a^-1@, @B=b^-1@, etc. +showGen :: Generator Int -> Char +showGen (Gen i) = chr (96+i) +showGen (Inv i) = chr (64+i) + +showWord :: Word Int -> String +showWord = map showGen + +-------------------------------------------------------------------------------- + +instance Functor Generator where + fmap f g = case g of + Gen x -> Gen (f x) + Inv y -> Inv (f y) + +-------------------------------------------------------------------------------- + +-- | The inverse of a generator +inverseGen :: Generator a -> Generator a +inverseGen g = case g of + Gen x -> Inv x + Inv x -> Gen x + +-- | The inverse of a word +inverseWord :: Word a -> Word a +inverseWord = map inverseGen . reverse + +-- | Lists all words of the given length (total number will be @(2g)^n@). +-- The numbering of the generators is @[1..g]@. +allWords + :: Int -- ^ @g@ = number of generators + -> Int -- ^ @n@ = length of the word + -> [Word Int] +allWords g = go where + go !0 = [[]] + go !n = [ x:xs | xs <- go (n-1) , x <- elems ] + elems = [ Gen a | a<-[1..g] ] + ++ [ Inv a | a<-[1..g] ] + +-- | Lists all words of the given length which do not contain inverse generators +-- (total number will be @g^n@). +-- The numbering of the generators is @[1..g]@. +allWordsNoInv + :: Int -- ^ @g@ = number of generators + -> Int -- ^ @n@ = length of the word + -> [Word Int] +allWordsNoInv g = go where + go !0 = [[]] + go !n = [ x:xs | xs <- go (n-1) , x <- elems ] + elems = [ Gen a | a<-[1..g] ] + +-------------------------------------------------------------------------------- +-- * Random words + +-- | A random group generator (or its inverse) between @1@ and @g@ +randomGenerator + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> g -> (Generator Int, g) +randomGenerator !d !g0 = (gen, g2) where + (b, !g1) = random g0 + (k, !g2) = randomR (1,d) g1 + gen = if b then Gen k else Inv k + +-- | A random group generator (but never its inverse) between @1@ and @g@ +randomGeneratorNoInv + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> g -> (Generator Int, g) +randomGeneratorNoInv !d !g0 = (Gen k, g1) where + (!k, !g1) = randomR (1,d) g0 + +-- | A random word of length @n@ using @g@ generators (or their inverses) +randomWord + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> Int -- ^ @n@ = length of the word + -> g -> (Word Int, g) +randomWord !d !n !g0 = (word,g1) where + (g1,word) = mapAccumL (\g _ -> swap (randomGenerator d g)) g0 [1..n] + +-- | A random word of length @n@ using @g@ generators (but not their inverses) +randomWordNoInv + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> Int -- ^ @n@ = length of the word + -> g -> (Word Int, g) +randomWordNoInv !d !n !g0 = (word,g1) where + (g1,word) = mapAccumL (\g _ -> swap (randomGeneratorNoInv d g)) g0 [1..n] + +-------------------------------------------------------------------------------- +-- * The free group on @g@ generators + +{-# SPECIALIZE multiplyFree :: Word Int -> Word Int -> Word Int #-} +{-# SPECIALIZE equivalentFree :: Word Int -> Word Int -> Bool #-} +{-# SPECIALIZE reduceWordFree :: Word Int -> Word Int #-} +{-# SPECIALIZE reduceWordFreeNaive :: Word Int -> Word Int #-} + +-- | Multiplication of the free group (returns the reduced result). It is true +-- for any two words w1 and w2 that +-- +-- > multiplyFree (reduceWordFree w1) (reduceWord w2) = multiplyFree w1 w2 +-- +multiplyFree :: Eq idx => Word idx -> Word idx -> Word idx +multiplyFree w1 w2 = reduceWordFree (w1 ++ w2) + +-- | Decides whether two words represent the same group element in the free group +equivalentFree :: Eq idx => Word idx -> Word idx -> Bool +equivalentFree w1 w2 = null $ reduceWordFree $ w1 ++ inverseWord w2 + +-- | Reduces a word in a free group by repeatedly removing @x*x^(-1)@ and +-- @x^(-1)*x@ pairs. The set of /reduced words/ forms the free group; the +-- multiplication is obtained by concatenation followed by reduction. +-- +reduceWordFree :: Eq idx => Word idx -> Word idx +reduceWordFree = loop where + + loop w = case reduceStep w of + Nothing -> w + Just w' -> loop w' + + reduceStep :: Eq a => Word a -> Maybe (Word a) + reduceStep = go False where + go !changed w = case w of + (Gen x : Inv y : rest) | x==y -> go True rest + (Inv x : Gen y : rest) | x==y -> go True rest + (this : rest) -> liftM (this:) $ go changed rest + _ -> if changed then Just w else Nothing + + +-- | Naive (but canonical) reduction algorithm for the free groups +reduceWordFreeNaive :: Eq idx => Word idx -> Word idx +reduceWordFreeNaive = loop where + loop w = let w' = step w in if w/=w' then loop w' else w + step = concatMap worker . groupBy (equating genIdx) where + worker gs + | s>0 = replicate s (Gen i) + | s<0 = replicate (abs s) (Inv i) + | otherwise = [] + where + i = genIdx (head gs) + s = sum' (map genSignValue gs) + +-------------------------------------------------------------------------------- + +-- | Counts the number of words of length @n@ which reduce to the identity element. +-- +-- Generating function is @Gf_g(u) = \\frac {2g-1} { g-1 + g \\sqrt{ 1 - (8g-4)u^2 } }@ +-- +countIdentityWordsFree + :: Int -- ^ g = number of generators in the free group + -> Int -- ^ n = length of the unreduced word + -> Integer +countIdentityWordsFree g n = countWordReductionsFree g n 0 + +-- | Counts the number of words of length @n@ whose reduced form has length @k@ +-- (clearly @n@ and @k@ must have the same parity for this to be nonzero): +-- +-- > countWordReductionsFree g n k == sum [ 1 | w <- allWords g n, k == length (reduceWordFree w) ] +-- +countWordReductionsFree + :: Int -- ^ g = number of generators in the free group + -> Int -- ^ n = length of the unreduced word + -> Int -- ^ k = length of the reduced word + -> Integer +countWordReductionsFree gens_ nn_ kk_ + | nn==0 = if k==0 then 1 else 0 + | even nn && kk == 0 = sum [ ( binomial (nn-i) (n -i) * gg^(i ) * (gg-1)^(n -i ) * ( i) ) `div` (nn-i) | i<-[0..n ] ] + | even nn && even kk = sum [ ( binomial (nn-i) (n-k-i) * gg^(i+1) * (gg-1)^(n+k-i-1) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] + | odd nn && odd kk = sum [ ( binomial (nn-i) (n-k-i) * gg^(i+1) * (gg-1)^(n+k-i ) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] + | otherwise = 0 + where + g = fromIntegral gens_ :: Integer + nn = fromIntegral nn_ :: Integer + kk = fromIntegral kk_ :: Integer + + gg = 2*g + n = div nn 2 + k = div kk 2 + +-------------------------------------------------------------------------------- +-- * Free powers of cyclic groups + +{-# SPECIALIZE multiplyZ2 :: Word Int -> Word Int -> Word Int #-} +{-# SPECIALIZE multiplyZ3 :: Word Int -> Word Int -> Word Int #-} +{-# SPECIALIZE multiplyZm :: Int -> Word Int -> Word Int -> Word Int #-} + +-- | Multiplication in free products of Z2's +multiplyZ2 :: Eq idx => Word idx -> Word idx -> Word idx +multiplyZ2 w1 w2 = reduceWordZ2 (w1 ++ w2) + +-- | Multiplication in free products of Z3's +multiplyZ3 :: Eq idx => Word idx -> Word idx -> Word idx +multiplyZ3 w1 w2 = reduceWordZ3 (w1 ++ w2) + +-- | Multiplication in free products of Zm's +multiplyZm :: Eq idx => Int -> Word idx -> Word idx -> Word idx +multiplyZm k w1 w2 = reduceWordZm k (w1 ++ w2) + +-------------------------------------------------------------------------------- + +{-# SPECIALIZE equivalentZ2 :: Word Int -> Word Int -> Bool #-} +{-# SPECIALIZE equivalentZ3 :: Word Int -> Word Int -> Bool #-} +{-# SPECIALIZE equivalentZm :: Int -> Word Int -> Word Int -> Bool #-} + +-- | Decides whether two words represent the same group element in free products of Z2 +equivalentZ2 :: Eq idx => Word idx -> Word idx -> Bool +equivalentZ2 w1 w2 = null $ reduceWordZ2 $ w1 ++ inverseWord w2 + +-- | Decides whether two words represent the same group element in free products of Z3 +equivalentZ3 :: Eq idx => Word idx -> Word idx -> Bool +equivalentZ3 w1 w2 = null $ reduceWordZ3 $ w1 ++ inverseWord w2 + +-- | Decides whether two words represent the same group element in free products of Zm +equivalentZm :: Eq idx => Int -> Word idx -> Word idx -> Bool +equivalentZm m w1 w2 = null $ reduceWordZm m $ w1 ++ inverseWord w2 + +-------------------------------------------------------------------------------- + +{-# SPECIALIZE reduceWordZ2 :: Word Int -> Word Int #-} +{-# SPECIALIZE reduceWordZ3 :: Word Int -> Word Int #-} +{-# SPECIALIZE reduceWordZm :: Int -> Word Int -> Word Int #-} + +-------------------------------------------------------------------------------- + +-- | Reduces a word, where each generator @x@ satisfies the additional relation @x^2=1@ +-- (that is, free products of Z2's) +reduceWordZ2 :: Eq idx => Word idx -> Word idx +reduceWordZ2 = loop where + loop w = case reduceStep w of + Nothing -> w + Just w' -> loop w' + + reduceStep :: Eq a => Word a -> Maybe (Word a) + reduceStep = go False where + go !changed w = case w of + (Gen x : Gen y : rest) | x==y -> go True rest + (Gen x : Inv y : rest) | x==y -> go True rest + (Inv x : Gen y : rest) | x==y -> go True rest + (Inv x : Inv y : rest) | x==y -> go True rest + (this : rest) -> liftM (absGen this:) $ go changed rest + _ -> if changed then Just w else Nothing + +-- | Reduces a word, where each generator @x@ satisfies the additional relation @x^3=1@ +-- (that is, free products of Z3's) +reduceWordZ3 :: Eq idx => Word idx -> Word idx +reduceWordZ3 = loop where + loop w = case reduceStep w of + Nothing -> w + Just w' -> loop w' + + reduceStep :: Eq a => Word a -> Maybe (Word a) + reduceStep = go False where + go !changed w = case w of + (Gen x : Inv y : rest) | x==y -> go True rest + (Inv x : Gen y : rest) | x==y -> go True rest + (Gen x : Gen y : Gen z : rest) | x==y && y==z -> go True rest + (Inv x : Inv y : Inv z : rest) | x==y && y==z -> go True rest + (Gen x : Gen y : rest) | x==y -> go True (Inv x : rest) -- !!! + (Inv x : Inv y : rest) | x==y -> go True (Gen x : rest) + (this : rest) -> liftM (this:) $ go changed rest + _ -> if changed then Just w else Nothing + +-- | Reduces a word, where each generator @x@ satisfies the additional relation @x^m=1@ +-- (that is, free products of Zm's) +reduceWordZm :: Eq idx => Int -> Word idx -> Word idx +reduceWordZm m = loop where + + loop w = case reduceStep w of + Nothing -> w + Just w' -> loop w' + + halfm = div m 2 -- if we encounter strictly more than m/2 equal elements in a row, we replace them by the inverses + + -- reduceStep :: Eq a => Word a -> Maybe (Word a) + reduceStep = go False where + go !changed w = case w of + (Gen x : Inv y : rest) | x==y -> go True rest + (Inv x : Gen y : rest) | x==y -> go True rest + something | Just (k,rest) <- dropIfMoreThanHalf w -> go True (replicate (m-k) (inverseGen (head w)) ++ rest) + (this : rest) -> liftM (this:) $ go changed rest + _ -> if changed then Just w else Nothing + + -- dropIfMoreThanHalf :: Eq a => Word a -> Maybe (Int, Word a) + dropIfMoreThanHalf w = + let (!k,rest) = dropWhileEqual w + in if k > halfm then Just (k,rest) + else Nothing + + -- dropWhileEqual :: Eq a => Word a -> (Int, Word a) + dropWhileEqual [] = (0,[]) + dropWhileEqual (x0:rest) = go 1 rest where + go !k [] = (k,[]) + go !k xxs@(x:xs) = if k==m then (m,xxs) + else if x==x0 then go (k+1) xs + else (k,xxs) + +{- + dropm :: Eq a => Word a -> Maybe (Word a) + dropm [] = Nothing + dropm (x:xs) = go (m-1) xs where + go 0 rest = Just rest + go j (y:ys) = if y==x + then go (j-1) ys + else Nothing + go j [] = Nothing +-} + +-------------------------------------------------------------------------------- + +{-# SPECIALIZE reduceWordZ2Naive :: Word Int -> Word Int #-} +{-# SPECIALIZE reduceWordZ3Naive :: Word Int -> Word Int #-} +{-# SPECIALIZE reduceWordZmNaive :: Int -> Word Int -> Word Int #-} + +-- | Reduces a word, where each generator @x@ satisfies the additional relation @x^2=1@ +-- (that is, free products of Z2's). Naive (but canonical) algorithm. +reduceWordZ2Naive :: Eq idx => Word idx -> Word idx +reduceWordZ2Naive = loop where + loop w = let w' = step w in if w/=w' then loop w' else w + step = concatMap worker . groupBy (equating genIdx) where + worker gs = + case mod s 2 of + 1 -> [Gen i] + 0 -> [] + _ -> error "reduceWordZ2: fatal error, shouldn't happen" + where + i = genIdx (head gs) + s = sum' (map genSignValue gs) + +-- | Reduces a word, where each generator @x@ satisfies the additional relation @x^3=1@ +-- (that is, free products of Z3's). Naive (but canonical) algorithm. +reduceWordZ3Naive :: Eq idx => Word idx -> Word idx +reduceWordZ3Naive = loop where + loop w = let w' = step w in if w/=w' then loop w' else w + step = concatMap worker . groupBy (equating genIdx) where + worker gs = + case mod s 3 of + 0 -> [] + 1 -> [Gen i] + 2 -> [Inv i] + _ -> error "reduceWordZ3: fatal error, shouldn't happen" + where + i = genIdx (head gs) + s = sum' (map genSignValue gs) + +-- | Reduces a word, where each generator @x@ satisfies the additional relation @x^m=1@ +-- (that is, free products of Zm's). Naive (but canonical) algorithm. +reduceWordZmNaive :: Eq idx => Int -> Word idx -> Word idx +reduceWordZmNaive m = loop where + loop w = let w' = step w in if w/=w' then loop w' else w + step = concatMap worker . groupBy (equating genIdx) where + halfm1 = div (m+1) 2 + worker gs + | mods <= halfm1 = replicate mods (Gen i) + | otherwise = replicate (m-mods) (Inv i) + where + i = genIdx (head gs) + s = sum' (map genSignValue gs) + mods = mod s m + +-------------------------------------------------------------------------------- + +-- | Counts the number of words (without inverse generators) of length @n@ +-- which reduce to the identity element, using the relations @x^2=1@. +-- +-- Generating function is @Gf_g(u) = \\frac {2g-2} { g-2 + g \\sqrt{ 1 - (4g-4)u^2 } }@ +-- +-- The first few @g@ cases: +-- +-- > A000984 = [ countIdentityWordsZ2 2 (2*n) | n<-[0..] ] = [1,2,6,20,70,252,924,3432,12870,48620,184756...] +-- > A089022 = [ countIdentityWordsZ2 3 (2*n) | n<-[0..] ] = [1,3,15,87,543,3543,23823,163719,1143999,8099511,57959535...] +-- > A035610 = [ countIdentityWordsZ2 4 (2*n) | n<-[0..] ] = [1,4,28,232,2092,19864,195352,1970896,20275660,211823800,2240795848...] +-- > A130976 = [ countIdentityWordsZ2 5 (2*n) | n<-[0..] ] = [1,5,45,485,5725,71445,925965,12335685,167817405,2321105525,32536755565...] +-- +countIdentityWordsZ2 + :: Int -- ^ g = number of generators in the free group + -> Int -- ^ n = length of the unreduced word + -> Integer +countIdentityWordsZ2 g n = countWordReductionsZ2 g n 0 + +-- | Counts the number of words (without inverse generators) of length @n@ whose +-- reduced form in the product of Z2-s (that is, for each generator @x@ we have @x^2=1@) +-- has length @k@ +-- (clearly @n@ and @k@ must have the same parity for this to be nonzero): +-- +-- > countWordReductionsZ2 g n k == sum [ 1 | w <- allWordsNoInv g n, k == length (reduceWordZ2 w) ] +-- +countWordReductionsZ2 + :: Int -- ^ g = number of generators in the free group + -> Int -- ^ n = length of the unreduced word + -> Int -- ^ k = length of the reduced word + -> Integer +countWordReductionsZ2 gens_ nn_ kk_ + | nn==0 = if k==0 then 1 else 0 + | even nn && kk == 0 = sum [ ( binomial (nn-i) (n -i) * g^(i ) * (g-1)^(n -i ) * ( i) ) `div` (nn-i) | i<-[0..n ] ] + | even nn && even kk = sum [ ( binomial (nn-i) (n-k-i) * g^(i+1) * (g-1)^(n+k-i-1) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] + | odd nn && odd kk = sum [ ( binomial (nn-i) (n-k-i) * g^(i+1) * (g-1)^(n+k-i ) * (kk+i) ) `div` (nn-i) | i<-[0..n-k] ] + | otherwise = 0 + where + g = fromIntegral gens_ :: Integer + nn = fromIntegral nn_ :: Integer + kk = fromIntegral kk_ :: Integer + + n = div nn 2 + k = div kk 2 + +-- | Counts the number of words (without inverse generators) of length @n@ +-- which reduce to the identity element, using the relations @x^3=1@. +-- +-- > countIdentityWordsZ3NoInv g n == sum [ 1 | w <- allWordsNoInv g n, 0 == length (reduceWordZ2 w) ] +-- +-- In mathematica, the formula is: @Sum[ g^k * (g-1)^(n-k) * k/n * Binomial[3*n-k-1, n-k] , {k, 1,n} ]@ +-- +countIdentityWordsZ3NoInv + :: Int -- ^ g = number of generators in the free group + -> Int -- ^ n = length of the unreduced word + -> Integer +countIdentityWordsZ3NoInv gens_ nn_ + | nn==0 = 1 + | mod nn 3 == 0 = sum [ ( binomial (3*n-i-1) (n-i) * g^i * (g-1)^(n-i) * i ) `div` n | i<-[1..n] ] + | otherwise = 0 + where + g = fromIntegral gens_ :: Integer + nn = fromIntegral nn_ :: Integer + + n = div nn 3 + +-------------------------------------------------------------------------------- + +{- + +-- some basic testing. TODO: QuickCheck tests + +import Math.Combinat.Helper +import Math.Combinat.Groups.Free + +g = 3 :: Int +maxn = 8 :: Int + +bad_free = [ w | n<-[0..maxn] , w <- allWords g n , not (reduceWordFree w `equivalentFree` reduceWordFreeNaive w) ] +bad_z2 = [ w | n<-[0..maxn] , w <- allWords g n , not (reduceWordZ2 w `equivalentZ2` reduceWordZ2Naive w) ] +bad_z3 = [ w | n<-[0..maxn] , w <- allWords g n , not (reduceWordZ3 w `equivalentZ3` reduceWordZ3Naive w) ] +bad_zm m = [ w | n<-[0..maxn] , w <- allWords g n , not (equivalentZm m (reduceWordZm m w) (reduceWordZmNaive m w)) ] + +speed_free = sum' [ length (reduceWordFree w) | n<-[0..maxn] , w <- allWords g n ] +speed_z2 = sum' [ length (reduceWordZ2 w) | n<-[0..maxn] , w <- allWords g n ] +speed_z3 = sum' [ length (reduceWordZ3 w) | n<-[0..maxn] , w <- allWords g n ] +speed_zm m = sum' [ length (reduceWordZm m w) | n<-[0..maxn] , w <- allWords g n ] + +naive_speed_free = sum' [ length (reduceWordFreeNaive w) | n<-[0..maxn] , w <- allWords g n ] +naive_speed_z2 = sum' [ length (reduceWordZ2Naive w) | n<-[0..maxn] , w <- allWords g n ] +naive_speed_z3 = sum' [ length (reduceWordZ3Naive w) | n<-[0..maxn] , w <- allWords g n ] +naive_speed_zm m = sum' [ length (reduceWordZmNaive m w) | n<-[0..maxn] , w <- allWords g n ] + +-} + +-------------------------------------------------------------------------------- + + diff --git a/Math/Combinat/Groups/Thompson/F.hs b/Math/Combinat/Groups/Thompson/F.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Groups/Thompson/F.hs @@ -0,0 +1,404 @@ + +-- | Thompson's group F. +-- +-- See eg. +-- +-- Based mainly on James Michael Belk's PhD thesis \"THOMPSON'S GROUP F\"; +-- see +-- + +{-# LANGUAGE TypeSynonymInstances, FlexibleInstances, BangPatterns, PatternSynonyms, DeriveFunctor #-} +module Math.Combinat.Groups.Thompson.F where + +-------------------------------------------------------------------------------- + +import Data.List + +import Math.Combinat.Classes +import Math.Combinat.ASCII + +import Math.Combinat.Trees.Binary ( BinTree ) +import qualified Math.Combinat.Trees.Binary as B + +-------------------------------------------------------------------------------- +-- * Tree diagrams + +-- | A tree diagram, consisting of two binary trees with the same number of leaves, +-- representing an element of the group F. +data TDiag = TDiag + { _width :: !Int -- ^ the width is the number of leaves, minus 1, of both diagrams + , _domain :: !T -- ^ the top diagram correspond to the /domain/ + , _range :: !T -- ^ the bottom diagram corresponds to the /range/ + } + deriving (Eq,Ord,Show) + +instance DrawASCII TDiag where + ascii = asciiTDiag + +instance HasWidth TDiag where + width = _width + +-- | Creates a tree diagram from two trees +mkTDiag :: T -> T -> TDiag +mkTDiag d1 d2 = reduce $ mkTDiagDontReduce d1 d2 + +-- | Creates a tree diagram, but does not reduce it. +mkTDiagDontReduce :: T -> T -> TDiag +mkTDiagDontReduce top bot = + if w1 == w2 + then TDiag w1 top bot + else error "mkTDiag: widths do not match" + where + w1 = treeWidth top + w2 = treeWidth bot + + +isValidTDiag :: TDiag -> Bool +isValidTDiag (TDiag w top bot) = (treeWidth top == w && treeWidth bot == w) + +isPositive :: TDiag -> Bool +isPositive (TDiag w top bot) = (bot == rightVine w) + +isReduced :: TDiag -> Bool +isReduced diag = (reduce diag == diag) + +-- | The generator x0 +x0 :: TDiag +x0 = TDiag 2 top bot where + top = branch caret leaf + bot = branch leaf caret + +-- | The generator x1 +x1 :: TDiag +x1 = xk 1 + +-- | The generators x0, x1, x2 ... +xk :: Int -> TDiag +xk = go where + go k | k< 0 = error "xk: negative indexed generator" + | k==0 = x0 + | otherwise = let TDiag _ t b = go (k-1) + in TDiag (k+2) (branch leaf t) (branch leaf b) + +-- | The identity element in the group F +identity :: TDiag +identity = TDiag 0 Lf Lf + +-- | A /positive diagram/ is a diagram whose bottom tree (the range) is a right vine. +positive :: T -> TDiag +positive t = TDiag w t (rightVine w) where w = treeWidth t + +-- | Swaps the top and bottom of a tree diagram. This is the inverse in the group F. +-- (Note: we don't do reduction here, as this operation keeps the reducedness) +inverse :: TDiag -> TDiag +inverse (TDiag w top bot) = TDiag w bot top + +-- | Decides whether two (possibly unreduced) tree diagrams represents the same group element in F. +equivalent :: TDiag -> TDiag -> Bool +equivalent diag1 diag2 = (identity == reduce (compose diag1 (inverse diag2))) + +-------------------------------------------------------------------------------- +-- * Reduction of tree diagrams + +-- | Reduces a diagram. The result is a normal form of an element in the group F. +reduce :: TDiag -> TDiag +reduce = worker where + + worker :: TDiag -> TDiag + worker diag = case step diag of + Nothing -> diag + Just diag' -> worker diag' + + step :: TDiag -> Maybe TDiag + step (TDiag w top bot) = + if null idxs + then Nothing + else Just $ TDiag w' top' bot' + where + cs1 = treeCaretList top + cs2 = treeCaretList bot + idxs = sortedIntersect cs1 cs2 + w' = w - length idxs + top' = removeCarets idxs top + bot' = removeCarets idxs bot + + -- | Intersects sorted lists + sortedIntersect :: [Int] -> [Int] -> [Int] + sortedIntersect = go where + go [] _ = [] + go _ [] = [] + go xxs@(x:xs) yys@(y:ys) = case compare x y of + LT -> go xs yys + EQ -> x : go xs ys + GT -> go xxs ys + +-- | List of carets at the bottom of the tree, indexed by their left edge position +treeCaretList :: T -> [Int] +treeCaretList = snd . go 0 where + go !x t = case t of + Lf -> (x+1 , [] ) + Ct -> (x+2 , [x] ) + Br t1 t2 -> (x2 , cs1++cs2) where + (x1 , cs1) = go x t1 + (x2 , cs2) = go x1 t2 + +-- | Remove the carets with the given indices +-- (throws an error if there is no caret at the given index) +removeCarets :: [Int] -> T -> T +removeCarets idxs tree = if null rem then final else error ("removeCarets: some stuff remained: " ++ show rem) where + + (_,rem,final) = go 0 idxs tree where + + go :: Int -> [Int] -> T -> (Int,[Int],T) + go !x [] t = (x + treeWidth t , [] , t) + go !x iis@(i:is) t = case t of + Lf -> (x+1 , iis , t) + Ct -> if x==i then (x+2 , is , Lf) else (x+2 , iis , Ct) + Br t1 t2 -> (x2 , iis2 , Br t1' t2') where + (x1 , iis1 , t1') = go x iis t1 + (x2 , iis2 , t2') = go x1 iis1 t2 + +-------------------------------------------------------------------------------- +-- * Composition of tree diagrams + +-- | If @diag1@ corresponds to the PL function @f@, and @diag2@ to @g@, then @compose diag1 diag2@ +-- will correspond to @(g.f)@ (note that the order is opposite than normal function composition!) +-- +-- This is the multiplication in the group F. +-- +compose :: TDiag -> TDiag -> TDiag +compose d1 d2 = reduce (composeDontReduce d1 d2) + +-- | Compose two tree diagrams without reducing the result +composeDontReduce :: TDiag -> TDiag -> TDiag +composeDontReduce (TDiag w1 top1 bot1) (TDiag w2 top2 bot2) = new where + new = mkTDiagDontReduce top' bot' + (list1,list2) = extensionToCommonTree bot1 top2 + top' = listGraft list1 top1 + bot' = listGraft list2 bot2 + +-- | Given two binary trees, we return a pair of list of subtrees which, grafted the to leaves of +-- the first (resp. the second) tree, results in the same extended tree. +extensionToCommonTree :: T -> T -> ([T],[T]) +extensionToCommonTree t1 t2 = snd $ go (0,0) (t1,t2) where + go (!x1,!x2) (!t1,!t2) = + case (t1,t2) of + ( Lf , Lf ) -> ( (x1+n1 , x2+n2 ) , ( [Lf] , [Lf] ) ) + ( Lf , Br _ _ ) -> ( (x1+n1 , x2+n2 ) , ( [t2] , replicate n2 Lf ) ) + ( Br _ _ , Lf ) -> ( (x1+n1 , x2+n2 ) , ( replicate n1 Lf , [t1] ) ) + ( Br l1 r1 , Br l2 r2 ) + -> let ( (x1' ,x2' ) , (ps1,ps2) ) = go (x1 ,x2 ) (l1,l2) + ( (x1'',x2'') , (qs1,qs2) ) = go (x1',x2') (r1,r2) + in ( (x1'',x2'') , (ps1++qs1, ps2++qs2) ) + where + n1 = numberOfLeaves t1 + n2 = numberOfLeaves t2 + +-------------------------------------------------------------------------------- +-- * Subdivions + +-- | Returns the list of dyadic subdivision points +subdivision1 :: T -> [Rational] +subdivision1 = go 0 1 where + go !a !b t = case t of + Leaf _ -> [a,b] + Branch l r -> go a c l ++ tail (go c b r) where c = (a+b)/2 + +-- | Returns the list of dyadic intervals +subdivision2 :: T -> [(Rational,Rational)] +subdivision2 = go 0 1 where + go !a !b t = case t of + Leaf _ -> [(a,b)] + Branch l r -> go a c l ++ go c b r where c = (a+b)/2 + + +-------------------------------------------------------------------------------- +-- * Binary trees + +-- | A (strict) binary tree with labelled leaves (but unlabelled nodes) +data Tree a + = Branch !(Tree a) !(Tree a) + | Leaf !a + deriving (Eq,Ord,Show,Functor) + +-- | The monadic join operation of binary trees +graft :: Tree (Tree a) -> Tree a +graft = go where + go (Branch l r) = Branch (go l) (go r) + go (Leaf t ) = t + +-- | A list version of 'graft' +listGraft :: [Tree a] -> Tree b -> Tree a +listGraft subs big = snd $ go subs big where + go ggs@(g:gs) t = case t of + Leaf _ -> (gs,g) + Branch l r -> (gs2, Branch l' r') where + (gs1,l') = go ggs l + (gs2,r') = go gs1 r + +-- | A completely unlabelled binary tree +type T = Tree () + +instance DrawASCII T where + ascii = asciiT + +instance HasNumberOfLeaves (Tree a) where + numberOfLeaves = treeNumberOfLeaves + +instance HasWidth (Tree a) where + width = treeWidth + +leaf :: T +leaf = Leaf () + +branch :: T -> T -> T +branch = Branch + +caret :: T +caret = branch leaf leaf + +treeNumberOfLeaves :: Tree a -> Int +treeNumberOfLeaves = go where + go (Branch l r) = go l + go r + go (Leaf _ ) = 1 + +-- | The width of the tree is the number of leaves minus 1. +treeWidth :: Tree a -> Int +treeWidth t = numberOfLeaves t - 1 + +-- | Enumerates the leaves a tree, starting from 0 +enumerate_ :: Tree a -> Tree Int +enumerate_ = snd . enumerate + +-- | Enumerates the leaves a tree, and also returns the number of leaves +enumerate :: Tree a -> (Int, Tree Int) +enumerate = go 0 where + go !k t = case t of + Leaf _ -> (k+1 , Leaf k) + Branch l r -> let (k' ,l') = go k l + (k'',r') = go k' r + in (k'', Branch l' r') + +-- | \"Right vine\" of the given width +rightVine :: Int -> T +rightVine k + | k< 0 = error "rightVine: negative width" + | k==0 = leaf + | otherwise = branch leaf (rightVine (k-1)) + +-- | \"Left vine\" of the given width +leftVine :: Int -> T +leftVine k + | k< 0 = error "leftVine: negative width" + | k==0 = leaf + | otherwise = branch (leftVine (k-1)) leaf + +-- | Flips each node of a binary tree +flipTree :: Tree a -> Tree a +flipTree = go where + go t = case t of + Leaf _ -> t + Branch l r -> Branch (go r) (go l) + +-------------------------------------------------------------------------------- +-- * Conversion to\/from BinTree + +-- | 'Tree' and 'BinTree' are the same type, except that 'Tree' is strict. +-- +-- TODO: maybe unify these two types? Until that, you can convert between the two +-- with these functions if necessary. +-- +toBinTree :: Tree a -> B.BinTree a +toBinTree = go where + go (Branch l r) = B.Branch (go l) (go r) + go (Leaf y ) = B.Leaf y + +fromBinTree :: B.BinTree a -> Tree a +fromBinTree = go where + go (B.Branch l r) = Branch (go l) (go r) + go (B.Leaf y ) = Leaf y + +-------------------------------------------------------------------------------- +-- * Pattern synonyms + +pattern Lf = Leaf () +pattern Br l r = Branch l r +pattern Ct = Br Lf Lf +pattern X0 = TDiag 2 (Br Ct Lf) (Br Lf Ct) +pattern X1 = TDiag 3 (Br Lf (Br Ct Lf)) (Br Lf (Br Lf Ct)) + +-------------------------------------------------------------------------------- +-- * ASCII + +-- | Draws a binary tree, with all leaves at the same (bottom) row +asciiT :: T -> ASCII +asciiT = asciiT' False + +-- | Draws a binary tree; when the boolean flag is @True@, we draw upside down +asciiT' :: Bool -> T -> ASCII +asciiT' inv = go where + + go t = case t of + Leaf _ -> emptyRect + Branch l r -> + if yl >= yr + then pasteOnto (yl+yr+1,if inv then yr else 0) (rs $ yl+1) $ + vcat HCenter + (bc $ yr+1) + (hcat bot al ar) + else pasteOnto (yl, if inv then yl else 0) (ls $ yr+1) $ + vcat HCenter + (bc $ yl+1) + (hcat bot al ar) + where + al = go l + ar = go r + yl = asciiYSize al + yr = asciiYSize ar + + bot = if inv then VTop else VBottom + hcat align p q = hCatWith align (HSepString " ") [p,q] + vcat align p q = vCatWith align VSepEmpty $ if inv then [q,p] else [p,q] + bc = if inv then asciiBigInvCaret else asciiBigCaret + ls = if inv then asciiBigRightSlope else asciiBigLeftSlope + rs = if inv then asciiBigLeftSlope else asciiBigRightSlope + + asciiBigCaret :: Int -> ASCII + asciiBigCaret k = hCatWith VTop HSepEmpty [ asciiBigLeftSlope k , asciiBigRightSlope k ] + + asciiBigInvCaret :: Int -> ASCII + asciiBigInvCaret k = hCatWith VTop HSepEmpty [ asciiBigRightSlope k , asciiBigLeftSlope k ] + + asciiBigLeftSlope :: Int -> ASCII + asciiBigLeftSlope k = if k>0 + then asciiFromLines [ replicate l ' ' ++ "/" | l<-[k-1,k-2..0] ] + else emptyRect + + asciiBigRightSlope :: Int -> ASCII + asciiBigRightSlope k = if k>0 + then asciiFromLines [ replicate l ' ' ++ "\\" | l<-[0..k-1] ] + else emptyRect + +-- | Draws a binary tree, with all leaves at the same (bottom) row, and labelling +-- the leaves starting with 0 (continuing with letters after 9) +asciiTLabels :: T -> ASCII +asciiTLabels = asciiTLabels' False + +-- | When the flag is true, we draw upside down +asciiTLabels' :: Bool -> T -> ASCII +asciiTLabels' inv t = + if inv + then vCatWith HLeft VSepEmpty [ labels , asciiT' inv t ] + else vCatWith HLeft VSepEmpty [ asciiT' inv t , labels ] + where + w = treeWidth t + labels = asciiFromString $ intersperse ' ' $ take (w+1) allLabels + allLabels = ['0'..'9'] ++ ['a'..'z'] + +-- | Draws a tree diagram +asciiTDiag :: TDiag -> ASCII +asciiTDiag (TDiag _ top bot) = vCatWith HLeft (VSepString " ") [asciiT' False top , asciiT' True bot] + +-------------------------------------------------------------------------------- + + diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs --- a/Math/Combinat/Helper.hs +++ b/Math/Combinat/Helper.hs @@ -1,6 +1,7 @@ -- | Miscellaneous helper functions +{-# LANGUAGE BangPatterns, PolyKinds #-} module Math.Combinat.Helper where -------------------------------------------------------------------------------- @@ -9,6 +10,7 @@ import Data.List import Data.Ord +import Data.Proxy import Data.Set (Set) ; import qualified Data.Set as Set import Data.Map (Map) ; import qualified Data.Map as Map @@ -22,6 +24,24 @@ debug x y = trace ("-- " ++ show x ++ "\n") y -------------------------------------------------------------------------------- +-- * proxy + +proxyUndef :: Proxy a -> a +proxyUndef _ = error "proxyUndef" + +proxyOf :: a -> Proxy a +proxyOf _ = Proxy + +proxyOf1 :: f a -> Proxy a +proxyOf1 _ = Proxy + +proxyOf2 :: g (f a) -> Proxy a +proxyOf2 _ = Proxy + +asProxyTypeOf1 :: f a -> Proxy a -> f a +asProxyTypeOf1 y _ = y + +-------------------------------------------------------------------------------- -- * pairs swap :: (a,b) -> (b,a) @@ -73,6 +93,41 @@ | otherwise = x : worker (Set.insert x s) xs -------------------------------------------------------------------------------- +-- * increasing \/ decreasing sequences + +{-# SPECIALIZE isWeaklyIncreasing :: [Int] -> Bool #-} +isWeaklyIncreasing :: Ord a => [a] -> Bool +isWeaklyIncreasing = go where + go xs = case xs of + (a:rest@(b:_)) -> a <= b && go rest + [_] -> True + [] -> True + +{-# SPECIALIZE isStrictlyIncreasing :: [Int] -> Bool #-} +isStrictlyIncreasing :: Ord a => [a] -> Bool +isStrictlyIncreasing = go where + go xs = case xs of + (a:rest@(b:_)) -> a < b && go rest + [_] -> True + [] -> True + +{-# SPECIALIZE isWeaklyDecreasing :: [Int] -> Bool #-} +isWeaklyDecreasing :: Ord a => [a] -> Bool +isWeaklyDecreasing = go where + go xs = case xs of + (a:rest@(b:_)) -> a >= b && go rest + [_] -> True + [] -> True + +{-# SPECIALIZE isStrictlyDecreasing :: [Int] -> Bool #-} +isStrictlyDecreasing :: Ord a => [a] -> Bool +isStrictlyDecreasing = go where + go xs = case xs of + (a:rest@(b:_)) -> a > b && go rest + [_] -> True + [] -> True + +-------------------------------------------------------------------------------- -- * first \/ last -- | The boolean argument will @True@ only for the last element @@ -150,8 +205,8 @@ -- iterated function application nest :: Int -> (a -> a) -> a -> a -nest 0 _ x = x -nest n f x = nest (n-1) f (f x) +nest !0 _ x = x +nest !n f x = nest (n-1) f (f x) unfold1 :: (a -> Maybe a) -> a -> [a] unfold1 f x = case f x of diff --git a/Math/Combinat/LatticePaths.hs b/Math/Combinat/LatticePaths.hs --- a/Math/Combinat/LatticePaths.hs +++ b/Math/Combinat/LatticePaths.hs @@ -15,6 +15,7 @@ import Data.List import System.Random +import Math.Combinat.Classes import Math.Combinat.Numbers import Math.Combinat.Trees.Binary import Math.Combinat.ASCII as ASCII @@ -90,6 +91,12 @@ go !h !y (t:ts) = case t of UpStep -> go (max h (y+1)) (y+1) ts DownStep -> go h (y-1) ts + +instance HasHeight LatticePath where + height = pathHeight + +instance HasWidth LatticePath where + width = length -- | Endpoint of a lattice path, which starts from @(0,0)@. pathEndpoint :: LatticePath -> (Int,Int) diff --git a/Math/Combinat/Numbers/Primes.hs b/Math/Combinat/Numbers/Primes.hs --- a/Math/Combinat/Numbers/Primes.hs +++ b/Math/Combinat/Numbers/Primes.hs @@ -22,6 +22,8 @@ , powerMod -- * Prime testing , millerRabinPrimalityTest + , isProbablyPrime + , isVeryProbablyPrime ) where @@ -32,6 +34,8 @@ import Data.List ( group , sort ) import Data.Bits +import System.Random + -------------------------------------------------------------------------------- -- List of prime numbers @@ -286,5 +290,65 @@ {-# SPECIALIZE powMod :: Integer -> Integer -> Integer -> Integer #-} powMod :: Integral a => a -> a -> a -> a powMod m = pow' (mulMod m) (squareMod m) + +-------------------------------------------------------------------------------- + +-- | For very small numbers, we use trial division, for larger numbers, we apply the +-- Miller-Rabin primality test @log4(n)@ times, with candidate witnesses derived +-- deterministically from @n@ using a pseudo-random sequence +-- (which /should be/ based on a cryptographic hash function, but isn\'t, yet). +-- +-- Thus the candidate witnesses should behave essentially like random, but the +-- resulting function is still a deterministic, pure function. +-- +-- TODO: implement the hash sequence, at the moment we use 'System.Random' instead... +-- +isProbablyPrime :: Integer -> Bool +isProbablyPrime n + | n < 2 = False + | even n = (n==2) + | n < 1000 = length (integerFactorsTrialDivision n) == 1 + | otherwise = and [ millerRabinPrimalityTest n a | a <- witnessList ] + where + log2n = integerLog2 n + nchecks = 1 + fromInteger (div log2n 2) :: Int + witnessList = take nchecks pseudoRnds + pseudoRnds = 2 : [ a | a <- integerRndSequence n , a > 1 && a < (n-1) ] + +-- | A more exhaustive version of 'isProbablyPrime', this one tests candidate +-- witnesses both the first log4(n) prime numbers and then log4(n) pseudo-random +-- numbers +isVeryProbablyPrime :: Integer -> Bool +isVeryProbablyPrime n + | n < 2 = False + | even n = (n==2) + | n < 1000 = length (integerFactorsTrialDivision n) == 1 + | otherwise = and [ millerRabinPrimalityTest n a | a <- witnessList ] + where + log2n = integerLog2 n + nchecks = 1 + fromInteger (div log2n 2) :: Int + witnessList = take nchecks primes ++ take nchecks pseudoRnds + pseudoRnds = [ a | a <- integerRndSequence (n+3) , a > 1 && a < (n-1) ] + +-------------------------------------------------------------------------------- + +{- +-- | Given an integer @n@, we return an infinite sequence of pseudo-random integers +-- between @0..n-1@, generated using a crypographic hash function. +-- +integerHashSequence :: Integer -> [Integer] +integerHashSequence = error "integerHashSequence: not implemented yet" +-} + +-- | Given an integer @n@, we initialize a system random generator with using a +-- seed derived from @n@ (note that this uses at most 32 or 64 bits), and generate +-- an infinite sequence of pseudo-random integers between @0..n-1@, generated by +-- that random generator. +-- +-- Note that this is not really a preferred way of generating such sequences! +-- +integerRndSequence :: Integer -> [Integer] +integerRndSequence n = randomRs (0,n-1) gen where + gen = mkStdGen $ fromInteger (n + 17 * integerLog2 n) -------------------------------------------------------------------------------- diff --git a/Math/Combinat/Numbers/Series.hs b/Math/Combinat/Numbers/Series.hs --- a/Math/Combinat/Numbers/Series.hs +++ b/Math/Combinat/Numbers/Series.hs @@ -143,7 +143,7 @@ [ b m * product [ (a i)^j | (i,j)<-es ] * fromInteger (multinomial (map snd es)) | p <- partitions n , let es = toExponentialForm p - , let m = width p + , let m = partitionWidth p ] -------------------------------------------------------------------------------- @@ -154,8 +154,8 @@ lagrangeCoeff p = div numer denom where numer = (-1)^m * product (map fromIntegral [n+1..n+m]) denom = fromIntegral (n+1) * product (map (factorial . snd) es) - m = width p - n = weight p + m = partitionWidth p + n = partitionWeight p es = toExponentialForm p -- | We expect the input series to match @(0:1:_)@. The following is true for the result (at least with exact arithmetic): diff --git a/Math/Combinat/Partitions/Integer.hs b/Math/Combinat/Partitions/Integer.hs --- a/Math/Combinat/Partitions/Integer.hs +++ b/Math/Combinat/Partitions/Integer.hs @@ -27,9 +27,10 @@ -- import Data.Map (Map) -- import qualified Data.Map as Map -import Math.Combinat.Helper +import Math.Combinat.Classes import Math.Combinat.ASCII as ASCII import Math.Combinat.Numbers (factorial,binomial,multinomial) +import Math.Combinat.Helper -------------------------------------------------------------------------------- -- * Type and basic stuff @@ -38,11 +39,6 @@ -- are monotone decreasing sequences of /positive/ integers. The @Ord@ instance is lexicographical. newtype Partition = Partition [Int] deriving (Eq,Ord,Show,Read) ---------------------------------------------------------------------------------- - -class HasNumberOfParts p where - numberOfParts :: p -> Int - instance HasNumberOfParts Partition where numberOfParts (Partition p) = length p @@ -70,31 +66,53 @@ isPartition [x] = x > 0 isPartition (x:xs@(y:_)) = (x >= y) && isPartition xs +isEmptyPartition :: Partition -> Bool +isEmptyPartition (Partition p) = null p + +emptyPartition :: Partition +emptyPartition = Partition [] + +instance CanBeEmpty Partition where + empty = emptyPartition + isEmpty = isEmptyPartition + fromPartition :: Partition -> [Int] fromPartition (Partition part) = part -- | The first element of the sequence. -height :: Partition -> Int -height (Partition part) = case part of +partitionHeight :: Partition -> Int +partitionHeight (Partition part) = case part of (p:_) -> p - [] -> 0 + [] -> 0 -- | The length of the sequence (that is, the number of parts). -width :: Partition -> Int -width (Partition part) = length part +partitionWidth :: Partition -> Int +partitionWidth (Partition part) = length part +instance HasHeight Partition where + height = partitionHeight + +instance HasWidth Partition where + width = partitionWidth + heightWidth :: Partition -> (Int,Int) heightWidth part = (height part, width part) -- | The weight of the partition -- (that is, the sum of the corresponding sequence). -weight :: Partition -> Int -weight (Partition part) = sum' part +partitionWeight :: Partition -> Int +partitionWeight (Partition part) = sum' part +instance HasWeight Partition where + weight = partitionWeight + -- | The dual (or conjugate) partition. dualPartition :: Partition -> Partition dualPartition (Partition part) = Partition (_dualPartition part) +instance HasDuality Partition where + dual = dualPartition + data Pair = Pair !Int !Int _dualPartition :: [Int] -> [Int] @@ -375,14 +393,21 @@ go !h !n = [ a:as | a<-[1..min n h], as <- go (a-1) (n-a) ] -------------------------------------------------------------------------------- --- * Sub-partitions of a given partition +-- * Sub- and super-partitions of a given partition -- | Returns @True@ of the first partition is a subpartition (that is, fit inside) of the second. -- This includes equality isSubPartitionOf :: Partition -> Partition -> Bool isSubPartitionOf (Partition ps) (Partition qs) = and $ zipWith (<=) ps (qs ++ repeat 0) +-- | This is provided for convenience\/completeness only, as: +-- +-- > isSuperPartitionOf q p == isSubPartitionOf p q +-- +isSuperPartitionOf :: Partition -> Partition -> Bool +isSuperPartitionOf (Partition qs) (Partition ps) = and $ zipWith (<=) ps (qs ++ repeat 0) + -- | Sub-partitions of a given partition with the given weight: -- -- > sort (subPartitions d q) == sort [ p | p <- partitions d, isSubPartitionOf p q ] @@ -421,6 +446,31 @@ | h==0 = [] | otherwise = [] : [ this:rest | this <- [1..min h b] , rest <- go this bs ] +---------------------------------------- + +-- | Super-partitions of a given partition with the given weight: +-- +-- > sort (superPartitions d p) == sort [ q | q <- partitions d, isSubPartitionOf p q ] +-- +superPartitions :: Int -> Partition -> [Partition] +superPartitions d (Partition ps) = map Partition (_superPartitions d ps) + +_superPartitions :: Int -> [Int] -> [[Int]] +_superPartitions dd small + | dd < w0 = [] + | null small = _partitions dd + | otherwise = go dd w1 dd (small ++ repeat 0) + where + w0 = sum' small + w1 = w0 - head small + -- d = remaining weight of the outer partition we are constructing + -- w = remaining weight of the inner partition (we need to reserve at least this amount) + -- h = max height (decreasing) + go !d !w !h (!a:as@(b:_)) + | d < 0 = [] + | d == 0 = if a == 0 then [[]] else [] + | otherwise = [ this:rest | this <- [max 1 a .. min h (d-w)] , rest <- go (d-this) (w-b) this as ] + -------------------------------------------------------------------------------- -- * The Pieri rule @@ -452,7 +502,10 @@ dualPieriRule :: Partition -> Int -> [Partition] dualPieriRule lam n = map dualPartition $ pieriRule (dualPartition lam) n -{- + +{- +-- moved to "Math.Combinat.Tableaux.GelfandTsetlin" + -- | Computes the Schur expansion of @h[n1]*h[n2]*h[n3]*...*h[nk]@ via iterating the Pieri rule iteratedPieriRule :: Num coeff => [Int] -> Map Partition coeff iteratedPieriRule = iteratedPieriRule' (Partition []) @@ -531,7 +584,7 @@ -------------------------------------------------------------------------------- -{- +#ifdef QUICKCHECK -- * Tests @@ -567,6 +620,6 @@ prop_dominating_list :: Partition -> Bool prop_dominating_list mu = (dominatingPartitions mu == [ lam | lam <- partitions (weight mu ), lam `dominates` mu ]) --} +#endif -------------------------------------------------------------------------------- diff --git a/Math/Combinat/Partitions/NonCrossing.hs b/Math/Combinat/Partitions/NonCrossing.hs --- a/Math/Combinat/Partitions/NonCrossing.hs +++ b/Math/Combinat/Partitions/NonCrossing.hs @@ -30,7 +30,7 @@ import Math.Combinat.LatticePaths import Math.Combinat.Helper import Math.Combinat.Partitions.Set -import Math.Combinat.Partitions ( HasNumberOfParts(..) ) +import Math.Combinat.Classes -------------------------------------------------------------------------------- -- * The type of non-crossing partitions diff --git a/Math/Combinat/Partitions/Plane.hs b/Math/Combinat/Partitions/Plane.hs --- a/Math/Combinat/Partitions/Plane.hs +++ b/Math/Combinat/Partitions/Plane.hs @@ -24,6 +24,7 @@ import Data.List import Data.Array +import Math.Combinat.Classes import Math.Combinat.Partitions import Math.Combinat.Tableaux as Tableaux import Math.Combinat.Helper @@ -49,6 +50,10 @@ y = length pps x = maximum (map length pps) +instance CanBeEmpty PlanePart where + isEmpty = null . fromPlanePart + empty = PlanePart [] + -- | Throws an exception if the input is not a plane partition toPlanePart :: [[Int]] -> PlanePart toPlanePart pps = if isValidPlanePart pps @@ -57,7 +62,7 @@ -- | The XY projected shape of a plane partition, as an integer partition planePartShape :: PlanePart -> Partition -planePartShape = Tableaux.shape . fromPlanePart +planePartShape = Tableaux.tableauShape . fromPlanePart -- | The Z height of a plane partition planePartZHeight :: PlanePart -> Int @@ -68,6 +73,9 @@ planePartWeight :: PlanePart -> Int planePartWeight (PlanePart xs) = sum' (map sum' xs) + +instance HasWeight PlanePart where + weight = planePartWeight -------------------------------------------------------------------------------- -- * constructing plane partitions diff --git a/Math/Combinat/Partitions/Set.hs b/Math/Combinat/Partitions/Set.hs --- a/Math/Combinat/Partitions/Set.hs +++ b/Math/Combinat/Partitions/Set.hs @@ -17,7 +17,7 @@ import Math.Combinat.Sets import Math.Combinat.Numbers import Math.Combinat.Helper -import Math.Combinat.Partitions ( HasNumberOfParts(..) ) +import Math.Combinat.Classes -------------------------------------------------------------------------------- -- * The type of set partitions diff --git a/Math/Combinat/Partitions/Skew.hs b/Math/Combinat/Partitions/Skew.hs --- a/Math/Combinat/Partitions/Skew.hs +++ b/Math/Combinat/Partitions/Skew.hs @@ -3,20 +3,30 @@ -- -- Skew partitions are the difference of two integer partitions, denoted by @lambda/mu@. -- +-- For example +-- +-- > mkSkewPartition (Partition [9,7,3,2,2,1] , Partition [5,3,2,1]) +-- +-- creates the skew partition @(9,7,3,2,2,1) / (5,3,2,1)@, which looks like +-- +-- <> +-- -{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP, BangPatterns #-} module Math.Combinat.Partitions.Skew where -------------------------------------------------------------------------------- import Data.List +import Math.Combinat.Classes import Math.Combinat.Partitions.Integer import Math.Combinat.ASCII -------------------------------------------------------------------------------- +-- * Basics --- | A skew partition @lambda/mu@ is represented by the list @[ (mu_i , lambda_i-mu_i) | i<-[1..n] ]@ +-- | A skew partition @lambda/mu@ is internally represented by the list @[ (mu_i , lambda_i-mu_i) | i<-[1..n] ]@ newtype SkewPartition = SkewPartition [(Int,Int)] deriving (Eq,Ord,Show) -- | @mkSkewPartition (lambda,mu)@ creates the skew partition @lambda/mu@. @@ -32,9 +42,14 @@ then Just $ SkewPartition $ zipWith (\b a -> (a,b-a)) bs (as ++ repeat 0) else Nothing +-- | The weight of a skew partition is the weight of the outer partition minus the +-- the weight of the inner partition (that is, the number of boxes present). skewPartitionWeight :: SkewPartition -> Int skewPartitionWeight (SkewPartition abs) = foldl' (+) 0 (map snd abs) +instance HasWeight SkewPartition where + weight = skewPartitionWeight + -- | This function \"cuts off\" the \"uninteresting parts\" of a skew partition normalizeSkewPartition :: SkewPartition -> SkewPartition normalizeSkewPartition (SkewPartition abs) = SkewPartition abs' where @@ -43,23 +58,61 @@ k = length (takeWhile (==0) bs) abs' = zip [ a-a0 | a <- drop k as ] (drop k bs) --- | Returns the outer and inner partition of a skew partition, respectively +-- | Returns the outer and inner partition of a skew partition, respectively: +-- +-- > mkSkewPartition . fromSkewPartition == id +-- fromSkewPartition :: SkewPartition -> (Partition,Partition) fromSkewPartition (SkewPartition list) = (toPartition (zipWith (+) as bs) , toPartition (filter (>0) as)) where (as,bs) = unzip list +-- | The @lambda@ part of @lambda/mu@ outerPartition :: SkewPartition -> Partition outerPartition = fst . fromSkewPartition +-- | The @mu@ part of @lambda/mu@ innerPartition :: SkewPartition -> Partition innerPartition = snd . fromSkewPartition +-- | The dual skew partition (that is, the mirror image to the main diagonal) dualSkewPartition :: SkewPartition -> SkewPartition dualSkewPartition = mkSkewPartition . f . fromSkewPartition where f (lam,mu) = (dualPartition lam, dualPartition mu) +instance HasDuality SkewPartition where + dual = dualSkewPartition + -------------------------------------------------------------------------------- +-- * Listing skew partitions +-- | Lists all skew partitions with the given outer shape and given (skew) weight +skewPartitionsWithOuterShape :: Partition -> Int -> [SkewPartition] +skewPartitionsWithOuterShape outer skewWeight + | innerWeight < 0 || innerWeight > outerWeight = [] + | otherwise = [ mkSkewPartition (outer,inner) | inner <- subPartitions innerWeight outer ] + where + outerWeight = weight outer + innerWeight = outerWeight - skewWeight + +-- | Lists all skew partitions with the given outer shape and any (skew) weight +allSkewPartitionsWithOuterShape :: Partition -> [SkewPartition] +allSkewPartitionsWithOuterShape outer + = concat [ skewPartitionsWithOuterShape outer w | w<-[0..outerWeight] ] + where + outerWeight = weight outer + +-- | Lists all skew partitions with the given inner shape and given (skew) weight +skewPartitionsWithInnerShape :: Partition -> Int -> [SkewPartition] +skewPartitionsWithInnerShape inner skewWeight + | innerWeight > outerWeight = [] + | otherwise = [ mkSkewPartition (outer,inner) | outer <- superPartitions outerWeight inner ] + where + outerWeight = innerWeight + skewWeight + innerWeight = weight inner + +-------------------------------------------------------------------------------- +-- * ASCII + asciiSkewFerrersDiagram :: SkewPartition -> ASCII asciiSkewFerrersDiagram = asciiSkewFerrersDiagram' ('@','.') EnglishNotation @@ -80,3 +133,38 @@ -------------------------------------------------------------------------------- +#ifdef QUICKCHECK + +prop_dual_dual :: SkewPartition -> Bool +prop_dual_dual sp = (dualSkewPartition (dualSkewPartition sp) == sp) + +prop_dual_from :: SkewPartition -> Bool +prop_dual_from sp = (p==p' && q==q') where + (p,q) = fromSkewPartition sp + sp' = dualSkewPartition sp + (p',q') = fromSkewPartition sp' + +prop_from_to :: SkewPartition -> Bool +prop_from_to sp = (mkSkewPartition (fromSkewPartition sp) == sp) + +prop_to_from :: (Partition,Partition) -> Bool +prop_to_from (p,q) = + case mb of + Nothing -> True + Just sp -> fromSkewPartition sp == (p,q) + where + mb = safeSkewPartition (p,q) + +prop_from_to_from :: SkewPartition -> Bool +prop_from_to_from sp = (pq == pq') where + pq = fromSkewPartition sp + sp' = mkSkewPartition pq + pq' = fromSkewPartition sp' + +prop_weight :: SkewPartition -> Bool +prop_weight sp = (skewPartitionWeight sp == weight p - weight q) where + (p,q) = fromSkewPartition sp + +#endif + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs --- a/Math/Combinat/Permutations.hs +++ b/Math/Combinat/Permutations.hs @@ -1,36 +1,67 @@ --- | Permutations. See: --- Donald E. Knuth: The Art of Computer Programming, vol 4, pre-fascicle 2B. +-- | Permutations. -- -{-# OPTIONS_GHC -fno-warn-name-shadowing #-} -{-# LANGUAGE CPP, ScopedTypeVariables, GeneralizedNewtypeDeriving #-} +-- See eg.: +-- Donald E. Knuth: The Art of Computer Programming, vol 4, pre-fascicle 2B. +-- +-- WARNING: As of version 0.2.8.0, I changed the convention of how permutations +-- are represented internally. Also now they act on the /right/ by default! +-- + +{-# LANGUAGE CPP, ScopedTypeVariables, GeneralizedNewtypeDeriving, FlexibleContexts #-} module Math.Combinat.Permutations - ( -- * Types - Permutation - , DisjointCycles + ( -- * The Permutation type + Permutation (..) , fromPermutation , permutationArray - , toPermutationUnsafe - , arrayToPermutationUnsafe + , permutationUArray + , uarrayToPermutationUnsafe , isPermutation + , maybePermutation , toPermutation + , toPermutationUnsafe , permutationSize -- * Disjoint cycles + , DisjointCycles (..) , fromDisjointCycles , disjointCyclesUnsafe , permutationToDisjointCycles , disjointCyclesToPermutation + , numberOfCycles + -- * Queries + , isIdentityPermutation + , isReversePermutation , isEvenPermutation , isOddPermutation - , signOfPermutation -- , Sign(..) + , signOfPermutation + , signValueOfPermutation + , module Math.Combinat.Sign -- , Sign(..) , isCyclicPermutation + -- * Some concrete permutations + , transposition + , transpositions + , adjacentTransposition + , adjacentTranspositions + , cycleLeft + , cycleRight + , reversePermutation -- * Permutation groups - , permute - , permuteList - , multiply - , inverse , identity - -- * Simple permutations + , inverse + , multiply + , multiplyMany + , multiplyMany' + -- * Action of the permutation group + , permute + , permuteList + , permuteLeft , permuteRight + , permuteLeftList , permuteRightList + -- * ASCII drawing + , asciiPermutation + , asciiDisjointCycles + , twoLineNotation + , twoLineNotation' + -- * List of permutations , permutations , _permutations , permutationsNaive @@ -60,16 +91,17 @@ import Control.Monad import Control.Monad.ST -#if BASE_VERSION < 4 -import Data.List -#else import Data.List hiding (permutations) -#endif -import Data.Array +import Data.Array (Array) import Data.Array.ST +import Data.Array.Unboxed +import Data.Array.IArray +import Data.Array.MArray import Data.Array.Unsafe +import Math.Combinat.ASCII +import Math.Combinat.Classes import Math.Combinat.Helper import Math.Combinat.Sign import Math.Combinat.Numbers (factorial,binomial) @@ -83,48 +115,146 @@ -------------------------------------------------------------------------------- -- * Types --- | Standard notation for permutations. Internally it is an array of the integers @[1..n]@. -newtype Permutation = Permutation (Array Int Int) deriving (Eq,Ord,Show,Read) +-- | A permutation. Internally it is an (unboxed) array of the integers @[1..n]@. +-- +-- If this array of integers is @[p1,p2,...,pn]@, then in two-line +-- notations, that represents the permutation +-- +-- > ( 1 2 3 ... n ) +-- > ( p1 p2 p3 ... pn ) +-- +-- That is, it is the permutation @sigma@ whose (right) action on the set @[1..n]@ is +-- +-- > sigma(1) = p1 +-- > sigma(2) = p2 +-- > ... +-- +-- (NOTE: this changed at version 0.2.8.0!) +-- +newtype Permutation = Permutation (UArray Int Int) 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 + +instance Read Permutation where + readsPrec d r = readParen (d > 10) fun r where + fun r = [ (toPermutation p,t) + | ("toPermutation",s) <- lex r + , (p,t) <- readsPrec 11 s -- app_prec = 10 + ] + +instance DrawASCII Permutation where + ascii = asciiPermutation + -- | Disjoint cycle notation for permutations. Internally it is @[[Int]]@. +-- +-- The cycles are to be understood as follows: a cycle @[c1,c2,...,ck]@ means +-- the permutation +-- +-- > ( c1 c2 c3 ... ck ) +-- > ( c2 c3 c4 ... c1 ) +-- newtype DisjointCycles = DisjointCycles [[Int]] deriving (Eq,Ord,Show,Read) fromPermutation :: Permutation -> [Int] fromPermutation (Permutation ar) = elems ar +permutationUArray :: Permutation -> UArray Int Int +permutationUArray (Permutation ar) = ar + +-- | Note: this is slower than 'permutationUArray' permutationArray :: Permutation -> Array Int Int -permutationArray (Permutation ar) = ar +permutationArray (Permutation ar) = listArray (1,n) (elems ar) where + (1,n) = bounds ar -- | Assumes that the input is a permutation of the numbers @[1..n]@. toPermutationUnsafe :: [Int] -> Permutation toPermutationUnsafe xs = Permutation perm where - n = length xs + n = length xs perm = listArray (1,n) xs --- Indexing starts from 1. -arrayToPermutationUnsafe :: Array Int Int -> Permutation -arrayToPermutationUnsafe = Permutation +-- | Note: Indexing starts from 1. +uarrayToPermutationUnsafe :: UArray Int Int -> Permutation +uarrayToPermutationUnsafe = Permutation -- | Checks whether the input is a permutation of the numbers @[1..n]@. isPermutation :: [Int] -> Bool isPermutation xs = (ar!0 == 0) && and [ ar!j == 1 | j<-[1..n] ] where n = length xs -- the zero index is an unidiomatic hack - ar = accumArray (+) 0 (0,n) $ map f xs + ar = (accumArray (+) 0 (0,n) $ map f xs) :: UArray Int Int f :: Int -> (Int,Int) f j = if j<1 || j>n then (0,1) else (j,1) +-- | Checks whether the input is a permutation of the numbers @[1..n]@. +maybePermutation :: [Int] -> Maybe Permutation +maybePermutation input = runST action where + n = length input + 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) + go (j:js) = if j<1 || j>n + then return Nothing + else do + z <- readArray ar j + writeArray ar j (z+1) + if z==0 then go js + else return Nothing + go input + -- | Checks the input. toPermutation :: [Int] -> Permutation -toPermutation xs = if isPermutation xs - then toPermutationUnsafe xs - else error "toPermutation: not a permutation" +toPermutation xs = case maybePermutation xs of + Just p -> p + Nothing -> error "toPermutation: not a permutation" -- | Returns @n@, where the input is a permutation of the numbers @[1..n]@ permutationSize :: Permutation -> Int permutationSize (Permutation ar) = snd $ bounds ar +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 + -------------------------------------------------------------------------------- +-- * ASCII + +-- | Synonym for 'twoLineNotation' +asciiPermutation :: Permutation -> ASCII +asciiPermutation = twoLineNotation + +asciiDisjointCycles :: DisjointCycles -> ASCII +asciiDisjointCycles (DisjointCycles cycles) = final where + final = hCatWith VTop (HSepSpaces 1) boxes + boxes = [ twoLineNotation' (f cyc) | cyc <- cycles ] + f cyc = pairs (cyc ++ [head cyc]) + +-- | The standard two-line notation, moving the element indexed by the top row into +-- the place indexed by the corresponding element in the bottom row. +twoLineNotation :: Permutation -> ASCII +twoLineNotation (Permutation arr) = twoLineNotation' $ zip [1..] (elems arr) + +twoLineNotation' :: [(Int,Int)] -> ASCII +twoLineNotation' xys = asciiFromLines [ topLine, botLine ] where + topLine = "( " ++ intercalate " " us ++ " )" + botLine = "( " ++ intercalate " " vs ++ " )" + pairs = [ (show x, show y) | (x,y) <- xys ] + (us,vs) = unzip (map f pairs) + f (s,t) = (s',t') where + a = length s + b = length t + c = max a b + s' = replicate (c-a) ' ' ++ s + t' = replicate (c-b) ' ' ++ t + +-------------------------------------------------------------------------------- -- * Disjoint cycles fromDisjointCycles :: DisjointCycles -> [[Int]] @@ -132,20 +262,42 @@ disjointCyclesUnsafe :: [[Int]] -> DisjointCycles disjointCyclesUnsafe = DisjointCycles + +instance DrawASCII DisjointCycles where + ascii = asciiDisjointCycles + +instance HasNumberOfCycles DisjointCycles where + numberOfCycles (DisjointCycles cycles) = length cycles + +instance HasNumberOfCycles Permutation where + numberOfCycles = numberOfCycles . permutationToDisjointCycles disjointCyclesToPermutation :: Int -> DisjointCycles -> Permutation disjointCyclesToPermutation n (DisjointCycles cycles) = Permutation perm where + pairs :: [Int] -> [(Int,Int)] pairs xs@(x:_) = worker (xs++[x]) where worker (x:xs@(y:_)) = (x,y):worker xs worker _ = [] - perm = runST $ do + pairs [] = error "disjointCyclesToPermutation: empty cycle" + + perm = runSTUArray $ do ar <- newArray_ (1,n) :: ST s (STUArray s Int Int) forM_ [1..n] $ \i -> writeArray ar i i forM_ cycles $ \cyc -> forM_ (pairs cyc) $ \(i,j) -> writeArray ar i j - freeze ar + return ar -- freeze ar --- | This is compatible with Maple's @convert(perm,\'disjcyc\')@. +-- | Convert to disjoint cycle notation. +-- +-- This is compatible with Maple's @convert(perm,\'disjcyc\')@ +-- and also with Mathematica's @PermutationCycles[perm]@ +-- +-- Note however, that for example Mathematica uses the +-- /top row/ to represent a permutation, while we use the +-- /bottom row/ - thus even though this function looks +-- identical, the /meaning/ of both the input and output +-- is different! +-- permutationToDisjointCycles :: Permutation -> DisjointCycles permutationToDisjointCycles (Permutation perm) = res where @@ -165,7 +317,7 @@ step tag k = do cyc <- worker tag k k [k] m <- next tag (k+1) - return (reverse cyc,m) + return (reverse cyc, m) next :: STUArray s Int Bool -> Int -> ST s (Maybe Int) next tag k = if k > n @@ -221,6 +373,8 @@ False -> Minus -- | Plus 1 or minus 1. +{-# SPECIALIZE signValueOfPermutation :: Permutation -> Int #-} +{-# SPECIALIZE signValueOfPermutation :: Permutation -> Integer #-} signValueOfPermutation :: Num a => Permutation -> a signValueOfPermutation = signValue . signOfPermutation @@ -233,53 +387,123 @@ where n = permutationSize perm DisjointCycles cycles = permutationToDisjointCycles perm - + -------------------------------------------------------------------------------- --- * Permutation groups - --- | Action of a permutation on a set. If our permutation is --- encoded with the sequence @[p1,p2,...,pn]@, then in the --- two-line notation we have +-- * Some concrete permutations + +-- | 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] + +-- | 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 + +-- | A transposition (swapping two elements). -- --- > ( 1 2 3 ... n ) --- > ( p1 p2 p3 ... pn ) +-- @transposition n (i,j)@ is the permutation of size @n@ which swaps @i@\'th and @j@'th elements. -- --- We adopt the convention that permutations act /on the left/ --- (as opposed to Knuth, where they act on the right). --- Thus, --- --- > permute pi1 (permute pi2 set) == permute (pi1 `multiply` pi2) set --- --- The second argument should be an array with bounds @(1,n)@. --- The function checks the array bounds. -permute :: Permutation -> Array Int a -> Array Int a -permute pi@(Permutation perm) ar = - if (a==1) && (b==n) - then listArray (1,n) [ ar!(perm!i) | i <- [1..n] ] - else error "permute: array bounds do not match" +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] ] + else error "transposition: index out of range" where - (_,n) = bounds perm - (a,b) = bounds ar + f k | k == i = j + | k == j = i + | otherwise = k --- | The list should be of length @n@. -permuteList :: Permutation -> [a] -> [a] -permuteList perm xs = elems $ permute perm $ listArray (1,n) xs where - n = permutationSize perm +-- | Product of transpositions. +-- +-- > transpositions n list == multiplyMany [ transposition n pair | pair <- list ] +-- +transpositions :: Int -> [(Int,Int)] -> Permutation +transpositions n list = Permutation (runSTUArray action) where --- | Multiplies two permutations together. See 'permute' for our --- conventions. + action :: ST s (STUArray s Int Int) + action = do + arr <- newArray_ (1,n) + forM_ [1..n] $ \i -> writeArray arr i i + let doSwap (i,j) = do + u <- readArray arr i + v <- readArray arr j + writeArray arr i v + writeArray arr j u + mapM_ doSwap list + return arr + +-- | @adjacentTransposition n k@ swaps the elements @k@ and @(k+1)@. +adjacentTransposition :: Int -> Int -> Permutation +adjacentTransposition n k + | k>0 && k adjacentTranspositions n list == multiplyMany [ adjacentTransposition n idx | idx <- list ] +-- +adjacentTranspositions :: Int -> [Int] -> Permutation +adjacentTranspositions n list = Permutation (runSTUArray action) where + + action :: ST s (STUArray s Int Int) + action = do + arr <- newArray_ (1,n) + forM_ [1..n] $ \i -> writeArray arr i i + let doSwap i + | i<0 || i>=n = error "adjacentTranspositions: index out of range" + | otherwise = do + u <- readArray arr i + v <- readArray arr (i+1) + writeArray arr i v + writeArray arr (i+1) u + mapM_ doSwap list + return arr + +-- | The permutation which cycles a list left by one step: +-- +-- > permuteList (cycleLeft 5) "abcde" == "bcdea" +-- +-- Or in two-line notation: +-- +-- > ( 1 2 3 4 5 ) +-- > ( 2 3 4 5 1 ) +-- +cycleLeft :: Int -> Permutation +cycleLeft n = Permutation $ listArray (1,n) $ [2..n] ++ [1] + +-- | The permutation which cycles a list right by one step: +-- +-- > permuteList (cycleRight 5) "abcde" == "eabcd" +-- +-- Or in two-line notation: +-- +-- > ( 1 2 3 4 5 ) +-- > ( 5 1 2 3 4 ) +-- +cycleRight :: Int -> Permutation +cycleRight n = Permutation $ listArray (1,n) $ n : [1..n-1] + +-------------------------------------------------------------------------------- +-- * Permutation groups + +-- | Multiplies two permutations together: @p `multiply` 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. +-- multiply :: Permutation -> Permutation -> Permutation -multiply pi1@(Permutation perm1) (Permutation perm2) = +multiply pi1@(Permutation perm1) pi2@(Permutation perm2) = if (n==m) then Permutation result else error "multiply: permutations of different sets" where (_,n) = bounds perm1 (_,m) = bounds perm2 - result = permute pi1 perm2 + result = permute pi2 perm1 infixr 7 `multiply` - + -- | The inverse permutation. inverse :: Permutation -> Permutation inverse (Permutation perm1) = Permutation result @@ -287,11 +511,109 @@ result = array (1,n) $ map swap $ assocs perm1 (_,n) = bounds perm1 --- | The trivial permutation. +-- | The identity (or trivial) permutation. identity :: Int -> Permutation identity n = Permutation $ listArray (1,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 + +-- | 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)" + -------------------------------------------------------------------------------- +-- * Action of the permutation group + +-- | /Right/ action of a permutation on a set. If our permutation is +-- encoded with the sequence @[p1,p2,...,pn]@, then in the +-- two-line notation we have +-- +-- > ( 1 2 3 ... n ) +-- > ( p1 p2 p3 ... pn ) +-- +-- We adopt the convention that permutations act /on the right/ +-- (as in Knuth): +-- +-- > permute pi2 (permute pi1 set) == permute (pi1 `multiply` pi2) set +-- +-- Synonym to 'permuteRight' +-- +{-# 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 + +-- | Right action on lists. Synonym to 'permuteListRight' +-- +permuteList :: Permutation -> [a] -> [a] +permuteList = permuteRightList + +-- | The right (standard) action of permutations on sets. +-- +-- > permuteRight pi2 (permuteRight pi1 set) == permuteRight (pi1 `multiply` 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 = + if (a==1) && (b==n) + then listArray (1,n) [ ar!(perm!i) | i <- [1..n] ] + else error "permuteRight: array bounds do not match" + where + (_,n) = bounds 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] +-- +permuteRightList :: forall a . Permutation -> [a] -> [a] +permuteRightList perm xs = elems $ permuteRight 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 +-- +-- It is related to 'permuteLeft' via: +-- +-- > permuteLeft pi arr == permuteRight (inverse pi) arr +-- > permuteRight pi arr == permuteLeft (inverse 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 = + -- 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" + where + (_,n) = bounds 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] +-- +permuteLeftList :: forall a. Permutation -> [a] -> [a] +permuteLeftList perm xs = elems $ permuteLeft perm $ arr where + arr = listArray (1,n) xs :: Array Int a + n = permutationSize perm + +-------------------------------------------------------------------------------- -- * Permutations of distinct elements -- | A synonym for 'permutationsNaive' @@ -301,7 +623,7 @@ _permutations :: Int -> [[Int]] _permutations = _permutationsNaive --- | Permutations of @[1..n]@ in lexicographic order, naive algorithm. +-- | All permutations of @[1..n]@ in lexicographic order, naive algorithm. permutationsNaive :: Int -> [Permutation] permutationsNaive n = map toPermutationUnsafe $ _permutations n @@ -422,6 +744,10 @@ naturalSet perm = listArray (1,n) [ Elem i | i<-[1..n] ] where n = permutationSize perm +permInternalSet :: Permutation -> Array Int Elem +permInternalSet perm@(Permutation arr) = listArray (1,n) [ Elem (arr!i) | i<-[1..n] ] where + n = permutationSize perm + sameSize :: Permutation -> Permutation -> Bool sameSize perm1 perm2 = ( permutationSize perm1 == permutationSize perm2) @@ -467,39 +793,86 @@ checkAll = do let f :: Testable a => a -> IO () f = quickCheck - f prop_disjcyc1 - f prop_disjcyc2 + + f prop_disjcyc_1 + f prop_disjcyc_2 + + f prop_disjcyc_Mathematica + f prop_randCyclic f prop_inverse + f prop_mulPerm + f prop_mulPermLeft + f prop_mulPermRight + + f prop_perm + f prop_permLeft + f prop_permRight + f prop_permLeftRight + + f prop_cycleLeft + f prop_cycleRight + f prop_mulSign f prop_invMul f prop_cyclSign f prop_permIsPerm f prop_isEven -prop_disjcyc1 perm = ( perm == disjointCyclesToPermutation n (permutationToDisjointCycles perm) ) +prop_disjcyc_1 perm = ( perm == disjointCyclesToPermutation n (permutationToDisjointCycles perm) ) where n = permutationSize perm -prop_disjcyc2 k dcyc = ( dcyc == permutationToDisjointCycles (disjointCyclesToPermutation n dcyc) ) + +prop_disjcyc_2 k dcyc = ( dcyc == permutationToDisjointCycles (disjointCyclesToPermutation n dcyc) ) where n = fromNat k + m m = case fromDisjointCycles dcyc of - [] -> 1 + [] -> 1 xxs -> maximum (concat xxs) +-- PermutationCycles[ { 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 } ] +-- Cycles [ {{1, 12, 18, 19, 14, 21, 23, 13, 22}, {2, 15, 8, 9, 20, 16, 10, 3, 5}, {4, 6, 7, 17}} ] +prop_disjcyc_Mathematica = (permutationToDisjointCycles perm == disjcyc) + && (disjointCyclesToPermutation n disjcyc == perm) + where + n = permutationSize perm + perm = toPermutation [ 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 ] + disjcyc = DisjointCycles [ [1, 12, 18, 19, 14, 21, 23, 13, 22], [2, 15, 8, 9, 20, 16, 10, 3, 5], [4, 6, 7, 17] ] + +xperm = toPermutation [ 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 ] +xdisjcyc = DisjointCycles [ [1, 12, 18, 19, 14, 21, 23, 13, 22], [2, 15, 8, 9, 20, 16, 10, 3, 5], [4, 6, 7, 17] ] + prop_randCyclic cycl = ( isCyclicPermutation (fromCyclic cycl) ) prop_inverse perm = ( perm == inverse (inverse perm) ) prop_mulPerm (SameSize perm1 perm2) = - ( permute perm1 (permute perm2 set) == permute (perm1 `multiply` perm2) set ) + ( permute perm2 (permute perm1 set) == permute (perm1 `multiply` perm2) set ) where set = naturalSet perm1 +prop_mulPermRight (SameSize perm1 perm2) = + ( permuteRight perm2 (permuteRight perm1 set) == permuteRight (perm1 `multiply` perm2) set ) + where + set = naturalSet perm1 + +prop_mulPermLeft (SameSize perm1 perm2) = + ( permuteLeft perm2 (permuteLeft perm1 set) == permuteLeft (perm2 `multiply` perm1) set ) + where + set = naturalSet perm1 + +prop_perm perm = permute perm (naturalSet perm) == permInternalSet perm +prop_permLeft perm = permuteLeft perm (permInternalSet perm) == naturalSet perm +prop_permRight perm = permuteRight perm (naturalSet perm) == permInternalSet perm +prop_permLeftRight perm = permuteLeft (inverse perm) (naturalSet perm) == permuteRight (perm) (naturalSet perm) + +prop_cycleLeft = permuteList (cycleLeft 5) "abcde" == "bcdea" +prop_cycleRight = permuteList (cycleRight 5) "abcde" == "eabcd" + prop_mulSign (SameSize perm1 perm2) = ( sgn perm1 * sgn perm2 == sgn (perm1 `multiply` perm2) ) where - sgn = signOfPermutation :: Permutation -> Int + sgn = signValue . signOfPermutation :: Permutation -> Int prop_invMul (SameSize perm1 perm2) = ( inverse perm2 `multiply` inverse perm1 == inverse (perm1 `multiply` perm2) ) diff --git a/Math/Combinat/Sign.hs b/Math/Combinat/Sign.hs --- a/Math/Combinat/Sign.hs +++ b/Math/Combinat/Sign.hs @@ -20,11 +20,17 @@ mappend = mulSign mconcat = productOfSigns +isPlus, isMinus :: Sign -> Bool +isPlus s = case s of { Plus -> True ; _ -> False } +isMinus s = case s of { Minus -> True ; _ -> False } + +-- | @+1@ or @-1@ signValue :: Num a => Sign -> a signValue s = case s of Plus -> 1 Minus -> -1 +-- | 'Plus' if even, 'Minus' if odd paritySign :: Integral a => a -> Sign paritySign x = if even x then Plus else Minus diff --git a/Math/Combinat/Tableaux.hs b/Math/Combinat/Tableaux.hs --- a/Math/Combinat/Tableaux.hs +++ b/Math/Combinat/Tableaux.hs @@ -1,5 +1,6 @@ -- | Young tableaux and similar gadgets. +-- -- See e.g. William Fulton: Young Tableaux, with Applications to -- Representation theory and Geometry (CUP 1997). -- @@ -22,17 +23,18 @@ -- > ] -- -{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} +{-# LANGUAGE CPP, BangPatterns, FlexibleInstances, TypeSynonymInstances, MultiParamTypeClasses #-} module Math.Combinat.Tableaux where -------------------------------------------------------------------------------- import Data.List -import Math.Combinat.Helper +import Math.Combinat.Classes import Math.Combinat.Numbers (factorial,binomial) import Math.Combinat.Partitions import Math.Combinat.ASCII +import Math.Combinat.Helper import Data.Map.Strict (Map) import qualified Data.Map.Strict as Map @@ -40,28 +42,51 @@ -------------------------------------------------------------------------------- -- * Basic stuff +-- | A tableau is simply represented as a list of lists. type Tableau a = [[a]] +-- | ASCII diagram of a tableau asciiTableau :: Show a => Tableau a -> ASCII asciiTableau t = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) $ (map . map) asciiShow $ t +instance CanBeEmpty (Tableau a) where + empty = [] + isEmpty = null + instance Show a => DrawASCII (Tableau a) where ascii = asciiTableau -_shape :: Tableau a -> [Int] -_shape t = map length t +_tableauShape :: Tableau a -> [Int] +_tableauShape t = map length t -shape :: Tableau a -> Partition -shape t = toPartition (_shape t) +-- | The shape of a tableau +tableauShape :: Tableau a -> Partition +tableauShape t = toPartition (_tableauShape t) +instance HasShape (Tableau a) Partition where + shape = tableauShape + +-- | Number of entries +tableauWeight :: Tableau a -> Int +tableauWeight = sum' . map length + +instance HasWeight (Tableau a) where + weight = tableauWeight + +-- | The dual of the tableau is the mirror image to the main diagonal. dualTableau :: Tableau a -> Tableau a dualTableau = transpose -content :: Tableau a -> [a] -content = concat +instance HasDuality (Tableau a) where + dual = dualTableau +-- | The content of a tableau is the list of its entries. The ordering is from the left to the right and +-- then from the top to the bottom +tableauContent :: Tableau a -> [a] +tableauContent = concat + -- | An element @(i,j)@ of the resulting tableau (which has shape of the -- given partition) means that the vertical part of the hook has length @i@, -- and the horizontal part @j@. The /hook length/ is thus @i+j-1@. @@ -85,9 +110,12 @@ -------------------------------------------------------------------------------- -- * Row and column words +-- | The /row word/ of a tableau is the list of its entry read from the right to the left and then +-- from the top to the bottom. rowWord :: Tableau a -> [a] rowWord = concat . reverse +-- | /Semistandard/ tableaux can be reconstructed from their row words rowWordToTableau :: Ord a => [a] -> Tableau a rowWordToTableau xs = reverse rows where rows = break xs @@ -97,9 +125,11 @@ then [x] : break xs else let (h:t) = break xs in (x:h):t +-- | The /column word/ of a tableau is the list of its entry read from the bottom to the top and then from the left to the right columnWord :: Tableau a -> [a] columnWord = rowWord . transpose +-- | /Standard/ tableaux can be reconstructed from either their column or row words columnWordToTableau :: Ord a => [a] -> Tableau a columnWordToTableau = transpose . rowWordToTableau @@ -121,10 +151,48 @@ cnt j = case Map.lookup j table' of Just k -> k Nothing -> 0 - + -------------------------------------------------------------------------------- +-- * Semistandard Young tableaux + +-- | A tableau is /semistandard/ if its entries are weekly increasing horizontally +-- and strictly increasing vertically +isSemiStandardTableau :: Tableau Int -> Bool +isSemiStandardTableau t = weak && strict where + weak = and [ isWeaklyIncreasing xs | xs <- t ] + strict = and [ isStrictlyIncreasing ys | ys <- dt ] + dt = dualTableau t + +-- | Semistandard Young tableaux of given shape, \"naive\" algorithm +semiStandardYoungTableaux :: Int -> Partition -> [Tableau Int] +semiStandardYoungTableaux n part = worker (repeat 0) shape where + shape = fromPartition part + worker _ [] = [[]] + worker prevRow (s:ss) + = [ (r:rs) | r <- row n s 1 prevRow, rs <- worker (map (+1) r) ss ] + + -- weekly increasing lists of length @len@, pointwise at least @xs@, + -- maximum value @n@, minimum value @prev@. + row :: Int -> Int -> Int -> [Int] -> [[Int]] + row _ 0 _ _ = [[]] + row n len prev (x:xs) = [ (a:as) | a <- [max x prev..n] , as <- row n (len-1) a xs ] + +-- | Stanley's hook formula (cf. Fulton page 55) +countSemiStandardYoungTableaux :: Int -> Partition -> Integer +countSemiStandardYoungTableaux n shape = k `div` h where + h = product $ map fromIntegral $ concat $ hookLengths shape + k = product [ fromIntegral (n+j-i) | (i,j) <- elements shape ] + + +-------------------------------------------------------------------------------- -- * Standard Young tableaux +-- | A tableau is /standard/ if it is semistandard and its content is exactly @[1..n]@, +-- where @n@ is the weight. +isStandardTableau :: Tableau Int -> Bool +isStandardTableau t = isSemiStandardTableau t && sort (concat t) == [1..n] where + n = sum [ length xs | xs <- t ] + -- | Standard Young tableaux of a given shape. -- Adapted from John Stembridge, -- . @@ -167,29 +235,7 @@ factorial n `div` h where h = product $ map fromIntegral $ concat $ hookLengths part n = weight part - --------------------------------------------------------------------------------- --- * Semistandard Young tableaux - --- | Semistandard Young tableaux of given shape, \"naive\" algorithm -semiStandardYoungTableaux :: Int -> Partition -> [Tableau Int] -semiStandardYoungTableaux n part = worker (repeat 0) shape where - shape = fromPartition part - worker _ [] = [[]] - worker prevRow (s:ss) - = [ (r:rs) | r <- row n s 1 prevRow, rs <- worker (map (+1) r) ss ] - -- weekly increasing lists of length @len@, pointwise at least @xs@, - -- maximum value @n@, minimum value @prev@. - row :: Int -> Int -> Int -> [Int] -> [[Int]] - row _ 0 _ _ = [[]] - row n len prev (x:xs) = [ (a:as) | a <- [max x prev..n] , as <- row n (len-1) a xs ] - --- | Stanley's hook formula (cf. Fulton page 55) -countSemiStandardYoungTableaux :: Int -> Partition -> Integer -countSemiStandardYoungTableaux n shape = k `div` h where - h = product $ map fromIntegral $ concat $ hookLengths shape - k = product [ fromIntegral (n+j-i) | (i,j) <- elements shape ] - -------------------------------------------------------------------------------- + diff --git a/Math/Combinat/Tableaux/LittlewoodRichardson.hs b/Math/Combinat/Tableaux/LittlewoodRichardson.hs --- a/Math/Combinat/Tableaux/LittlewoodRichardson.hs +++ b/Math/Combinat/Tableaux/LittlewoodRichardson.hs @@ -2,8 +2,10 @@ -- | The Littlewood-Richardson rule module Math.Combinat.Tableaux.LittlewoodRichardson - ( lrRule , _lrRule - , lrRuleNaive + ( lrCoeff , lrCoeff' + , lrMult + , lrRule , _lrRule , lrRuleNaive + , lrScalar , _lrScalar ) where @@ -16,14 +18,15 @@ import Math.Combinat.Partitions.Skew import Math.Combinat.Tableaux import Math.Combinat.Tableaux.Skew +import Math.Combinat.Helper import Data.Map.Strict (Map) import qualified Data.Map.Strict as Map -------------------------------------------------------------------------------- --- | Naive, reference implementation of the Littlewood-Richardson rule, based on the definition --- "count the semistandard skew tableaux whose row content is a lattice word" +-- | Naive (very slow) reference implementation of the Littlewood-Richardson rule, based +-- on the definition "count the semistandard skew tableaux whose row content is a lattice word" -- lrRuleNaive :: SkewPartition -> Map Partition Int lrRuleNaive skew = final where @@ -33,23 +36,30 @@ f old nu = Map.insertWith (+) nu 1 old -------------------------------------------------------------------------------- +-- SKEW EXPANSION -- | @lrRule@ computes the expansion of a skew Schur function --- @s[lambda/mu]@ via the Littlewood-Richardson rule. +-- @s[lambda\/mu]@ via the Littlewood-Richardson rule. -- -- Adapted from John Stembridge's Maple code: -- -- +-- > lrRule (mkSkewPartition (lambda,nu)) == Map.fromList list where +-- > muw = weight lambda - weight nu +-- > list = [ (mu, coeff) +-- > | mu <- partitions muw +-- > , let coeff = lrCoeff lambda (mu,nu) +-- > , coeff /= 0 +-- > ] +-- lrRule :: SkewPartition -> Map Partition Int lrRule skew = _lrRule lam mu where (lam,mu) = fromSkewPartition skew -- | @_lrRule lambda mu@ computes the expansion of the skew --- Schur function @s[lambda/mu]@ via the Littlewood-Richardson rule. ---- -{-# SPECIALIZE _lrRule :: Partition -> Partition -> Map Partition Int #-} -{-# SPECIALIZE _lrRule :: Partition -> Partition -> Map Partition Integer #-} -_lrRule :: Num coeff => Partition -> Partition -> Map Partition coeff +-- Schur function @s[lambda\/mu]@ via the Littlewood-Richardson rule. +-- +_lrRule :: Partition -> Partition -> Map Partition Int _lrRule plam@(Partition lam) pmu@(Partition mu0) = if not (pmu `isSubPartitionOf` plam) then Map.empty @@ -58,7 +68,7 @@ f old nu = Map.insertWith (+) (Partition nu) 1 old diagram = [ (i,j) | (i,a,b) <- reverse (zip3 [1..] lam mu) , j <- [b+1..a] ] mu = mu0 ++ repeat 0 - n = sum $ zipWith (-) lam mu -- n == length diagram + n = sum' $ zipWith (-) lam mu -- n == length diagram {- LR_rule:=proc(lambda) local l,mu,alpha,beta,i,j,dgrm; @@ -168,5 +178,222 @@ -} -------------------------------------------------------------------------------- +-- COEFF +-- | @lrCoeff lam (mu,nu)@ computes the coressponding Littlewood-Richardson coefficients. +-- This is also the coefficient of @s[lambda]@ in the product @s[mu]*s[nu]@ +-- +-- Note: This is much slower than using 'lrRule' or 'lrMult' to compute several coefficients +-- at the same time! +lrCoeff :: Partition -> (Partition,Partition) -> Int +lrCoeff lam (mu,nu) = + if nu `isSubPartitionOf` lam + then lrScalar (mkSkewPartition (lam,nu)) (mkSkewPartition (mu,emptyPartition)) + else 0 +-- | @lrCoeff (lam\/nu) mu@ computes the coressponding Littlewood-Richardson coefficients. +-- This is also the coefficient of @s[mu]@ in the product @s[lam\/nu]@ +-- +-- Note: This is much slower than using 'lrRule' or 'lrMult' to compute several coefficients +-- at the same time! +lrCoeff' :: SkewPartition -> Partition -> Int +lrCoeff' skew p = lrScalar skew (mkSkewPartition (p,emptyPartition)) + +-------------------------------------------------------------------------------- +-- SCALAR PRODUCT + +-- | @lrScalar (lambda\/mu) (alpha\/beta)@ computes the scalar product of the two skew +-- Schur functions @s[lambda\/mu]@ and @s[alpha\/beta]@ via the Littlewood-Richardson rule. +-- +-- Adapted from John Stembridge Maple code: +-- +-- +lrScalar :: SkewPartition -> SkewPartition -> Int +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)) = + if not (pmu `isSubPartitionOf` plam ) + || not (pbeta `isSubPartitionOf` palpha) + || (sum' lam + sum' beta) /= (sum' alpha + sum' mu0) -- equivalent to (lambda-mu) /= (alpha-beta) + then 0 + else length $ fillings' n diagram (alpha,beta) + where + f old nu = Map.insertWith (+) (Partition nu) 1 old + diagram = [ (i,j) | (i,a,b) <- reverse (zip3 [1..] lam mu) , j <- [b+1..a] ] + mu = mu0 ++ repeat 0 + n = sum' $ zipWith (-) lam mu -- n == length diagram + +{- +LR_rule:=proc(lambda) local l,mu,alpha,beta,i,j,dgrm; + if not `LR_rule/fit`(lambda,args[2]) then RETURN(0) fi; + l:=nops(lambda); mu:=[op(args[2]),0$l]; + dgrm:=[seq(seq([i,-j],j=-lambda[i]..-1-mu[i]),i=1..l)]; + if nargs>2 then alpha:=args[3]; + if nargs>3 then beta:=args[4] else beta:=[] fi; + if not `LR_rule/fit`(alpha,beta) then RETURN(0) fi; + l:=convert([op(lambda),op(beta)],`+`); + if l<>convert([op(alpha),op(mu)],`+`) then RETURN(0) fi; + nops(LR_fillings(dgrm,[alpha,beta])) + else + convert([seq(s[op(i[1])],i=LR_fillings(dgrm))],`+`) + fi +end; +-} + +-------------------------------------------------------------------------------- + +-- | Note: we use reverse ordering in Diagrams compared the Stembridge's code. +-- Also, for performance reasons, we need the length of the diagram +fillings' :: Int -> Diagram -> ([Int],[Int]) -> [Filling] +fillings' _ [] (alpha,beta) = [ (beta,[]) ] +fillings' n diagram@((x,y):rest) alphaBeta@(alpha,beta) = stuff where + stuff = concatMap (nextLetter' lower upper alpha) (fillings' (n-1) rest alphaBeta) + upper = case findIndex (==(x ,y+1)) diagram of { Just j -> n-j ; Nothing -> 0 } + lower = case findIndex (==(x-1,y )) diagram of { Just j -> n-j ; Nothing -> 0 } + +{- +LR_fillings:=proc(dgrm) local n,x,upper,lower; + if dgrm=[] then + if nargs=1 then x:=[] else x:=args[2][2] fi; + RETURN([[x,[]]]) + fi; + n:=nops(dgrm); x:=dgrm[n]; + if not member([x[1],x[2]+1],dgrm,'upper') then upper:=0 fi; + if not member([x[1]-1,x[2]],dgrm,'lower') then lower:=0 fi; + if nargs=1 then + map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)]),lower,upper) + else + map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)],args[2]), + lower,upper,args[2][1]) + fi; +end: +-} + +-------------------------------------------------------------------------------- + +nextLetter' :: Int -> Int -> [Int] -> Filling -> [Filling] +nextLetter' lower upper alpha (nu,lpart) = stuff where + stuff = [ ( incr i shape , lpart++[i] ) | i<-nlist ] + shape = nu ++ [0] + lb = if lower>0 + then lpart !! (lower-1) + else 0 + ub1 = if upper>0 + then min (length shape) (lpart !! (upper-1)) + else length shape + ub = min (length alpha) ub1 + nlist = filter (>0) $ map f [lb+1..ub] + f j = if ( shape!!(j-1) < alpha!!(j-1)) && + (j==1 || shape!!(j-2) > shape!!(j-1)) + then j + else 0 + + -- increments the i-th element (starting from 1) + incr :: Int -> [Int] -> [Int] + incr i (x:xs) = case i of + 0 -> finish (x:xs) + 1 -> (x+1) : finish xs + _ -> x : incr (i-1) xs + incr _ [] = [] + + -- removes tailing zeros + finish :: [Int] -> [Int] + finish (x:xs) = if x>0 then x : finish xs else [] + finish [] = [] + +{- +`LR/nextletter`:=proc(T) local shape,lp,lb,ub,i,nl; + shape:=[op(T[1]),0]; lp:=T[2]; ub:=nops(shape); + if nargs>3 then ub:=min(ub,nops(args[4])) fi; + if args[2]=0 then lb:=0 else lb:=lp[args[2]] fi; + if args[3]>0 then ub:=min(lp[args[3]],ub) fi; + if nargs<4 then + nl:=map(proc(x,y) if x=1 or y[x-1]>y[x] then x fi end,[$lb+1..ub],shape) + else + nl:=map(proc(x,y,z) if y[x]y[x]) then x fi end, + [$lb+1..ub],shape,args[4]) + fi; + nl:=[seq([subsop(i=shape[i]+1,shape),[op(lp),i]],i=nl)]; + op(subs(0=NULL,nl)) +end: +-} + +-------------------------------------------------------------------------------- +-- MULTIPLICATION + +type Part = [Int] + +-- | Computes the expansion of the product of Schur polynomials @s[mu]*s[nu]@ using the +-- Littlewood-Richardson rule. Note: this is symmetric in the two arguments. +-- +-- Based on the wikipedia article +-- +-- > lrMult mu nu == Map.fromList list where +-- > lamw = weight nu + weight mu +-- > list = [ (lambda, coeff) +-- > | lambda <- partitions lamw +-- > , let coeff = lrCoeff lambda (mu,nu) +-- > , coeff /= 0 +-- > ] +-- +lrMult :: Partition -> Partition -> Map Partition Int +lrMult pmu@(Partition mu) pnu@(Partition nu) = result where + result = foldl' add Map.empty (addMu mu nu) where + add !old lambda = Map.insertWith (+) (Partition lambda) 1 old + +-- | This basically lists all the outer shapes (with multiplicities) which can be result from the LR rule +addMu :: Part -> Part -> [Part] +addMu mu part = go ubs0 mu dmu part where + + go _ [] _ part = [part] + go ubs (m:ms) (d:ds) part = concat [ go (drop d ubs') ms ds part' | (ubs',part') <- addRowOf ubs part ] + + ubs0 = take (headOrZero mu) [headOrZero part + 1..] + dmu = diffSeq mu + + +-- | Adds a full row of @(length pcols)@ boxes containing to a partition, with +-- pcols being the upper bounds of the columns, respectively. We also return the +-- newly added columns +addRowOf :: [Int] -> Part -> [([Int],Part)] +addRowOf pcols part = go 0 pcols part [] where + go !lb [] p ncols = [(reverse ncols , p)] + go !lb (!ub:ubs) p ncols = concat [ go col ubs (addBox ij p) (col:ncols) | ij@(row,col) <- newBoxes (lb+1) ub p ] + +-- | Returns the (row,column) pairs of the new boxes which +-- can be added to the given partition with the given column bounds +-- and the 1-Rieri rule +newBoxes :: Int -> Int -> Part -> [(Int,Int)] +newBoxes lb ub part = reverse $ go [1..] part (headOrZero part + 1) where + go (!i:_ ) [] !lp + | lb <= 1 && 1 <= ub && lp > 0 = [(i,1)] + | otherwise = [] + go (!i:is) (!j:js) !lp + | j1 < lb = [] + | j1 <= ub && lp > j = (i,j1) : go is js j + | otherwise = go is js j + where + j1 = j+1 + +-- | Adds a box to a partition +addBox :: (Int,Int) -> Part -> Part +addBox (k,_) part = go 1 part where + go !i (p:ps) = if i==k then (p+1):ps else p : go (i+1) ps + go !i [] = if i==k then [1] else error "addBox: shouldn't happen" + +-- | Safe head defaulting to zero +headOrZero :: [Int] -> Int +headOrZero xs = case xs of + (!x:_) -> x + [] -> 0 + +-- | Computes the sequence of differences from a partition (including the last difference to zero) +diffSeq :: Part -> [Int] +diffSeq = go where + go (p:ps@(q:_)) = (p-q) : go ps + go [p] = [p] + go [] = [] + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Tableaux/Skew.hs b/Math/Combinat/Tableaux/Skew.hs --- a/Math/Combinat/Tableaux/Skew.hs +++ b/Math/Combinat/Tableaux/Skew.hs @@ -1,7 +1,12 @@ -- | Skew tableaux are skew partitions filled with numbers. +-- +-- For example: +-- +-- <> +-- -{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} +{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, MultiParamTypeClasses #-} module Math.Combinat.Tableaux.Skew where @@ -9,16 +14,18 @@ import Data.List +import Math.Combinat.Classes import Math.Combinat.Partitions.Integer import Math.Combinat.Partitions.Skew import Math.Combinat.Tableaux import Math.Combinat.ASCII +import Math.Combinat.Helper import Data.Map.Strict (Map) import qualified Data.Map.Strict as Map -------------------------------------------------------------------------------- - +-- * Basics -- | A skew tableau is represented by a list of offsets and entries newtype SkewTableau a = SkewTableau [(Int,[a])] deriving (Eq,Ord,Show) @@ -27,12 +34,24 @@ instance Functor SkewTableau where fmap f (SkewTableau t) = SkewTableau [ (a, map f xs) | (a,xs) <- t ] - -skewShape :: SkewTableau a -> SkewPartition -skewShape (SkewTableau list) = SkewPartition [ (o,length xs) | (o,xs) <- list ] +-- | The shape of a skew tableau +skewTableauShape :: SkewTableau a -> SkewPartition +skewTableauShape (SkewTableau list) = SkewPartition [ (o,length xs) | (o,xs) <- list ] + +instance HasShape (SkewTableau a) SkewPartition where + shape = skewTableauShape + +-- | The weight of a tableau is the weight of its shape, or the number of entries +skewTableauWeight :: SkewTableau a -> Int +skewTableauWeight = skewPartitionWeight . skewTableauShape + +instance HasWeight (SkewTableau a) where + weight = skewTableauWeight + -------------------------------------------------------------------------------- +-- | The dual of a skew tableau, that is, its mirror image to the main diagonal dualSkewTableau :: forall a. SkewTableau a -> SkewTableau a dualSkewTableau (SkewTableau axs) = SkewTableau (go axs) where @@ -71,9 +90,29 @@ ] -} +instance HasDuality (SkewTableau a) where + dual = dualSkewTableau + -------------------------------------------------------------------------------- +-- * Semistandard tableau --- | Semi-standard skew tableaux filled with numbers @[1..n]@ +-- | A tableau is /semistandard/ if its entries are weekly increasing horizontally +-- and strictly increasing vertically +isSemiStandardSkewTableau :: SkewTableau Int -> Bool +isSemiStandardSkewTableau st@(SkewTableau axs) = weak && strict where + weak = and [ isWeaklyIncreasing xs | (a,xs) <- axs ] + strict = and [ isStrictlyIncreasing ys | (b,ys) <- bys ] + SkewTableau bys = dualSkewTableau st + +-- | A tableau is /standard/ if it is semistandard and its content is exactly @[1..n]@, +-- where @n@ is the weight. +isStandardSkewTableau :: SkewTableau Int -> Bool +isStandardSkewTableau st = isSemiStandardSkewTableau st && sort (skewTableauRowWord st) == [1..n] where + n = skewTableauWeight st + +-------------------------------------------------------------------------------- + +-- | All semi-standard skew tableaux filled with the numbers @[1..n]@ semiStandardSkewTableaux :: Int -> SkewPartition -> [SkewTableau Int] semiStandardSkewTableaux n (SkewPartition abs) = map SkewTableau stuff where @@ -105,14 +144,16 @@ -} -------------------------------------------------------------------------------- +-- * ASCII +-- | ASCII drawing of a skew tableau (using the English notation) asciiSkewTableau :: Show a => SkewTableau a -> ASCII asciiSkewTableau = asciiSkewTableau' "." EnglishNotation asciiSkewTableau' :: Show a - => String -- ^ string representing the elements of the inner (unfilled) partition - -> PartitionConvention -- Orientation + => String -- ^ string representing the elements of the inner (unfilled) partition + -> PartitionConvention -- ^ orientation -> SkewTableau a -> ASCII asciiSkewTableau' innerstr orient (SkewTableau axs) = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) stuff where @@ -127,12 +168,13 @@ ascii = asciiSkewTableau -------------------------------------------------------------------------------- +-- * Row \/ column words --- | The reversed rows, concatenated +-- | The reversed (right-to-left) rows, concatenated skewTableauRowWord :: SkewTableau a -> [a] skewTableauRowWord (SkewTableau axs) = concatMap (reverse . snd) axs --- | The reversed rows, concatenated +-- | The reversed (bottom-to-top) columns, concatenated skewTableauColumnWord :: SkewTableau a -> [a] skewTableauColumnWord = skewTableauRowWord . dualSkewTableau @@ -176,5 +218,38 @@ cnt j = case Map.lookup j table' of Just k -> k Nothing -> 0 + +-------------------------------------------------------------------------------- + +#ifdef QUICKCHECK + +prop_dual_dual :: SkewTableau Int -> Bool +prop_dual_dual st = (dualSkewTableau (dualSkewTableau st) == st) + +prop_rowWord :: SkewTableau Int -> Bool +prop_rowWord st = (fillSkewPartitionWithRowWord shape content == st) where + shape = skewShape st + content = skewTableauRowWord st + +prop_columnWord :: SkewTableau Int -> Bool +prop_columnWord st = (fillSkewPartitionWithColumnWord shape content == st) where + shape = skewShape st + content = skewTableauColumnWord st + +prop_fill_shape :: SkewPartition -> Bool +prop_fill_shape shape = (shape == shape') where + tableau = fillSkewPartitionWithColumnWord shape [1..] + shape' = skewShape tableau + +prop_semistandard :: SkewPartition -> Bool +prop_semistandard shape = and + [ isSemiStandardSkewTableau st + | n <- [1..nn] + , st <- semiStandardSkewTableaux n shape + ] + where + nn = skewPartitionWeight shape + +#endif -------------------------------------------------------------------------------- diff --git a/Math/Combinat/Trees/Binary.hs b/Math/Combinat/Trees/Binary.hs --- a/Math/Combinat/Trees/Binary.hs +++ b/Math/Combinat/Trees/Binary.hs @@ -11,15 +11,22 @@ module Math.Combinat.Trees.Binary ( -- * Types BinTree(..) - , leaf + , leaf + , graft , BinTree'(..) , forgetNodeDecorations , Paren(..) , parenthesesToString - , stringToParentheses - -- * Conversion to rose trees (@Data.Tree@) + , stringToParentheses + , numberOfNodes + , numberOfLeaves + -- * Conversion to rose trees (@Data.Tree@) , toRoseTree , toRoseTree' , module Data.Tree + -- * Enumerate leaves + , enumerateLeaves_ + , enumerateLeaves + , enumerateLeaves' -- * Nested parentheses , nestedParentheses , randomNestedParentheses @@ -78,9 +85,9 @@ import Math.Combinat.Trees.Graphviz ( Dot , graphvizDotBinTree , graphvizDotBinTree' - , graphvizDotForest , graphvizDotTree + , graphvizDotForest , graphvizDotTree ) - +import Math.Combinat.Classes import Math.Combinat.Helper import Math.Combinat.ASCII as ASCII @@ -96,6 +103,14 @@ leaf :: BinTree () leaf = Leaf () +-- | The monadic join operation of binary trees +graft :: BinTree (BinTree a) -> BinTree a +graft = go where + go (Branch l r) = Branch (go l) (go r) + go (Leaf t ) = t + +-------------------------------------------------------------------------------- + -- | A binary tree with leaves and internal nodes decorated -- with types @a@ and @b@, respectively. data BinTree' a b @@ -104,11 +119,59 @@ deriving (Eq,Ord,Show,Read) forgetNodeDecorations :: BinTree' a b -> BinTree a -forgetNodeDecorations (Branch' left _ right) = - Branch (forgetNodeDecorations left) (forgetNodeDecorations right) -forgetNodeDecorations (Leaf' decor) = Leaf decor +forgetNodeDecorations = go where + go (Branch' left _ right) = Branch (go left) (go right) + go (Leaf' decor ) = Leaf decor -------------------------------------------------------------------------------- + +instance HasNumberOfNodes (BinTree a) where + numberOfNodes = go where + go (Leaf _ ) = 0 + go (Branch l r) = go l + go r + 1 + +instance HasNumberOfLeaves (BinTree a) where + numberOfLeaves = go where + go (Leaf _ ) = 1 + go (Branch l r) = go l + go r + + +instance HasNumberOfNodes (BinTree' a b) where + numberOfNodes = go where + go (Leaf' _ ) = 0 + go (Branch' l _ r) = go l + go r + 1 + +instance HasNumberOfLeaves (BinTree' a b) where + numberOfLeaves = go where + go (Leaf' _ ) = 1 + go (Branch' l _ r) = go l + go r + +-------------------------------------------------------------------------------- +-- * Enumerate leaves + +-- | Enumerates the leaves a tree, starting from 0, ignoring old labels +enumerateLeaves_ :: BinTree a -> BinTree Int +enumerateLeaves_ = snd . go 0 where + go !k t = case t of + Leaf _ -> (k+1 , Leaf k) + Branch l r -> (k'', Branch l' r') where + (k' ,l') = go k l + (k'',r') = go k' r + +-- | Enumerates the leaves a tree, starting from zero, and also returns the number of leaves +enumerateLeaves' :: BinTree a -> (Int, BinTree (a,Int)) +enumerateLeaves' = go 0 where + go !k t = case t of + Leaf y -> (k+1 , Leaf (y,k)) + Branch l r -> (k'', Branch l' r') where + (k' ,l') = go k l + (k'',r') = go k' r + +-- | Enumerates the leaves a tree, starting from zero +enumerateLeaves :: BinTree a -> BinTree (a,Int) +enumerateLeaves = snd . enumerateLeaves' + +-------------------------------------------------------------------------------- -- * conversion to 'Data.Tree' -- | Convert a binary tree to a rose tree (from "Data.Tree") @@ -140,6 +203,18 @@ go (Leaf x) = Leaf <$> f x go (Branch left right) = Branch <$> go left <*> go right +instance Applicative BinTree where + pure = Leaf + u <*> t = go u where + go (Branch l r) = Branch (go l) (go r) + go (Leaf f ) = fmap f t + +instance Monad BinTree where + return = Leaf + (>>=) t f = go t where + go (Branch l r) = Branch (go l) (go r) + go (Leaf y ) = f y + -------------------------------------------------------------------------------- -- * Nested parentheses @@ -324,7 +399,7 @@ -------------------------------------------------------------------------------- -- * Generating binary trees --- | Generates all binary trees with n nodes. +-- | Generates all binary trees with @n@ nodes. -- At the moment just a synonym for 'binaryTreesNaive'. binaryTrees :: Int -> [BinTree ()] binaryTrees = binaryTreesNaive diff --git a/Math/Combinat/Trees/Nary.hs b/Math/Combinat/Trees/Nary.hs --- a/Math/Combinat/Trees/Nary.hs +++ b/Math/Combinat/Trees/Nary.hs @@ -4,8 +4,11 @@ {-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} module Math.Combinat.Trees.Nary ( - -- * Regular trees - ternaryTrees + -- * Types + module Data.Tree + , Tree(..) + -- * Regular trees + , ternaryTrees , regularNaryTrees , semiRegularTrees , countTernaryTrees @@ -71,8 +74,23 @@ import Math.Combinat.Trees.Graphviz ( Dot , graphvizDotForest , graphvizDotTree ) +import Math.Combinat.Classes import Math.Combinat.ASCII as ASCII import Math.Combinat.Helper + +-------------------------------------------------------------------------------- + +instance HasNumberOfNodes (Tree a) where + numberOfNodes = go where + go (Node label subforest) = if null subforest + then 0 + else 1 + sum' (map go subforest) + +instance HasNumberOfLeaves (Tree a) where + numberOfLeaves = go where + go (Node label subforest) = if null subforest + then 1 + else sum' (map go subforest) -------------------------------------------------------------------------------- diff --git a/combinat.cabal b/combinat.cabal --- a/combinat.cabal +++ b/combinat.cabal @@ -1,10 +1,10 @@ Name: combinat -Version: 0.2.7.2 +Version: 0.2.8.0 Synopsis: Generate and manipulate various combinatorial objects. -Description: A collection of functions to generate, count and manipulate - all kinds of combinatorial objects like partitions, - compositions, permutations, Young tableaux, binary trees, - and so on. +Description: A collection of functions to generate, count, manipulate + and visualize all kinds of combinatorial objects like + partitions, compositions, trees, permutations, braids, + Young tableaux, and so on. License: BSD3 License-file: LICENSE Author: Balazs Komuves @@ -13,7 +13,7 @@ Homepage: http://code.haskell.org/~bkomuves/ Stability: Experimental Category: Math -Tested-With: GHC == 7.8.3 +Tested-With: GHC == 7.10.2 Cabal-Version: >= 1.18 Build-Type: Simple @@ -25,23 +25,17 @@ Flag withQuickCheck Description: Compile with the QuickCheck tests. default: False - -Flag base4 - Description: Base v4 Library - if flag(base4) - Build-Depends: base >= 4 && < 5, array >= 0.4, containers, random, transformers - cpp-options: -DBASE_VERSION=4 - else - Build-Depends: base >= 3 && < 4, array >= 0.4, containers, random, transformers - cpp-options: -DBASE_VERSION=3 + Build-Depends: base >= 4 && < 5, array >= 0.5, containers, random, transformers if flag(withQuickCheck) Build-Depends: QuickCheck + cpp-options: -DQUICKCHECK Exposed-Modules: Math.Combinat + Math.Combinat.Classes Math.Combinat.Numbers Math.Combinat.Numbers.Series Math.Combinat.Numbers.Primes @@ -49,6 +43,10 @@ Math.Combinat.Sets Math.Combinat.Tuples Math.Combinat.Compositions + Math.Combinat.Groups.Thompson.F + Math.Combinat.Groups.Free + Math.Combinat.Groups.Braid + Math.Combinat.Groups.Braid.NF Math.Combinat.Partitions Math.Combinat.Partitions.Integer Math.Combinat.Partitions.Skew @@ -68,20 +66,17 @@ Math.Combinat.Trees.Nary Math.Combinat.Trees.Graphviz Math.Combinat.LatticePaths - Math.Combinat.FreeGroups Math.Combinat.ASCII Math.Combinat.Helper Default-Extensions: CPP, BangPatterns Other-Extensions: MultiParamTypeClasses, ScopedTypeVariables, - GeneralizedNewtypeDeriving, BangPatterns + GeneralizedNewtypeDeriving, + DataKinds, KindSignatures Default-Language: Haskell2010 Hs-Source-Dirs: . - if flag(withQuickCheck) - cpp-options: -DQUICKCHECK - - ghc-options: -Wall -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports + ghc-options: -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports diff --git a/svg/dyck_path.svg b/svg/dyck_path.svg --- a/svg/dyck_path.svg +++ b/svg/dyck_path.svg @@ -1,4 +1,3 @@ - + "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"> diff --git a/svg/skew3.svg b/svg/skew3.svg new file mode 100644 --- /dev/null +++ b/svg/skew3.svg @@ -0,0 +1,3 @@ + + diff --git a/svg/skew_tableau.svg b/svg/skew_tableau.svg new file mode 100644 --- /dev/null +++ b/svg/skew_tableau.svg @@ -0,0 +1,3 @@ + +7512122111111 diff --git a/svg/src/gen_figures.hs b/svg/src/gen_figures.hs --- a/svg/src/gen_figures.hs +++ b/svg/src/gen_figures.hs @@ -9,14 +9,18 @@ import Math.Combinat.Partitions.Integer import Math.Combinat.Partitions.Plane import Math.Combinat.Partitions.NonCrossing +import Math.Combinat.Partitions.Skew import Math.Combinat.Tableaux +import Math.Combinat.Tableaux.Skew import Math.Combinat.LatticePaths import Math.Combinat.Trees.Binary import Math.Combinat.Diagrams.Partitions.Integer import Math.Combinat.Diagrams.Partitions.Plane import Math.Combinat.Diagrams.Partitions.NonCrossing +import Math.Combinat.Diagrams.Partitions.Skew import Math.Combinat.Diagrams.Tableaux +import Math.Combinat.Diagrams.Tableaux.Skew import Math.Combinat.Diagrams.LatticePaths import Math.Combinat.Diagrams.Trees.Binary @@ -36,17 +40,20 @@ go [] = [] go zs = take m zs : go (drop m zs) +padding fac diag = pad fac $ centerXY diag +margin siz diag = hcat [ strutX siz , vcat [ strutY siz , centerXY diag , strutY siz ] , strutX siz ] + -------------------------------------------------------------------------------- main = do - export "plane_partition.svg" (Width 320) $ drawPlanePartition3D $ + export "plane_partition.svg" (mkWidth 320) $ margin 0.05 $ drawPlanePartition3D $ PlanePart [[5,4,3,3,1],[4,4,2,1],[3,2],[2,1],[1],[1]] - export "noncrossing.svg" (Width 256) $ pad 1.10 $ drawNonCrossingCircleDiagram' orange True $ + export "noncrossing.svg" (mkWidth 256) $ padding 1.10 $ drawNonCrossingCircleDiagram' orange True $ NonCrossing [[3],[5,4,2],[7,6,1],[9,8]] - export "young_tableau.svg" (Width 256) $ drawTableau $ + export "young_tableau.svg" (mkWidth 256) $ margin 0.05 $ drawTableau $ [ [ 1 , 3 , 4 , 6 , 7 ] , [ 2 , 5 , 8 ,10 ] , [ 9 ] @@ -55,12 +62,20 @@ let u = UpStep d = DownStep path = [ u,u,d,u,u,u,d,u,d,d,u,d,u,u,u,d,d,d,d,d,u,d,u,u,d,d ] - export "dyck_path.svg" (Width 500) $ drawLatticePath $ path + export "dyck_path.svg" (mkWidth 500) $ margin 0.05 $ drawLatticePath $ path -- print (pathHeight path, pathNumberOfZeroTouches path, pathNumberOfPeaks path) - export "ferrers.svg" (Width 256) $ drawFerrersDiagram' EnglishNotation red True $ + export "ferrers.svg" (mkWidth 256) $ margin 0.05 $ drawFerrersDiagram' EnglishNotation red True $ Partition [8,6,3,3,1] - export "bintrees.svg" (Width 750) $ boxSep 7 $ map drawBinTree_ (binaryTrees 4) + export "bintrees.svg" (mkWidth 750) $ boxSep 7 $ map drawBinTree_ (binaryTrees 4) + + let skew = mkSkewPartition (Partition [9,7,3,2,2,1] , Partition [5,3,2,1]) + -- export "skew.svg" (mkWidth 256) $ margin 0.05 $ drawSkewFerrersDiagram skew + -- export "skew2.svg" (mkWidth 256) $ margin 0.05 $ drawSkewFerrersDiagram' EnglishNotation green True (True,True) skew + export "skew3.svg" (mkWidth 256) $ margin 0.05 $ drawSkewPartitionBoxes EnglishNotation skew + + let skewtableau = (semiStandardSkewTableaux 7 skew) !! 123 + export "skew_tableau.svg" (mkWidth 320) $ margin 0.05 $ drawSkewTableau' EnglishNotation blue True skewtableau --------------------------------------------------------------------------------