combinat 0.2.7.2 → 0.2.8.0
raw patch · 29 files changed
+3710/−645 lines, 29 filesdep ~arraydep ~base
Dependency ranges changed: array, base
Files
- Math/Combinat.hs +15/−10
- Math/Combinat/ASCII.hs +200/−62
- Math/Combinat/Classes.hs +66/−0
- Math/Combinat/FreeGroups.hs +0/−376
- Math/Combinat/Groups/Braid.hs +651/−0
- Math/Combinat/Groups/Braid/NF.hs +544/−0
- Math/Combinat/Groups/Free.hs +523/−0
- Math/Combinat/Groups/Thompson/F.hs +404/−0
- Math/Combinat/Helper.hs +57/−2
- Math/Combinat/LatticePaths.hs +7/−0
- Math/Combinat/Numbers/Primes.hs +64/−0
- Math/Combinat/Numbers/Series.hs +3/−3
- Math/Combinat/Partitions/Integer.hs +70/−17
- Math/Combinat/Partitions/NonCrossing.hs +1/−1
- Math/Combinat/Partitions/Plane.hs +9/−1
- Math/Combinat/Partitions/Set.hs +1/−1
- Math/Combinat/Partitions/Skew.hs +91/−3
- Math/Combinat/Permutations.hs +449/−76
- Math/Combinat/Sign.hs +6/−0
- Math/Combinat/Tableaux.hs +78/−32
- Math/Combinat/Tableaux/LittlewoodRichardson.hs +238/−11
- Math/Combinat/Tableaux/Skew.hs +85/−10
- Math/Combinat/Trees/Binary.hs +84/−9
- Math/Combinat/Trees/Nary.hs +20/−2
- combinat.cabal +16/−21
- svg/dyck_path.svg +1/−2
- svg/skew3.svg +3/−0
- svg/skew_tableau.svg +3/−0
- svg/src/gen_figures.hs +21/−6
Math/Combinat.hs view
@@ -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
Math/Combinat/ASCII.hs view
@@ -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
+ Math/Combinat/Classes.hs view
@@ -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++--------------------------------------------------------------------------------
− Math/Combinat/FreeGroups.hs
@@ -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- ----------------------------------------------------------------------------------
+ Math/Combinat/Groups/Braid.hs view
@@ -0,0 +1,651 @@++-- | Braids. See eg. <https://en.wikipedia.org/wiki/Braid_group>+--+--+-- Based on: +--+-- * Joan S. Birman, Tara E. Brendle: BRAIDS - A SURVEY+-- <https://www.math.columbia.edu/~jb/Handbook-21.pdf>+--+--+-- 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<t)is the generator @sigma_{s,t}@ in Birman-Ko-Lee's+-- \"new presentation\". It twistes the strands @s@ and @t@ while going over all+-- other strands. For @t==s+1@ we get back @sigma s@+-- +doubleSigma :: KnownNat n => 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<n then Sigma i else error "positiveWord: index out of range"+ +-- | The (positive) half-twist of all the braid strands, usually denoted by @Delta@.+halfTwist :: KnownNat n => 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++--------------------------------------------------------------------------------
+ Math/Combinat/Groups/Braid/NF.hs view
@@ -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+-- <https://www.math.columbia.edu/~jb/Handbook-21.pdf> (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<n = (P.reversePermutation n `P.multiply` P.adjacentTransposition n k)+ | otherwise = error "posLPerm: index out of range"++-- | The permutation @posR k :: Braid n@ is realizing+posRPerm :: Int -> Int -> Permutation+posRPerm n k + | k>0 && k<n = (P.adjacentTransposition n k `P.multiply` P.reversePermutation n )+ | otherwise = error "posRPerm: index out of range"++--------------------------------------------------------------------------------++-- | We recognize left-greedy factors which are @Delta@-s (easy, since they are the only ones+-- with length @(n choose 2)@), and move them to the left, returning their summed exponent+-- and the filtered new factors. We also filter trivial permutations (which should only happen +-- for the trivial braid, but it happens there?)+--+filterDeltaFactors :: Int -> [[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<a) -- Lemma 2.4 in [2]++-- | This satisfies+-- +-- > 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++--------------------------------------------------------------------------------
+ Math/Combinat/Groups/Free.hs view
@@ -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 ]++-}++--------------------------------------------------------------------------------++
+ Math/Combinat/Groups/Thompson/F.hs view
@@ -0,0 +1,404 @@++-- | Thompson's group F.+--+-- See eg. <https://en.wikipedia.org/wiki/Thompson_groups>+--+-- Based mainly on James Michael Belk's PhD thesis \"THOMPSON'S GROUP F\";+-- see <http://www.math.u-psud.fr/~breuilla/Belk.pdf>+--++{-# 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]++--------------------------------------------------------------------------------++
Math/Combinat/Helper.hs view
@@ -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
Math/Combinat/LatticePaths.hs view
@@ -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)
Math/Combinat/Numbers/Primes.hs view
@@ -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) --------------------------------------------------------------------------------
Math/Combinat/Numbers/Series.hs view
@@ -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):
Math/Combinat/Partitions/Integer.hs view
@@ -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 --------------------------------------------------------------------------------
Math/Combinat/Partitions/NonCrossing.hs view
@@ -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
Math/Combinat/Partitions/Plane.hs view
@@ -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
Math/Combinat/Partitions/Set.hs view
@@ -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
Math/Combinat/Partitions/Skew.hs view
@@ -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+--+-- <<svg/skew3.svg>>+-- -{-# 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++--------------------------------------------------------------------------------
Math/Combinat/Permutations.hs view
@@ -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<n = transposition n (k,k+1)+ | otherwise = error "adjacentTransposition: index out of range"++-- | Product of adjacent transpositions.+--+-- > 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) )
Math/Combinat/Sign.hs view
@@ -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
Math/Combinat/Tableaux.hs view
@@ -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, -- <http://www.math.lsa.umich.edu/~jrs/software/SFexamples/tableaux>.@@ -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 ]- --------------------------------------------------------------------------------+
Math/Combinat/Tableaux/LittlewoodRichardson.hs view
@@ -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: -- <http://www.math.lsa.umich.edu/~jrs/software/SFexamples/LR_rule> --+-- > 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: +-- <http://www.math.lsa.umich.edu/~jrs/software/SFexamples/LR_rule>+--+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]<z[x] and (x=1 or y[x-1]>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 <https://en.wikipedia.org/wiki/LittlewoodRichardson_rule>+--+-- > 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 [] = []++--------------------------------------------------------------------------------
Math/Combinat/Tableaux/Skew.hs view
@@ -1,7 +1,12 @@ -- | Skew tableaux are skew partitions filled with numbers.+--+-- For example:+--+-- <<svg/skew_tableau.svg>>+-- -{-# 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 --------------------------------------------------------------------------------
Math/Combinat/Trees/Binary.hs view
@@ -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
Math/Combinat/Trees/Nary.hs view
@@ -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) --------------------------------------------------------------------------------
combinat.cabal view
@@ -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
svg/dyck_path.svg view
@@ -1,4 +1,3 @@ <?xml version="1.0" encoding="UTF-8"?> <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"- "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">-<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" version="1.1" width="500.00000000000006" height="96.15384615384616" font-size="1" viewBox="0 0 500 96" stroke="rgb(0,0,0)" stroke-opacity="1"><g><g fill="rgb(0,0,0)" fill-opacity="0.0"><g stroke="rgb(0,0,0)" stroke-opacity="1.0" stroke-width="0.8770580193070295" stroke-linecap="butt" stroke-linejoin="miter" font-size="1.0em" stroke-miterlimit="10.0"><g stroke="rgb(128,128,128)" stroke-opacity="1.0" stroke-width="0.4370629370629371"><g><path d="M 454.54545454545456,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 437.06293706293707,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 419.58041958041963,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 402.09790209790214,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 384.61538461538464,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 367.13286713286715,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 349.65034965034965,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 332.16783216783216,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 314.6853146853147,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 297.2027972027972,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 279.72027972027973,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 262.23776223776224,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 244.75524475524477,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 227.27272727272728,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 209.79020979020981,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 192.30769230769232,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 174.82517482517483,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 157.34265734265736,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 139.86013986013987,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 122.37762237762239,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 104.89510489510491,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 87.41258741258741,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 69.93006993006993,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 52.447552447552454,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 34.96503496503497,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 17.482517482517483,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 0.0,96.15384615384616 v -87.41258741258741 " /></g><g><path d="M 0.0,8.741258741258747 h 454.54545454545456 " /></g><g><path d="M 0.0,26.223776223776227 h 454.54545454545456 " /></g><g><path d="M 0.0,43.70629370629371 h 454.54545454545456 " /></g><g><path d="M 0.0,61.188811188811194 h 454.54545454545456 " /></g><g><path d="M 0.0,78.67132867132868 h 454.54545454545456 " /></g><g><path d="M 0.0,96.15384615384616 h 454.54545454545456 " /></g></g><g stroke="rgb(255,0,0)" stroke-opacity="1.0" stroke-width="0.8741258741258742"><path d="M 0.0,96.15384615384616 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,-17.482517482517483 l 17.482517482517483,17.482517482517483 l 17.482517482517483,17.482517482517483 " /></g></g></g></g></svg>+ "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"><svg stroke="rgb(0,0,0)" version="1.1" xmlns="http://www.w3.org/2000/svg" height="97.7011" viewBox="0 0 500 98" stroke-opacity="1" font-size="1" xmlns:xlink="http://www.w3.org/1999/xlink" width="500.0000"><defs></defs><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 453.7182,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 436.3027,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 418.8871,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 401.4716,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 384.0561,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 366.6405,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 349.2250,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 331.8095,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 314.3939,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 296.9784,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 279.5629,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 262.1473,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 244.7318,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 227.3163,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 209.9007,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 192.4852,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 175.0697,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 157.6541,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 140.2386,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 122.8231,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 105.4075,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 87.9920,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 70.5765,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 53.1609,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 35.7454,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 18.3299,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,87.9920 v -87.0777 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,0.9143 h 452.8039 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,18.3299 h 452.8039 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,35.7454 h 452.8039 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,53.1609 h 452.8039 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,70.5765 h 452.8039 "/></g><defs></defs><g stroke="rgb(128,128,128)" fill="rgb(0,0,0)" stroke-width="0.4353883664228493" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,87.9920 h 452.8039 "/></g><defs></defs><g stroke="rgb(255,0,0)" fill="rgb(0,0,0)" stroke-width="0.8707767328456986" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 0.9143,87.9920 l 17.4155,-17.4155 l 17.4155,-17.4155 l 17.4155,17.4155 l 17.4155,-17.4155 l 17.4155,-17.4155 l 17.4155,-17.4155 l 17.4155,17.4155 l 17.4155,-17.4155 l 17.4155,17.4155 l 17.4155,17.4155 l 17.4155,-17.4155 l 17.4155,17.4155 l 17.4155,-17.4155 l 17.4155,-17.4155 l 17.4155,-17.4155 l 17.4155,17.4155 l 17.4155,17.4155 l 17.4155,17.4155 l 17.4155,17.4155 l 17.4155,17.4155 l 17.4155,-17.4155 l 17.4155,17.4155 l 17.4155,-17.4155 l 17.4155,-17.4155 l 17.4155,17.4155 l 17.4155,17.4155 "/></g></svg>
+ svg/skew3.svg view
@@ -0,0 +1,3 @@+<?xml version="1.0" encoding="UTF-8"?>+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"+ "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"><svg stroke="rgb(0,0,0)" version="1.1" xmlns="http://www.w3.org/2000/svg" height="171.6044" viewBox="0 0 256 172" stroke-opacity="1" font-size="1" xmlns:xlink="http://www.w3.org/1999/xlink" width="256.0000"><defs></defs><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 231.5125,1.3427 v 25.5744 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 205.9381,1.3427 v 25.5744 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 180.3636,1.3427 v 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 154.7892,1.3427 v 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 129.2148,1.3427 v 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 103.6404,26.9171 v 25.5744 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 78.0659,26.9171 v 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 52.4915,52.4915 v 76.7233 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 26.9171,78.0659 v 76.7233 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.3427,103.6404 v 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.3427,154.7892 h 25.5744 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.3427,129.2148 h 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.3427,103.6404 h 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 26.9171,78.0659 h 51.1489 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 52.4915,52.4915 h 127.8721 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 78.0659,26.9171 h 153.4466 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.2787212787212787" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 129.2148,1.3427 h 102.2977 "/></g></svg>
+ svg/skew_tableau.svg view
@@ -0,0 +1,3 @@+<?xml version="1.0" encoding="UTF-8"?>+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"+ "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"><svg stroke="rgb(0,0,0)" version="1.1" xmlns="http://www.w3.org/2000/svg" height="214.5055" viewBox="0 0 320 215" stroke-opacity="1" font-size="1" xmlns:xlink="http://www.w3.org/1999/xlink" width="320.0000"><defs></defs><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 161.5185,1.6783 v 31.9680 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 129.5504,1.6783 v 31.9680 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 97.5824,1.6783 v 63.9361 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 65.6144,1.6783 v 95.9041 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 33.6464,1.6783 v 127.8721 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,1.6783 v 127.8721 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,129.5504 h 31.9680 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,97.5824 h 63.9361 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,65.6144 h 95.9041 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,33.6464 h 159.8402 "/></g><defs></defs><g stroke="rgb(211,211,211)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,1.6783 h 159.8402 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 289.3906,1.6783 v 31.9680 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 257.4226,1.6783 v 31.9680 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 225.4545,1.6783 v 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 193.4865,1.6783 v 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 161.5185,1.6783 v 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 129.5504,33.6464 v 31.9680 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 97.5824,33.6464 v 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 65.6144,65.6144 v 95.9041 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 33.6464,97.5824 v 95.9041 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,129.5504 v 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,193.4865 h 31.9680 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,161.5185 h 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 1.6783,129.5504 h 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 33.6464,97.5824 h 63.9361 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 65.6144,65.6144 h 159.8402 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 97.5824,33.6464 h 191.8082 "/></g><defs></defs><g stroke="rgb(0,0,0)" fill="rgb(0,0,0)" stroke-width="1.5984015984015982" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0"><path d="M 161.5185,1.6783 h 127.8721 "/></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,17.6623,186.4535)" dominant-baseline="middle">7</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,49.6304,154.4855)" dominant-baseline="middle">5</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,17.6623,154.4855)" dominant-baseline="middle">1</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,49.6304,122.5175)" dominant-baseline="middle">2</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,81.5984,90.5495)" dominant-baseline="middle">1</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,209.4705,58.5814)" dominant-baseline="middle">2</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,177.5025,58.5814)" dominant-baseline="middle">2</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,145.5345,58.5814)" dominant-baseline="middle">1</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,113.5664,58.5814)" dominant-baseline="middle">1</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,273.4066,26.6134)" dominant-baseline="middle">1</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,241.4386,26.6134)" dominant-baseline="middle">1</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,209.4705,26.6134)" dominant-baseline="middle">1</text></g><defs></defs><g stroke="rgb(0,0,255)" fill="rgb(0,0,255)" stroke-width="0.0" stroke-linecap="butt" stroke-miterlimit="10.0" stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="27.172827172827166px"><text stroke="none" text-anchor="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,177.5025,26.6134)" dominant-baseline="middle">1</text></g></svg>
svg/src/gen_figures.hs view
@@ -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 --------------------------------------------------------------------------------