packages feed

combinat 0.2.7.1 → 0.2.7.2

raw patch · 24 files changed

+2024/−1817 lines, 24 filesdep ~base

Dependency ranges changed: base

Files

Math/Combinat/ASCII.hs view
@@ -1,296 +1,300 @@-
--- | A mini-DSL for ASCII drawing of structures.
---
---
--- From some structures there is also Graphviz and\/or @diagrams@ 
--- (<http://projects.haskell.org/diagrams>) visualization support 
--- (the latter in the separate libray @combinat-diagrams@).
---
-
-module Math.Combinat.ASCII where
-
---------------------------------------------------------------------------------
-
-import Data.List
-
-import Math.Combinat.Helper
-
---------------------------------------------------------------------------------
--- * The basic type
-
--- | The type of a (rectangular) ASCII figure. 
--- Internally it is a list of lines of the same length plus the size.
---
--- Note: The Show instance is pretty-printing, so that it\'s convenient in ghci.
---
-data ASCII = ASCII 
-  { asciiSize  :: (Int,Int) 
-  , asciiLines :: [String]
-  }
-
-instance Show ASCII where
-  show = asciiString
-
--- | An empty (0x0) rectangle
-emptyRect :: ASCII
-emptyRect = ASCII (0,0) []
-
-asciiXSize, asciiYSize :: ASCII -> Int
-asciiXSize = fst . asciiSize
-asciiYSize = snd . asciiSize
-
-asciiString :: ASCII -> String
-asciiString (ASCII sz ls) = unlines ls
-
-printASCII :: ASCII -> IO ()
-printASCII = putStrLn . asciiString
-
-asciiFromLines :: [String] -> ASCII
-asciiFromLines ls = ASCII (x,y) (map f ls) where
-  y   = length ls
-  x   = maximum (map length ls)
-  f l = l ++ replicate (x - length l) ' '
-
-asciiFromString :: String -> ASCII
-asciiFromString = asciiFromLines . lines
-
---------------------------------------------------------------------------------
--- * Alignment
-
--- | Horizontal alignment
-data HAlign 
-  = HLeft 
-  | HCenter 
-  | HRight 
-  deriving (Eq,Show)
-
--- | Vertical alignment
-data VAlign 
-  = VTop 
-  | VCenter 
-  | VBottom 
-  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
-
--- | Horizontal separator
-data HSep 
-  = HSepEmpty           -- ^ empty separator
-  | HSepSpaces Int      -- ^ @n@ spaces
-  | HSepString String   -- ^ some custom string, eg. @\" | \"@
-  deriving Show
-
-hSepSize :: HSep -> Int
-hSepSize hsep = case hsep of
-  HSepEmpty    -> 0
-  HSepSpaces k -> k
-  HSepString s -> length s
-
-hSepString :: HSep -> String
-hSepString hsep = case hsep of
-  HSepEmpty    -> ""
-  HSepSpaces k -> replicate k ' '
-  HSepString s -> s
-
--- | Vertical separator
-data VSep 
-  = VSepEmpty           -- ^ empty separator
-  | VSepSpaces Int      -- ^ @n@ spaces
-  | VSepString [Char]   -- ^ some custom list of characters, eg. @\" - \"@ (the characters are interpreted as below each other)
-  deriving Show
-
-vSepSize :: VSep -> Int
-vSepSize vsep = case vsep of
-  VSepEmpty    -> 0
-  VSepSpaces k -> k
-  VSepString s -> length s
-
-vSepString :: VSep -> [Char]
-vSepString vsep = case vsep of
-  VSepEmpty    -> []
-  VSepSpaces k -> replicate k ' '
-  VSepString s -> s
-
---------------------------------------------------------------------------------
--- * 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 
-
---------------------------------------------------------------------------------
--- * Concatenation
-
--- | Horizontal concatenation
-hCatWith :: VAlign -> HSep -> [ASCII] -> ASCII
-hCatWith valign hsep rects = ASCII (x',maxy) final where
-  n    = length rects
-  maxy = maximum [ y | ASCII (_,y) _ <- rects ]
-  xsz  =         [ x | ASCII (x,_) _ <- rects ]
-  sep   = hSepString hsep
-  sepx  = length sep
-  rects1 = map (vExtendTo valign maxy) rects
-  x' = sum' xsz + (n-1)*sepx
-  final = map (intercalate sep) $ transpose (map asciiLines rects1)
-
--- | Vertical concatenation
-vCatWith :: HAlign -> VSep -> [ASCII] -> ASCII
-vCatWith halign vsep rects = ASCII (maxx,y') final where
-  n    = length rects
-  maxx = maximum [ x | ASCII (x,_) _ <- rects ]
-  ysz  =         [ y | ASCII (_,y) _ <- rects ]
-  sepy    = vSepSize vsep
-  fullsep = transpose (replicate maxx $ vSepString vsep) :: [String]
-  rects1  = map (hExtendTo halign maxx) rects
-  y'    = sum' ysz + (n-1)*sepy
-  final = intercalate fullsep $ map asciiLines rects1
-
---------------------------------------------------------------------------------
--- * Tabulate
-
-tabulate :: (HAlign,VAlign) -> (HSep,VSep) -> [[ASCII]] -> ASCII
-tabulate (halign,valign) (hsep,vsep) rects0 = final where
-  n = length rects0
-  m = maximum (map length rects0)
-  rects1 = map (\rs -> rs ++ replicate (m - length rs) emptyRect) rects0
-  ys = map (\rs -> maximum (map asciiYSize rs)) rects1
-  xs = map (\rs -> maximum (map asciiXSize rs)) (transpose rects1)
-  rects2 = map (\rs -> [      hExtendTo halign x  r  | (x,r ) <- zip xs rs     ]) rects1
-  rects3 =             [ map (vExtendTo valign y) rs | (y,rs) <- zip ys rects2 ]  
-  final  = vCatWith HLeft vsep 
-         $ map (hCatWith VTop hsep) rects3
-
--- | Order of elements in a matrix
-data MatrixOrder 
-  = RowMajor
-  | ColMajor
-  deriving (Eq,Ord,Show,Read)
-
--- | Automatically tabulates ASCII rectangles.
---
-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
-  -> [ASCII]          -- ^ list of ASCII rectangles
-  -> ASCII
-autoTabulate mtxorder ei list = final where
-  
-  final = tabulate (HLeft,VBottom) (HSepSpaces 2,VSepSpaces 1) rects 
-
-  n = length list
-
-  rects = case ei of
-
-    Left y  -> case mtxorder of
-                 ColMajor -> transpose (parts y list)
-                 RowMajor -> invparts y list
-
-    Right x -> case mtxorder of
-                 ColMajor -> transpose (invparts x list)
-                 RowMajor -> parts x list
-
-  transposeIf b = if b then transpose else id
-
-  -- chops into parts (the last one can be smaller)
-  parts d = go where
-    go [] = []
-    go xs = take d xs : go (drop d xs)
-
-  invparts d xs = parts' ds xs where
-    (q,r) = divMod n d
-    ds = replicate r (q+1) ++ replicate (d-r) q
-
-  parts' ds xs = go ds xs where
-    go _  [] = []                                      
-    go [] _  = []
-    go (d:ds) xs = take d xs : go ds (drop d xs)
-
---------------------------------------------------------------------------------
--- * Captions
-
--- | Adds a caption to the bottom, with default settings.
-caption :: String -> ASCII -> ASCII
-caption = caption' False HLeft
-
--- | Adds a caption to the bottom. The @Bool@ flag specifies whether to add an empty between 
--- the caption and the figure
-caption' :: Bool -> HAlign -> String -> ASCII -> ASCII
-caption' emptyline halign str rect = vCatWith halign sep [rect,capt] where
-  sep  = if emptyline then VSepSpaces 1 else VSepEmpty 
-  capt = asciiFromString str
-
---------------------------------------------------------------------------------
--- * Testing \/ miscellanea
-
--- | An ASCII 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
-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) '-' ++ "/"
-
-asciiNumber :: Int -> ASCII
-asciiNumber = asciiShow
-
-asciiShow :: Show a => a -> ASCII
-asciiShow = asciiFromLines . (:[]) . show
-
---------------------------------------------------------------------------------
++-- | A mini-DSL for ASCII drawing of structures.+--+--+-- From some structures there is also Graphviz and\/or @diagrams@ +-- (<http://projects.haskell.org/diagrams>) visualization support +-- (the latter in the separate libray @combinat-diagrams@).+--++module Math.Combinat.ASCII where++--------------------------------------------------------------------------------++import Data.List++import Math.Combinat.Helper++--------------------------------------------------------------------------------+-- * The basic type++-- | The type of a (rectangular) ASCII figure. +-- Internally it is a list of lines of the same length plus the size.+--+-- Note: The Show instance is pretty-printing, so that it\'s convenient in ghci.+--+data ASCII = ASCII +  { asciiSize  :: (Int,Int) +  , asciiLines :: [String]+  }++-- | A type class to have a simple way to draw things +class DrawASCII a where+  ascii :: a -> ASCII++instance Show ASCII where+  show = asciiString++-- | An empty (0x0) rectangle+emptyRect :: ASCII+emptyRect = ASCII (0,0) []++asciiXSize, asciiYSize :: ASCII -> Int+asciiXSize = fst . asciiSize+asciiYSize = snd . asciiSize++asciiString :: ASCII -> String+asciiString (ASCII sz ls) = unlines ls++printASCII :: ASCII -> IO ()+printASCII = putStrLn . asciiString++asciiFromLines :: [String] -> ASCII+asciiFromLines ls = ASCII (x,y) (map f ls) where+  y   = length ls+  x   = maximum (map length ls)+  f l = l ++ replicate (x - length l) ' '++asciiFromString :: String -> ASCII+asciiFromString = asciiFromLines . lines++--------------------------------------------------------------------------------+-- * Alignment++-- | Horizontal alignment+data HAlign +  = HLeft +  | HCenter +  | HRight +  deriving (Eq,Show)++-- | Vertical alignment+data VAlign +  = VTop +  | VCenter +  | VBottom +  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++-- | Horizontal separator+data HSep +  = HSepEmpty           -- ^ empty separator+  | HSepSpaces Int      -- ^ @n@ spaces+  | HSepString String   -- ^ some custom string, eg. @\" | \"@+  deriving Show++hSepSize :: HSep -> Int+hSepSize hsep = case hsep of+  HSepEmpty    -> 0+  HSepSpaces k -> k+  HSepString s -> length s++hSepString :: HSep -> String+hSepString hsep = case hsep of+  HSepEmpty    -> ""+  HSepSpaces k -> replicate k ' '+  HSepString s -> s++-- | Vertical separator+data VSep +  = VSepEmpty           -- ^ empty separator+  | VSepSpaces Int      -- ^ @n@ spaces+  | VSepString [Char]   -- ^ some custom list of characters, eg. @\" - \"@ (the characters are interpreted as below each other)+  deriving Show++vSepSize :: VSep -> Int+vSepSize vsep = case vsep of+  VSepEmpty    -> 0+  VSepSpaces k -> k+  VSepString s -> length s++vSepString :: VSep -> [Char]+vSepString vsep = case vsep of+  VSepEmpty    -> []+  VSepSpaces k -> replicate k ' '+  VSepString s -> s++--------------------------------------------------------------------------------+-- * 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 ++--------------------------------------------------------------------------------+-- * Concatenation++-- | Horizontal concatenation+hCatWith :: VAlign -> HSep -> [ASCII] -> ASCII+hCatWith valign hsep rects = ASCII (x',maxy) final where+  n    = length rects+  maxy = maximum [ y | ASCII (_,y) _ <- rects ]+  xsz  =         [ x | ASCII (x,_) _ <- rects ]+  sep   = hSepString hsep+  sepx  = length sep+  rects1 = map (vExtendTo valign maxy) rects+  x' = sum' xsz + (n-1)*sepx+  final = map (intercalate sep) $ transpose (map asciiLines rects1)++-- | Vertical concatenation+vCatWith :: HAlign -> VSep -> [ASCII] -> ASCII+vCatWith halign vsep rects = ASCII (maxx,y') final where+  n    = length rects+  maxx = maximum [ x | ASCII (x,_) _ <- rects ]+  ysz  =         [ y | ASCII (_,y) _ <- rects ]+  sepy    = vSepSize vsep+  fullsep = transpose (replicate maxx $ vSepString vsep) :: [String]+  rects1  = map (hExtendTo halign maxx) rects+  y'    = sum' ysz + (n-1)*sepy+  final = intercalate fullsep $ map asciiLines rects1++--------------------------------------------------------------------------------+-- * Tabulate++tabulate :: (HAlign,VAlign) -> (HSep,VSep) -> [[ASCII]] -> ASCII+tabulate (halign,valign) (hsep,vsep) rects0 = final where+  n = length rects0+  m = maximum (map length rects0)+  rects1 = map (\rs -> rs ++ replicate (m - length rs) emptyRect) rects0+  ys = map (\rs -> maximum (map asciiYSize rs)) rects1+  xs = map (\rs -> maximum (map asciiXSize rs)) (transpose rects1)+  rects2 = map (\rs -> [      hExtendTo halign x  r  | (x,r ) <- zip xs rs     ]) rects1+  rects3 =             [ map (vExtendTo valign y) rs | (y,rs) <- zip ys rects2 ]  +  final  = vCatWith HLeft vsep +         $ map (hCatWith VTop hsep) rects3++-- | Order of elements in a matrix+data MatrixOrder +  = RowMajor+  | ColMajor+  deriving (Eq,Ord,Show,Read)++-- | Automatically tabulates ASCII rectangles.+--+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+  -> [ASCII]          -- ^ list of ASCII rectangles+  -> ASCII+autoTabulate mtxorder ei list = final where+  +  final = tabulate (HLeft,VBottom) (HSepSpaces 2,VSepSpaces 1) rects ++  n = length list++  rects = case ei of++    Left y  -> case mtxorder of+                 ColMajor -> transpose (parts y list)+                 RowMajor -> invparts y list++    Right x -> case mtxorder of+                 ColMajor -> transpose (invparts x list)+                 RowMajor -> parts x list++  transposeIf b = if b then transpose else id++  -- chops into parts (the last one can be smaller)+  parts d = go where+    go [] = []+    go xs = take d xs : go (drop d xs)++  invparts d xs = parts' ds xs where+    (q,r) = divMod n d+    ds = replicate r (q+1) ++ replicate (d-r) q++  parts' ds xs = go ds xs where+    go _  [] = []                                      +    go [] _  = []+    go (d:ds) xs = take d xs : go ds (drop d xs)++--------------------------------------------------------------------------------+-- * Captions++-- | Adds a caption to the bottom, with default settings.+caption :: String -> ASCII -> ASCII+caption = caption' False HLeft++-- | Adds a caption to the bottom. The @Bool@ flag specifies whether to add an empty between +-- the caption and the figure+caption' :: Bool -> HAlign -> String -> ASCII -> ASCII+caption' emptyline halign str rect = vCatWith halign sep [rect,capt] where+  sep  = if emptyline then VSepSpaces 1 else VSepEmpty +  capt = asciiFromString str++--------------------------------------------------------------------------------+-- * Testing \/ miscellanea++-- | An ASCII 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+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) '-' ++ "/"++asciiNumber :: Int -> ASCII+asciiNumber = asciiShow++asciiShow :: Show a => a -> ASCII+asciiShow = asciiFromLines . (:[]) . show++--------------------------------------------------------------------------------
Math/Combinat/Compositions.hs view
@@ -70,13 +70,13 @@   => a       -- ^ length   -> a       -- ^ sum   -> [[Int]]-compositions1 len' d' -  | len > d = []+compositions1 len d +  | len > d   = []   | otherwise = map plus1 $ compositions len (d-len)   where     plus1 = map (+1)-    len = fromIntegral len'-    d   = fromIntegral d'+    -- len = fromIntegral len'+    -- d   = fromIntegral d'  countCompositions1 :: Integral a => a -> a -> Integer countCompositions1 len d = countCompositions len (d-len)
Math/Combinat/FreeGroups.hs view
@@ -2,10 +2,19 @@ -- | Words in free groups (and free powers of cyclic groups). -- This module is not re-exported by "Math.Combinat" ---{-# LANGUAGE PatternGuards #-}+{-# 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 )
Math/Combinat/Helper.hs view
@@ -24,8 +24,6 @@ -------------------------------------------------------------------------------- -- * pairs -{-# SPECIALIZE swap :: (a  ,a  ) -> (a  ,a  ) #-}-{-# SPECIALIZE swap :: (Int,Int) -> (Int,Int) #-} swap :: (a,b) -> (b,a) swap (x,y) = (y,x)
Math/Combinat/LatticePaths.hs view
@@ -1,368 +1,379 @@-
--- | Dyck paths, lattice paths, etc
---
--- For example, the following figure represents a Dyck path of height 5 with 3 zero-touches (not counting the starting point,
--- but counting the endpoint) and 7 peaks:
---
--- <<svg/dyck_path.svg>>
---
-
-{-# LANGUAGE BangPatterns #-}
-module Math.Combinat.LatticePaths where
-
---------------------------------------------------------------------------------
-
-import Data.List
-import System.Random
-
-import Math.Combinat.Numbers
-import Math.Combinat.Trees.Binary
-import Math.Combinat.ASCII as ASCII
-
---------------------------------------------------------------------------------
--- * Types
-
--- | A step in a lattice path
-data Step 
-  = UpStep         -- ^ the step @(1,1)@
-  | DownStep       -- ^ the step @(1,-1)@
-  deriving (Eq,Ord,Show)
-
--- | A lattice path is a path using only the allowed steps, never going below the zero level line @y=0@. 
---
--- Note that if you rotate such a path by 45 degrees counterclockwise,
--- you get a path which uses only the steps @(1,0)@ and @(0,1)@, and stays
--- above the main diagonal (hence the name, we just use a different convention).
---
-type LatticePath = [Step]
-
---------------------------------------------------------------------------------
--- * ascii drawing of paths
-
--- | Draws the path into a list of lines. For example try:
---
--- > autotabulate RowMajor (Right 5) (map asciiPath $ dyckPaths 4)
---
-asciiPath :: LatticePath -> ASCII
-asciiPath p = asciiFromLines $ transpose (go 0 p) where
-
-  go !h [] = []
-  go !h (x:xs) = case x of
-    UpStep   -> ee  h    x : go (h+1) xs
-    DownStep -> ee (h-1) x : go (h-1) xs
-
-  maxh   = pathHeight p
-
-  ee h x = replicate (maxh-h-1) ' ' ++ [ch x] ++ replicate h ' '
-  ch x   = case x of 
-    UpStep   -> '/' 
-    DownStep -> '\\' 
-
---------------------------------------------------------------------------------
--- * elementary queries
-
--- | A lattice path is called \"valid\", if it never goes below the @y=0@ line.
-isValidPath :: LatticePath -> Bool
-isValidPath = go 0 where
-  go !y []     = y>=0
-  go !y (t:ts) = let y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }
-                 in  if y'<0 then False 
-                             else go y' ts
-
--- | A Dyck path is a lattice path whose last point lies on the @y=0@ line
-isDyckPath :: LatticePath -> Bool
-isDyckPath = go 0 where
-  go !y []     = y==0
-  go !y (t:ts) = let y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }
-                 in  if y'<0 then False 
-                             else go y' ts
-
--- | Maximal height of a lattice path
-pathHeight :: LatticePath -> Int
-pathHeight = go 0 0 where
-  go !h !y []     = h
-  go !h !y (t:ts) = case t of
-    UpStep   -> go (max h (y+1)) (y+1) ts
-    DownStep -> go      h        (y-1) ts
-
--- | Endpoint of a lattice path, which starts from @(0,0)@.
-pathEndpoint :: LatticePath -> (Int,Int)
-pathEndpoint = go 0 0 where
-  go !x !y []     = (x,y)
-  go !x !y (t:ts) = case t of                         
-    UpStep   -> go (x+1) (y+1) ts
-    DownStep -> go (x+1) (y-1) ts
-
--- | Returns the coordinates of the path (excluding the starting point @(0,0)@, but including
--- the endpoint)
-pathCoordinates :: LatticePath -> [(Int,Int)]
-pathCoordinates = go 0 0 where
-  go _  _  []     = []
-  go !x !y (t:ts) = let x' = x + 1
-                        y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }
-                    in  (x',y') : go x' y' ts
-
--- | Counts the up-steps
-pathNumberOfUpSteps :: LatticePath -> Int
-pathNumberOfUpSteps   = fst . pathNumberOfUpDownSteps
-
--- | Counts the down-steps
-pathNumberOfDownSteps :: LatticePath -> Int
-pathNumberOfDownSteps = snd . pathNumberOfUpDownSteps
-
--- | Counts both the up-steps and down-steps
-pathNumberOfUpDownSteps :: LatticePath -> (Int,Int)
-pathNumberOfUpDownSteps = go 0 0 where 
-  go !u !d (p:ps) = case p of 
-    UpStep   -> go (u+1)  d    ps  
-    DownStep -> go  u    (d+1) ps    
-  go !u !d []     = (u,d)
-
---------------------------------------------------------------------------------
--- * path-specific queries
-
--- | Number of peaks of a path (excluding the endpoint)
-pathNumberOfPeaks :: LatticePath -> Int
-pathNumberOfPeaks = go 0 where
-  go !k (x:xs@(y:_)) = go (if x==UpStep && y==DownStep then k+1 else k) xs
-  go !k [x] = k
-  go !k [ ] = k
-
--- | Number of points on the path which touch the @y=0@ zero level line
--- (excluding the starting point @(0,0)@, but including the endpoint; that is, for Dyck paths it this is always positive!).
-pathNumberOfZeroTouches :: LatticePath -> Int
-pathNumberOfZeroTouches = pathNumberOfTouches' 0
-
--- | Number of points on the path which touch the level line at height @h@
--- (excluding the starting point @(0,0)@, but including the endpoint).
-pathNumberOfTouches' 
-  :: Int       -- ^ @h@ = the touch level
-  -> LatticePath -> Int
-pathNumberOfTouches' h = go 0 0 0 where
-  go !cnt _  _  []     = cnt
-  go !cnt !x !y (t:ts) = let y'   = case t of { UpStep -> y+1 ; DownStep -> y-1 }
-                             cnt' = if y'==h then cnt+1 else cnt
-                         in  go cnt' (x+1) y' ts
-
---------------------------------------------------------------------------------
--- * Dyck paths
-
--- | @dyckPaths m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@. 
--- 
--- Remark: Dyck paths are obviously in bijection with nested parentheses, and thus
--- also with binary trees.
---
--- Order is reverse lexicographical:
---
--- > sort (dyckPaths m) == reverse (dyckPaths m)
--- 
-dyckPaths :: Int -> [LatticePath]
-dyckPaths = map nestedParensToDyckPath . nestedParentheses 
-
--- | @dyckPaths m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@. 
---
--- > sort (dyckPathsNaive m) == sort (dyckPaths m) 
---  
--- Naive recursive algorithm, order is ad-hoc
---
-dyckPathsNaive :: Int -> [LatticePath]
-dyckPathsNaive = worker where
-  worker  0 = [[]]
-  worker  m = as ++ bs where
-    as = [ bracket p      | p <- worker (m-1) ] 
-    bs = [ bracket p ++ q | k <- [1..m-1] , p <- worker (k-1) , q <- worker (m-k) ]
-  bracket p = UpStep : p ++ [DownStep]
-
--- | The number of Dyck paths from @(0,0)@ to @(2m,0)@ is simply the m\'th Catalan number.
-countDyckPaths :: Int -> Integer
-countDyckPaths m = catalan m
-
--- | The trivial bijection
-nestedParensToDyckPath :: [Paren] -> LatticePath
-nestedParensToDyckPath = map f where
-  f p = case p of { LeftParen -> UpStep ; RightParen -> DownStep }
-
--- | The trivial bijection in the other direction
-dyckPathToNestedParens :: LatticePath -> [Paren]
-dyckPathToNestedParens = map g where
-  g s = case s of { UpStep -> LeftParen ; DownStep -> RightParen }
-
---------------------------------------------------------------------------------
--- * Bounded Dyck paths
-
--- | @boundedDyckPaths h m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ whose height is at most @h@.
--- Synonym for 'boundedDyckPathsNaive'.
---
-boundedDyckPaths
-  :: Int   -- ^ @h@ = maximum height
-  -> Int   -- ^ @m@ = half-length
-  -> [LatticePath]
-boundedDyckPaths = boundedDyckPathsNaive 
-
--- | @boundedDyckPathsNaive h m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ whose height is at most @h@.
---
--- > sort (boundedDyckPaths h m) == sort [ p | p <- dyckPaths m , pathHeight p <= h ]
--- > sort (boundedDyckPaths m m) == sort (dyckPaths m) 
---
--- Naive recursive algorithm, resulting order is pretty ad-hoc.
---
-boundedDyckPathsNaive
-  :: Int   -- ^ @h@ = maximum height
-  -> Int   -- ^ @m@ = half-length
-  -> [LatticePath]
-boundedDyckPathsNaive = worker where
-  worker !h !m 
-    | h<0        = []
-    | m<0        = []
-    | m==0       = [[]]
-    | h<=0       = []
-    | otherwise  = as ++ bs 
-    where
-      bracket p = UpStep : p ++ [DownStep]
-      as = [ bracket p      |                 p <- boundedDyckPaths (h-1) (m-1)                                 ]
-      bs = [ bracket p ++ q | k <- [1..m-1] , p <- boundedDyckPaths (h-1) (k-1) , q <- boundedDyckPaths h (m-k) ]
-
---------------------------------------------------------------------------------
--- * More general lattice paths
-
--- | All lattice paths from @(0,0)@ to @(x,y)@. Clearly empty unless @x-y@ is even.
--- Synonym for 'latticePathsNaive'
---
-latticePaths :: (Int,Int) -> [LatticePath]
-latticePaths = latticePathsNaive
-
--- | All lattice paths from @(0,0)@ to @(x,y)@. Clearly empty unless @x-y@ is even.
---
--- Note that
---
--- > sort (dyckPaths n) == sort (latticePaths (0,2*n))
---
--- Naive recursive algorithm, resulting order is pretty ad-hoc.
---
-latticePathsNaive :: (Int,Int) -> [LatticePath]
-latticePathsNaive (x,y) = worker x y where
-  worker !x !y 
-    | odd (x-y)     = []
-    | x<0           = []
-    | y<0           = []
-    | y==0          = dyckPaths (div x 2)
-    | x==1 && y==1  = [[UpStep]]
-    | otherwise     = as ++ bs
-    where
-      bracket p = UpStep : p ++ [DownStep] 
-      as = [ UpStep : p     | p <- worker (x-1) (y-1) ]
-      bs = [ bracket p ++ q | k <- [1..(div x 2)] , p <- dyckPaths (k-1) , q <- worker (x-2*k) y ]
-
--- | Lattice paths are counted by the numbers in the Catalan triangle.
-countLatticePaths :: (Int,Int) -> Integer
-countLatticePaths (x,y) 
-  | even (x+y)  = catalanTriangle (div (x+y) 2) (div (x-y) 2)
-  | otherwise   = 0
-
---------------------------------------------------------------------------------
--- * Zero-level touches
-
--- | @touchingDyckPaths k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ which touch the 
--- zero level line @y=0@ exactly @k@ times (excluding the starting point, but including the endpoint;
--- thus, @k@ should be positive). Synonym for 'touchingDyckPathsNaive'.
-touchingDyckPaths
-  :: Int   -- ^ @k@ = number of zero-touches
-  -> Int   -- ^ @m@ = half-length
-  -> [LatticePath]
-touchingDyckPaths = touchingDyckPathsNaive
-
-
--- | @touchingDyckPathsNaive k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ which touch the 
--- zero level line @y=0@ exactly @k@ times (excluding the starting point, but including the endpoint;
--- thus, @k@ should be positive).
---
--- > sort (touchingDyckPathsNaive k m) == sort [ p | p <- dyckPaths m , pathNumberOfZeroTouches p == k ]
--- 
--- Naive recursive algorithm, resulting order is pretty ad-hoc.
---
-touchingDyckPathsNaive
-  :: Int   -- ^ @k@ = number of zero-touches
-  -> Int   -- ^ @m@ = half-length
-  -> [LatticePath]
-touchingDyckPathsNaive = worker where
-  worker !k !m 
-    | m == 0    = if k==0 then [[]] else []
-    | k <= 0    = []
-    | m <  0    = []
-    | k == 1    = [ bracket p      |                 p <- dyckPaths (m-1)                           ]
-    | otherwise = [ bracket p ++ q | l <- [1..m-1] , p <- dyckPaths (l-1) , q <- worker (k-1) (m-l) ]
-    where
-      bracket p = UpStep : p ++ [DownStep] 
-
-
--- | There is a bijection from the set of non-empty Dyck paths of length @2n@ which touch the zero lines @t@ times,
--- to lattice paths from @(0,0)@ to @(2n-t-1,t-1)@ (just remove all the down-steps just before touching
--- the zero line, and also the very first up-step). This gives us a counting formula.
-countTouchingDyckPaths 
-  :: Int   -- ^ @k@ = number of zero-touches
-  -> Int   -- ^ @m@ = half-length
-  -> Integer
-countTouchingDyckPaths t n
-  | t==0 && n==0   = 1
-  | otherwise      = countLatticePaths (2*n-t-1,t-1)
-
---------------------------------------------------------------------------------
--- * Dyck paths with given number of peaks
-
--- | @peakingDyckPaths k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ with exactly @k@ peaks.
---
--- Synonym for 'peakingDyckPathsNaive'
---
-peakingDyckPaths
-  :: Int      -- ^ @k@ = number of peaks
-  -> Int      -- ^ @m@ = half-length
-  -> [LatticePath]
-peakingDyckPaths = peakingDyckPathsNaive 
-
--- | @peakingDyckPathsNaive k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ with exactly @k@ peaks.
---
--- > sort (peakingDyckPathsNaive k m) = sort [ p | p <- dyckPaths m , pathNumberOfPeaks p == k ]
---  
--- Naive recursive algorithm, resulting order is pretty ad-hoc.
---
-peakingDyckPathsNaive 
-  :: Int      -- ^ @k@ = number of peaks
-  -> Int      -- ^ @m@ = half-length
-  -> [LatticePath]
-peakingDyckPathsNaive = worker where
-  worker !k !m
-    | m == 0    = if k==0 then [[]] else []       
-    | k <= 0    = []
-    | m <  0    = []
-    | k == 1    = [ singlePeak m ] 
-    | otherwise = as ++ bs ++ cs
-    where
-      as = [ bracket p      |                                 p <- worker k (m-1)                           ]
-      bs = [ smallHill ++ q |                                                       q <- worker (k-1) (m-1) ]
-      cs = [ bracket p ++ q | l <- [2..m-1] , a <- [1..k-1] , p <- worker a (l-1) , q <- worker (k-a) (m-l) ]
-      smallHill     = [ UpStep , DownStep ]
-      singlePeak !m = replicate m UpStep ++ replicate m DownStep 
-      bracket p = UpStep : p ++ [DownStep] 
-
--- | Dyck paths of length @2m@ with @k@ peaks are counted by the Narayana numbers @N(m,k) = \binom{m}{k} \binom{m}{k-1} / m@
-countPeakingDyckPaths
-  :: Int      -- ^ @k@ = number of peaks
-  -> Int      -- ^ @m@ = half-length
-  -> Integer
-countPeakingDyckPaths k m 
-  | m == 0    = if k==0 then 1 else 0
-  | k <= 0    = 0
-  | m <  0    = 0
-  | k == 1    = 1
-  | otherwise = div (binomial m k * binomial m (k-1)) (fromIntegral m)
-
---------------------------------------------------------------------------------
--- * Random lattice paths
-
--- | A uniformly random Dyck path of length @2m@
-randomDyckPath :: RandomGen g => Int -> g -> (LatticePath,g)
-randomDyckPath m g0 = (nestedParensToDyckPath parens, g1) where
-  (parens,g1) = randomNestedParentheses m g0
-
---------------------------------------------------------------------------------
-
++-- | Dyck paths, lattice paths, etc+--+-- For example, the following figure represents a Dyck path of height 5 with 3 zero-touches (not counting the starting point,+-- but counting the endpoint) and 7 peaks:+--+-- <<svg/dyck_path.svg>>+--++{-# LANGUAGE BangPatterns, FlexibleInstances, TypeSynonymInstances #-}+module Math.Combinat.LatticePaths where++--------------------------------------------------------------------------------++import Data.List+import System.Random++import Math.Combinat.Numbers+import Math.Combinat.Trees.Binary+import Math.Combinat.ASCII as ASCII++--------------------------------------------------------------------------------+-- * Types++-- | A step in a lattice path+data Step +  = UpStep         -- ^ the step @(1,1)@+  | DownStep       -- ^ the step @(1,-1)@+  deriving (Eq,Ord,Show)++-- | A lattice path is a path using only the allowed steps, never going below the zero level line @y=0@. +--+-- Note that if you rotate such a path by 45 degrees counterclockwise,+-- you get a path which uses only the steps @(1,0)@ and @(0,1)@, and stays+-- above the main diagonal (hence the name, we just use a different convention).+--+type LatticePath = [Step]++--------------------------------------------------------------------------------+-- * ascii drawing of paths++-- | Draws the path into a list of lines. For example try:+--+-- > autotabulate RowMajor (Right 5) (map asciiPath $ dyckPaths 4)+--+asciiPath :: LatticePath -> ASCII+asciiPath p = asciiFromLines $ transpose (go 0 p) where++  go !h [] = []+  go !h (x:xs) = case x of+    UpStep   -> ee  h    x : go (h+1) xs+    DownStep -> ee (h-1) x : go (h-1) xs++  maxh   = pathHeight p++  ee h x = replicate (maxh-h-1) ' ' ++ [ch x] ++ replicate h ' '+  ch x   = case x of +    UpStep   -> '/' +    DownStep -> '\\' ++instance DrawASCII LatticePath where +  ascii = asciiPath++--------------------------------------------------------------------------------+-- * elementary queries++-- | A lattice path is called \"valid\", if it never goes below the @y=0@ line.+isValidPath :: LatticePath -> Bool+isValidPath = go 0 where+  go :: Int -> LatticePath -> Bool+  go !y []     = y>=0+  go !y (t:ts) = let y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }+                 in  if y'<0 then False +                             else go y' ts++-- | A Dyck path is a lattice path whose last point lies on the @y=0@ line+isDyckPath :: LatticePath -> Bool+isDyckPath = go 0 where+  go :: Int -> LatticePath -> Bool+  go !y []     = y==0+  go !y (t:ts) = let y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }+                 in  if y'<0 then False +                             else go y' ts++-- | Maximal height of a lattice path+pathHeight :: LatticePath -> Int+pathHeight = go 0 0 where+  go :: Int -> Int -> LatticePath -> Int+  go !h !y []     = h+  go !h !y (t:ts) = case t of+    UpStep   -> go (max h (y+1)) (y+1) ts+    DownStep -> go      h        (y-1) ts++-- | Endpoint of a lattice path, which starts from @(0,0)@.+pathEndpoint :: LatticePath -> (Int,Int)+pathEndpoint = go 0 0 where+  go :: Int -> Int -> LatticePath -> (Int,Int)+  go !x !y []     = (x,y)+  go !x !y (t:ts) = case t of                         +    UpStep   -> go (x+1) (y+1) ts+    DownStep -> go (x+1) (y-1) ts++-- | Returns the coordinates of the path (excluding the starting point @(0,0)@, but including+-- the endpoint)+pathCoordinates :: LatticePath -> [(Int,Int)]+pathCoordinates = go 0 0 where+  go :: Int -> Int -> LatticePath -> [(Int,Int)]+  go _  _  []     = []+  go !x !y (t:ts) = let x' = x + 1+                        y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }+                    in  (x',y') : go x' y' ts++-- | Counts the up-steps+pathNumberOfUpSteps :: LatticePath -> Int+pathNumberOfUpSteps   = fst . pathNumberOfUpDownSteps++-- | Counts the down-steps+pathNumberOfDownSteps :: LatticePath -> Int+pathNumberOfDownSteps = snd . pathNumberOfUpDownSteps++-- | Counts both the up-steps and down-steps+pathNumberOfUpDownSteps :: LatticePath -> (Int,Int)+pathNumberOfUpDownSteps = go 0 0 where +  go :: Int -> Int -> LatticePath -> (Int,Int)+  go !u !d (p:ps) = case p of +    UpStep   -> go (u+1)  d    ps  +    DownStep -> go  u    (d+1) ps    +  go !u !d []     = (u,d)++--------------------------------------------------------------------------------+-- * path-specific queries++-- | Number of peaks of a path (excluding the endpoint)+pathNumberOfPeaks :: LatticePath -> Int+pathNumberOfPeaks = go 0 where+  go :: Int -> LatticePath -> Int+  go !k (x:xs@(y:_)) = go (if x==UpStep && y==DownStep then k+1 else k) xs+  go !k [x] = k+  go !k [ ] = k++-- | Number of points on the path which touch the @y=0@ zero level line+-- (excluding the starting point @(0,0)@, but including the endpoint; that is, for Dyck paths it this is always positive!).+pathNumberOfZeroTouches :: LatticePath -> Int+pathNumberOfZeroTouches = pathNumberOfTouches' 0++-- | Number of points on the path which touch the level line at height @h@+-- (excluding the starting point @(0,0)@, but including the endpoint).+pathNumberOfTouches' +  :: Int       -- ^ @h@ = the touch level+  -> LatticePath -> Int+pathNumberOfTouches' h = go 0 0 0 where+  go :: Int -> Int -> Int -> LatticePath -> Int+  go !cnt _  _  []     = cnt+  go !cnt !x !y (t:ts) = let y'   = case t of { UpStep -> y+1 ; DownStep -> y-1 }+                             cnt' = if y'==h then cnt+1 else cnt+                         in  go cnt' (x+1) y' ts++--------------------------------------------------------------------------------+-- * Dyck paths++-- | @dyckPaths m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@. +-- +-- Remark: Dyck paths are obviously in bijection with nested parentheses, and thus+-- also with binary trees.+--+-- Order is reverse lexicographical:+--+-- > sort (dyckPaths m) == reverse (dyckPaths m)+-- +dyckPaths :: Int -> [LatticePath]+dyckPaths = map nestedParensToDyckPath . nestedParentheses ++-- | @dyckPaths m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@. +--+-- > sort (dyckPathsNaive m) == sort (dyckPaths m) +--  +-- Naive recursive algorithm, order is ad-hoc+--+dyckPathsNaive :: Int -> [LatticePath]+dyckPathsNaive = worker where+  worker  0 = [[]]+  worker  m = as ++ bs where+    as = [ bracket p      | p <- worker (m-1) ] +    bs = [ bracket p ++ q | k <- [1..m-1] , p <- worker (k-1) , q <- worker (m-k) ]+  bracket p = UpStep : p ++ [DownStep]++-- | The number of Dyck paths from @(0,0)@ to @(2m,0)@ is simply the m\'th Catalan number.+countDyckPaths :: Int -> Integer+countDyckPaths m = catalan m++-- | The trivial bijection+nestedParensToDyckPath :: [Paren] -> LatticePath+nestedParensToDyckPath = map f where+  f p = case p of { LeftParen -> UpStep ; RightParen -> DownStep }++-- | The trivial bijection in the other direction+dyckPathToNestedParens :: LatticePath -> [Paren]+dyckPathToNestedParens = map g where+  g s = case s of { UpStep -> LeftParen ; DownStep -> RightParen }++--------------------------------------------------------------------------------+-- * Bounded Dyck paths++-- | @boundedDyckPaths h m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ whose height is at most @h@.+-- Synonym for 'boundedDyckPathsNaive'.+--+boundedDyckPaths+  :: Int   -- ^ @h@ = maximum height+  -> Int   -- ^ @m@ = half-length+  -> [LatticePath]+boundedDyckPaths = boundedDyckPathsNaive ++-- | @boundedDyckPathsNaive h m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ whose height is at most @h@.+--+-- > sort (boundedDyckPaths h m) == sort [ p | p <- dyckPaths m , pathHeight p <= h ]+-- > sort (boundedDyckPaths m m) == sort (dyckPaths m) +--+-- Naive recursive algorithm, resulting order is pretty ad-hoc.+--+boundedDyckPathsNaive+  :: Int   -- ^ @h@ = maximum height+  -> Int   -- ^ @m@ = half-length+  -> [LatticePath]+boundedDyckPathsNaive = worker where+  worker !h !m +    | h<0        = []+    | m<0        = []+    | m==0       = [[]]+    | h<=0       = []+    | otherwise  = as ++ bs +    where+      bracket p = UpStep : p ++ [DownStep]+      as = [ bracket p      |                 p <- boundedDyckPaths (h-1) (m-1)                                 ]+      bs = [ bracket p ++ q | k <- [1..m-1] , p <- boundedDyckPaths (h-1) (k-1) , q <- boundedDyckPaths h (m-k) ]++--------------------------------------------------------------------------------+-- * More general lattice paths++-- | All lattice paths from @(0,0)@ to @(x,y)@. Clearly empty unless @x-y@ is even.+-- Synonym for 'latticePathsNaive'+--+latticePaths :: (Int,Int) -> [LatticePath]+latticePaths = latticePathsNaive++-- | All lattice paths from @(0,0)@ to @(x,y)@. Clearly empty unless @x-y@ is even.+--+-- Note that+--+-- > sort (dyckPaths n) == sort (latticePaths (0,2*n))+--+-- Naive recursive algorithm, resulting order is pretty ad-hoc.+--+latticePathsNaive :: (Int,Int) -> [LatticePath]+latticePathsNaive (x,y) = worker x y where+  worker !x !y +    | odd (x-y)     = []+    | x<0           = []+    | y<0           = []+    | y==0          = dyckPaths (div x 2)+    | x==1 && y==1  = [[UpStep]]+    | otherwise     = as ++ bs+    where+      bracket p = UpStep : p ++ [DownStep] +      as = [ UpStep : p     | p <- worker (x-1) (y-1) ]+      bs = [ bracket p ++ q | k <- [1..(div x 2)] , p <- dyckPaths (k-1) , q <- worker (x-2*k) y ]++-- | Lattice paths are counted by the numbers in the Catalan triangle.+countLatticePaths :: (Int,Int) -> Integer+countLatticePaths (x,y) +  | even (x+y)  = catalanTriangle (div (x+y) 2) (div (x-y) 2)+  | otherwise   = 0++--------------------------------------------------------------------------------+-- * Zero-level touches++-- | @touchingDyckPaths k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ which touch the +-- zero level line @y=0@ exactly @k@ times (excluding the starting point, but including the endpoint;+-- thus, @k@ should be positive). Synonym for 'touchingDyckPathsNaive'.+touchingDyckPaths+  :: Int   -- ^ @k@ = number of zero-touches+  -> Int   -- ^ @m@ = half-length+  -> [LatticePath]+touchingDyckPaths = touchingDyckPathsNaive+++-- | @touchingDyckPathsNaive k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ which touch the +-- zero level line @y=0@ exactly @k@ times (excluding the starting point, but including the endpoint;+-- thus, @k@ should be positive).+--+-- > sort (touchingDyckPathsNaive k m) == sort [ p | p <- dyckPaths m , pathNumberOfZeroTouches p == k ]+-- +-- Naive recursive algorithm, resulting order is pretty ad-hoc.+--+touchingDyckPathsNaive+  :: Int   -- ^ @k@ = number of zero-touches+  -> Int   -- ^ @m@ = half-length+  -> [LatticePath]+touchingDyckPathsNaive = worker where+  worker !k !m +    | m == 0    = if k==0 then [[]] else []+    | k <= 0    = []+    | m <  0    = []+    | k == 1    = [ bracket p      |                 p <- dyckPaths (m-1)                           ]+    | otherwise = [ bracket p ++ q | l <- [1..m-1] , p <- dyckPaths (l-1) , q <- worker (k-1) (m-l) ]+    where+      bracket p = UpStep : p ++ [DownStep] +++-- | There is a bijection from the set of non-empty Dyck paths of length @2n@ which touch the zero lines @t@ times,+-- to lattice paths from @(0,0)@ to @(2n-t-1,t-1)@ (just remove all the down-steps just before touching+-- the zero line, and also the very first up-step). This gives us a counting formula.+countTouchingDyckPaths +  :: Int   -- ^ @k@ = number of zero-touches+  -> Int   -- ^ @m@ = half-length+  -> Integer+countTouchingDyckPaths t n+  | t==0 && n==0   = 1+  | otherwise      = countLatticePaths (2*n-t-1,t-1)++--------------------------------------------------------------------------------+-- * Dyck paths with given number of peaks++-- | @peakingDyckPaths k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ with exactly @k@ peaks.+--+-- Synonym for 'peakingDyckPathsNaive'+--+peakingDyckPaths+  :: Int      -- ^ @k@ = number of peaks+  -> Int      -- ^ @m@ = half-length+  -> [LatticePath]+peakingDyckPaths = peakingDyckPathsNaive ++-- | @peakingDyckPathsNaive k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ with exactly @k@ peaks.+--+-- > sort (peakingDyckPathsNaive k m) = sort [ p | p <- dyckPaths m , pathNumberOfPeaks p == k ]+--  +-- Naive recursive algorithm, resulting order is pretty ad-hoc.+--+peakingDyckPathsNaive +  :: Int      -- ^ @k@ = number of peaks+  -> Int      -- ^ @m@ = half-length+  -> [LatticePath]+peakingDyckPathsNaive = worker where+  worker !k !m+    | m == 0    = if k==0 then [[]] else []       +    | k <= 0    = []+    | m <  0    = []+    | k == 1    = [ singlePeak m ] +    | otherwise = as ++ bs ++ cs+    where+      as = [ bracket p      |                                 p <- worker k (m-1)                           ]+      bs = [ smallHill ++ q |                                                       q <- worker (k-1) (m-1) ]+      cs = [ bracket p ++ q | l <- [2..m-1] , a <- [1..k-1] , p <- worker a (l-1) , q <- worker (k-a) (m-l) ]+      smallHill     = [ UpStep , DownStep ]+      singlePeak !m = replicate m UpStep ++ replicate m DownStep +      bracket p = UpStep : p ++ [DownStep] ++-- | Dyck paths of length @2m@ with @k@ peaks are counted by the Narayana numbers @N(m,k) = \binom{m}{k} \binom{m}{k-1} / m@+countPeakingDyckPaths+  :: Int      -- ^ @k@ = number of peaks+  -> Int      -- ^ @m@ = half-length+  -> Integer+countPeakingDyckPaths k m +  | m == 0    = if k==0 then 1 else 0+  | k <= 0    = 0+  | m <  0    = 0+  | k == 1    = 1+  | otherwise = div (binomial m k * binomial m (k-1)) (fromIntegral m)++--------------------------------------------------------------------------------+-- * Random lattice paths++-- | A uniformly random Dyck path of length @2m@+randomDyckPath :: RandomGen g => Int -> g -> (LatticePath,g)+randomDyckPath m g0 = (nestedParensToDyckPath parens, g1) where+  (parens,g1) = randomNestedParentheses m g0++--------------------------------------------------------------------------------+
Math/Combinat/Numbers/Series.hs view
@@ -233,7 +233,7 @@  -- | Power series expansion of @(1-Sqrt[1-4x])/(2x)@ (the coefficients are the Catalan numbers) dyckSeries :: Num a => [a]-dyckSeries = [ fromInteger (catalan i) | i<-[0..] ]+dyckSeries = [ fromInteger (catalan i) | i<-[(0::Int)..] ]  -------------------------------------------------------------------------------- -- * \"Coin\" series
Math/Combinat/Partitions/Integer.hs view
@@ -526,6 +526,9 @@           EnglishNotationCCW -> reverse $ fromPartition $ dualPartition part           FrenchNotation     -> reverse $ fromPartition $ part +instance DrawASCII Partition where+  ascii = asciiFerrersDiagram+ --------------------------------------------------------------------------------  {-
Math/Combinat/Partitions/Multiset.hs view
@@ -1,24 +1,24 @@-
--- | Partitions of a multiset
-module Math.Combinat.Partitions.Multiset where
-
---------------------------------------------------------------------------------
-
-import Data.Array.Unboxed
-import Data.List
-
-import Math.Combinat.Partitions.Vector
-
---------------------------------------------------------------------------------
-                              
--- | Partitions of a multiset. Internally, this uses the vector partition algorithm
-partitionMultiset :: (Eq a, Ord a) => [a] -> [[[a]]]
-partitionMultiset xs = parts where
-  parts = (map . map) (f . elems) temp
-  f ns = concat (zipWith replicate ns zs)
-  temp = fasc3B_algorithm_M counts
-  counts = map length ys
-  ys = group (sort xs) 
-  zs = map head ys
-
++-- | Partitions of a multiset+module Math.Combinat.Partitions.Multiset where++--------------------------------------------------------------------------------++import Data.Array.Unboxed+import Data.List++import Math.Combinat.Partitions.Vector++--------------------------------------------------------------------------------+                              +-- | Partitions of a multiset. Internally, this uses the vector partition algorithm+partitionMultiset :: (Eq a, Ord a) => [a] -> [[[a]]]+partitionMultiset xs = parts where+  parts = (map . map) (f . elems) temp+  f ns = concat (zipWith replicate ns zs)+  temp = fasc3B_algorithm_M counts+  counts = map length ys+  ys = group (sort xs) +  zs = map head ys+ --------------------------------------------------------------------------------
Math/Combinat/Partitions/NonCrossing.hs view
@@ -1,205 +1,205 @@-
--- | Non-crossing partitions.
---
--- See eg. <http://en.wikipedia.org/wiki/Noncrossing_partition>
---
--- Non-crossing partitions of the set @[1..n]@ are encoded as lists of lists
--- in standard form: Entries decreasing in each block  and blocks listed in increasing order of their first entries.
--- For example the partition in the diagram
---
--- <<svg/noncrossing.svg>>
---
--- is represented as
---
--- > NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]
---
-
-{-# LANGUAGE BangPatterns #-}
-module Math.Combinat.Partitions.NonCrossing where
-
---------------------------------------------------------------------------------
-
-import Control.Applicative
-
-import Data.List
-import Data.Ord
-
-import System.Random
-
-import Math.Combinat.Numbers
-import Math.Combinat.LatticePaths
-import Math.Combinat.Helper
-import Math.Combinat.Partitions.Set
-import Math.Combinat.Partitions ( HasNumberOfParts(..) )
-
---------------------------------------------------------------------------------
--- * The type of non-crossing partitions
-
--- | A non-crossing partition of the set @[1..n]@ in standard form: 
--- entries decreasing in each block  and blocks listed in increasing order of their first entries.
-newtype NonCrossing = NonCrossing [[Int]] deriving (Eq,Ord,Show,Read)
-
--- | Checks whether a set partition is noncrossing.
---
--- Implementation method: we convert to a Dyck path and then back again, and finally compare. 
--- Probably not very efficient, but should be better than a naive check for crosses...)
---
-_isNonCrossing :: [[Int]] -> Bool
-_isNonCrossing zzs0 = _isNonCrossingUnsafe (_standardizeNonCrossing zzs0)
-
--- | Warning: This function assumes the standard ordering!
-_isNonCrossingUnsafe :: [[Int]] -> Bool
-_isNonCrossingUnsafe zzs = 
-  case _nonCrossingPartitionToDyckPathMaybe zzs of
-    Nothing   -> False
-    Just dyck -> case dyckPathToNonCrossingPartitionMaybe dyck of
-      Nothing                -> False
-      Just (NonCrossing yys) -> yys == zzs
-
--- | Convert to standard form: entries decreasing in each block 
--- and blocks listed in increasing order of their first entries.
-_standardizeNonCrossing :: [[Int]] -> [[Int]]
-_standardizeNonCrossing = sortBy (comparing myhead) . map reverseSort where
-  myhead xs = case xs of
-    (x:xs) -> x
-    []     -> error "_standardizeNonCrossing: empty subset"
-
-fromNonCrossing :: NonCrossing -> [[Int]]
-fromNonCrossing (NonCrossing xs) = xs
-
-toNonCrossingUnsafe :: [[Int]] -> NonCrossing
-toNonCrossingUnsafe = NonCrossing
-
--- | Throws an error if the input is not a non-crossing partition
-toNonCrossing :: [[Int]] -> NonCrossing
-toNonCrossing xxs = case toNonCrossingMaybe xxs of
-  Just nc -> nc
-  Nothing -> error "toNonCrossing: not a non-crossing partition"
-
-toNonCrossingMaybe :: [[Int]] -> Maybe NonCrossing
-toNonCrossingMaybe xxs0 = 
-  if _isNonCrossingUnsafe xxs
-    then Just $ NonCrossing xxs
-    else Nothing
-  where 
-    xxs = _standardizeNonCrossing xxs0
-
--- | If a set partition is actually non-crossing, then we can convert it
-setPartitionToNonCrossing :: SetPartition -> Maybe NonCrossing
-setPartitionToNonCrossing (SetPartition zzs0) =
-  if _isNonCrossingUnsafe zzs
-    then Just $ NonCrossing zzs
-    else Nothing
-  where
-    zzs = _standardizeNonCrossing zzs0
-
-instance HasNumberOfParts NonCrossing where
-  numberOfParts (NonCrossing p) = length p
-
---------------------------------------------------------------------------------
--- * Bijection to Dyck paths
-
--- | Bijection between Dyck paths and noncrossing partitions
---
--- Based on: David Callan: /Sets, Lists and Noncrossing Partitions/
---
--- Fails if the input is not a Dyck path.
-dyckPathToNonCrossingPartition :: LatticePath -> NonCrossing
-dyckPathToNonCrossingPartition = NonCrossing . go 0 [] [] [] where
-  go :: Int -> [Int] -> [Int] -> [[Int]] -> LatticePath -> [[Int]] 
-  go !cnt stack small big path =
-    case path of
-      (x:xs) -> case x of 
-        UpStep   -> let cnt' = cnt + 1 in case xs of
-          (y:ys)   -> case y of
-            UpStep   -> go cnt' (cnt':stack) small                  big  xs  
-            DownStep -> go cnt' (cnt':stack) []    (reverse small : big) xs
-          []       -> error "dyckPathToNonCrossingPartition: last step is an UpStep (thus input was not a Dyck path)"
-        DownStep -> case stack of
-          (k:ks)   -> go cnt ks (k:small) big xs
-          []       -> error "dyckPathToNonCrossingPartition: empty stack, shouldn't happen (thus input was not a Dyck path)"
-      [] -> tail $ reverse (reverse small : big)
-
--- | Safe version of 'dyckPathToNonCrossingPartition'
-dyckPathToNonCrossingPartitionMaybe :: LatticePath -> Maybe NonCrossing
-dyckPathToNonCrossingPartitionMaybe = fmap NonCrossing . go 0 [] [] [] where
-  go :: Int -> [Int] -> [Int] -> [[Int]] -> LatticePath -> Maybe [[Int]] 
-  go !cnt stack small big path =
-    case path of
-      (x:xs) -> case x of 
-        UpStep   -> let cnt' = cnt + 1 in case xs of
-          (y:ys)   -> case y of
-            UpStep   -> go cnt' (cnt':stack) small                  big  xs  
-            DownStep -> go cnt' (cnt':stack) []    (reverse small : big) xs
-          []       -> Nothing
-        DownStep -> case stack of
-          (k:ks)   -> go cnt ks (k:small) big xs
-          []       -> Nothing
-      [] -> Just $ tail $ reverse (reverse small : big)
-
--- | The inverse bijection (should never fail proper 'NonCrossing'-s)
-nonCrossingPartitionToDyckPath :: NonCrossing -> LatticePath
-nonCrossingPartitionToDyckPath (NonCrossing zzs) = go 0 zzs where
-  go !k (ys@(y:_):yys) = replicate (y-k) UpStep ++ replicate (length ys) DownStep ++ go y yys
-  go !k []             = []
-  go _  _              = error "nonCrossingPartitionToDyckPath: shouldnt't happen"
-
--- | Safe version 'nonCrossingPartitionToDyckPath'
-_nonCrossingPartitionToDyckPathMaybe :: [[Int]] -> Maybe LatticePath
-_nonCrossingPartitionToDyckPathMaybe = go 0 where
-  go !k (ys@(y:_):yys) = fmap (\zs -> replicate (y-k) UpStep ++ replicate (length ys) DownStep ++ zs) (go y yys)
-  go !k []             = Just []
-  go _  _              = Nothing
-
---------------------------------------------------------------------------------
-
-{- 
--- this should be mapped to NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]
-testpath = [u,u,u,d,u,u,d,d,d,u,u,d,d,d,u,u,d,d] where
-  u = UpStep
-  d = DownStep
-
-testnc = NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]
--}
-
---------------------------------------------------------------------------------
--- * Generating non-crossing partitions
-
--- | Lists all non-crossing partitions of @[1..n]@
---
--- Equivalent to (but orders of magnitude faster than) filtering out the non-crossing ones:
---
--- > (sort $ catMaybes $ map setPartitionToNonCrossing $ setPartitions n) == sort (nonCrossingPartitions n)
---
-nonCrossingPartitions :: Int -> [NonCrossing]
-nonCrossingPartitions = map dyckPathToNonCrossingPartition . dyckPaths
-
--- | Lists all non-crossing partitions of @[1..n]@ into @k@ parts.
---
--- > sort (nonCrossingPartitionsWithKParts k n) == sort [ p | p <- nonCrossingPartitions n , numberOfParts p == k ]
---
-nonCrossingPartitionsWithKParts 
-  :: Int   -- ^ @k@ = number of parts 
-  -> Int   -- ^ @n@ = size of the set
-  -> [NonCrossing]
-nonCrossingPartitionsWithKParts k n = map dyckPathToNonCrossingPartition $ peakingDyckPaths k n
-
--- | Non-crossing partitions are counted by the Catalan numbers
-countNonCrossingPartitions :: Int -> Integer
-countNonCrossingPartitions = countDyckPaths
-
--- | Non-crossing partitions with @k@ parts are counted by the Naranaya numbers
-countNonCrossingPartitionsWithKParts 
-  :: Int   -- ^ @k@ = number of parts 
-  -> Int   -- ^ @n@ = size of the set
-  -> Integer
-countNonCrossingPartitionsWithKParts = countPeakingDyckPaths
-
---------------------------------------------------------------------------------
-
--- | Uniformly random non-crossing partition
-randomNonCrossingPartition :: RandomGen g => Int -> g -> (NonCrossing,g)
-randomNonCrossingPartition n g0 = (dyckPathToNonCrossingPartition dyck, g1) where
-  (dyck,g1) = randomDyckPath n g0
-
---------------------------------------------------------------------------------
++-- | Non-crossing partitions.+--+-- See eg. <http://en.wikipedia.org/wiki/Noncrossing_partition>+--+-- Non-crossing partitions of the set @[1..n]@ are encoded as lists of lists+-- in standard form: Entries decreasing in each block  and blocks listed in increasing order of their first entries.+-- For example the partition in the diagram+--+-- <<svg/noncrossing.svg>>+--+-- is represented as+--+-- > NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]+--++{-# LANGUAGE BangPatterns #-}+module Math.Combinat.Partitions.NonCrossing where++--------------------------------------------------------------------------------++import Control.Applicative++import Data.List+import Data.Ord++import System.Random++import Math.Combinat.Numbers+import Math.Combinat.LatticePaths+import Math.Combinat.Helper+import Math.Combinat.Partitions.Set+import Math.Combinat.Partitions ( HasNumberOfParts(..) )++--------------------------------------------------------------------------------+-- * The type of non-crossing partitions++-- | A non-crossing partition of the set @[1..n]@ in standard form: +-- entries decreasing in each block  and blocks listed in increasing order of their first entries.+newtype NonCrossing = NonCrossing [[Int]] deriving (Eq,Ord,Show,Read)++-- | Checks whether a set partition is noncrossing.+--+-- Implementation method: we convert to a Dyck path and then back again, and finally compare. +-- Probably not very efficient, but should be better than a naive check for crosses...)+--+_isNonCrossing :: [[Int]] -> Bool+_isNonCrossing zzs0 = _isNonCrossingUnsafe (_standardizeNonCrossing zzs0)++-- | Warning: This function assumes the standard ordering!+_isNonCrossingUnsafe :: [[Int]] -> Bool+_isNonCrossingUnsafe zzs = +  case _nonCrossingPartitionToDyckPathMaybe zzs of+    Nothing   -> False+    Just dyck -> case dyckPathToNonCrossingPartitionMaybe dyck of+      Nothing                -> False+      Just (NonCrossing yys) -> yys == zzs++-- | Convert to standard form: entries decreasing in each block +-- and blocks listed in increasing order of their first entries.+_standardizeNonCrossing :: [[Int]] -> [[Int]]+_standardizeNonCrossing = sortBy (comparing myhead) . map reverseSort where+  myhead xs = case xs of+    (x:xs) -> x+    []     -> error "_standardizeNonCrossing: empty subset"++fromNonCrossing :: NonCrossing -> [[Int]]+fromNonCrossing (NonCrossing xs) = xs++toNonCrossingUnsafe :: [[Int]] -> NonCrossing+toNonCrossingUnsafe = NonCrossing++-- | Throws an error if the input is not a non-crossing partition+toNonCrossing :: [[Int]] -> NonCrossing+toNonCrossing xxs = case toNonCrossingMaybe xxs of+  Just nc -> nc+  Nothing -> error "toNonCrossing: not a non-crossing partition"++toNonCrossingMaybe :: [[Int]] -> Maybe NonCrossing+toNonCrossingMaybe xxs0 = +  if _isNonCrossingUnsafe xxs+    then Just $ NonCrossing xxs+    else Nothing+  where +    xxs = _standardizeNonCrossing xxs0++-- | If a set partition is actually non-crossing, then we can convert it+setPartitionToNonCrossing :: SetPartition -> Maybe NonCrossing+setPartitionToNonCrossing (SetPartition zzs0) =+  if _isNonCrossingUnsafe zzs+    then Just $ NonCrossing zzs+    else Nothing+  where+    zzs = _standardizeNonCrossing zzs0++instance HasNumberOfParts NonCrossing where+  numberOfParts (NonCrossing p) = length p++--------------------------------------------------------------------------------+-- * Bijection to Dyck paths++-- | Bijection between Dyck paths and noncrossing partitions+--+-- Based on: David Callan: /Sets, Lists and Noncrossing Partitions/+--+-- Fails if the input is not a Dyck path.+dyckPathToNonCrossingPartition :: LatticePath -> NonCrossing+dyckPathToNonCrossingPartition = NonCrossing . go 0 [] [] [] where+  go :: Int -> [Int] -> [Int] -> [[Int]] -> LatticePath -> [[Int]] +  go !cnt stack small big path =+    case path of+      (x:xs) -> case x of +        UpStep   -> let cnt' = cnt + 1 in case xs of+          (y:ys)   -> case y of+            UpStep   -> go cnt' (cnt':stack) small                  big  xs  +            DownStep -> go cnt' (cnt':stack) []    (reverse small : big) xs+          []       -> error "dyckPathToNonCrossingPartition: last step is an UpStep (thus input was not a Dyck path)"+        DownStep -> case stack of+          (k:ks)   -> go cnt ks (k:small) big xs+          []       -> error "dyckPathToNonCrossingPartition: empty stack, shouldn't happen (thus input was not a Dyck path)"+      [] -> tail $ reverse (reverse small : big)++-- | Safe version of 'dyckPathToNonCrossingPartition'+dyckPathToNonCrossingPartitionMaybe :: LatticePath -> Maybe NonCrossing+dyckPathToNonCrossingPartitionMaybe = fmap NonCrossing . go 0 [] [] [] where+  go :: Int -> [Int] -> [Int] -> [[Int]] -> LatticePath -> Maybe [[Int]] +  go !cnt stack small big path =+    case path of+      (x:xs) -> case x of +        UpStep   -> let cnt' = cnt + 1 in case xs of+          (y:ys)   -> case y of+            UpStep   -> go cnt' (cnt':stack) small                  big  xs  +            DownStep -> go cnt' (cnt':stack) []    (reverse small : big) xs+          []       -> Nothing+        DownStep -> case stack of+          (k:ks)   -> go cnt ks (k:small) big xs+          []       -> Nothing+      [] -> Just $ tail $ reverse (reverse small : big)++-- | The inverse bijection (should never fail proper 'NonCrossing'-s)+nonCrossingPartitionToDyckPath :: NonCrossing -> LatticePath+nonCrossingPartitionToDyckPath (NonCrossing zzs) = go 0 zzs where+  go !k (ys@(y:_):yys) = replicate (y-k) UpStep ++ replicate (length ys) DownStep ++ go y yys+  go !k []             = []+  go _  _              = error "nonCrossingPartitionToDyckPath: shouldnt't happen"++-- | Safe version 'nonCrossingPartitionToDyckPath'+_nonCrossingPartitionToDyckPathMaybe :: [[Int]] -> Maybe LatticePath+_nonCrossingPartitionToDyckPathMaybe = go 0 where+  go !k (ys@(y:_):yys) = fmap (\zs -> replicate (y-k) UpStep ++ replicate (length ys) DownStep ++ zs) (go y yys)+  go !k []             = Just []+  go _  _              = Nothing++--------------------------------------------------------------------------------++{- +-- this should be mapped to NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]+testpath = [u,u,u,d,u,u,d,d,d,u,u,d,d,d,u,u,d,d] where+  u = UpStep+  d = DownStep++testnc = NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]+-}++--------------------------------------------------------------------------------+-- * Generating non-crossing partitions++-- | Lists all non-crossing partitions of @[1..n]@+--+-- Equivalent to (but orders of magnitude faster than) filtering out the non-crossing ones:+--+-- > (sort $ catMaybes $ map setPartitionToNonCrossing $ setPartitions n) == sort (nonCrossingPartitions n)+--+nonCrossingPartitions :: Int -> [NonCrossing]+nonCrossingPartitions = map dyckPathToNonCrossingPartition . dyckPaths++-- | Lists all non-crossing partitions of @[1..n]@ into @k@ parts.+--+-- > sort (nonCrossingPartitionsWithKParts k n) == sort [ p | p <- nonCrossingPartitions n , numberOfParts p == k ]+--+nonCrossingPartitionsWithKParts +  :: Int   -- ^ @k@ = number of parts +  -> Int   -- ^ @n@ = size of the set+  -> [NonCrossing]+nonCrossingPartitionsWithKParts k n = map dyckPathToNonCrossingPartition $ peakingDyckPaths k n++-- | Non-crossing partitions are counted by the Catalan numbers+countNonCrossingPartitions :: Int -> Integer+countNonCrossingPartitions = countDyckPaths++-- | Non-crossing partitions with @k@ parts are counted by the Naranaya numbers+countNonCrossingPartitionsWithKParts +  :: Int   -- ^ @k@ = number of parts +  -> Int   -- ^ @n@ = size of the set+  -> Integer+countNonCrossingPartitionsWithKParts = countPeakingDyckPaths++--------------------------------------------------------------------------------++-- | Uniformly random non-crossing partition+randomNonCrossingPartition :: RandomGen g => Int -> g -> (NonCrossing,g)+randomNonCrossingPartition n g0 = (dyckPathToNonCrossingPartition dyck, g1) where+  (dyck,g1) = randomDyckPath n g0++--------------------------------------------------------------------------------
Math/Combinat/Partitions/Plane.hs view
@@ -1,116 +1,116 @@-
--- | Plane partitions. See eg. <http://en.wikipedia.org/wiki/Plane_partition>
---
--- Plane partitions are encoded as lists of lists of Z heights. For example the plane 
--- partition in the picture
--- 
--- <<svg/plane_partition.svg>>
---
--- is encoded as
---
--- > PlanePart [ [5,4,3,3,1]
--- >           , [4,4,2,1]
--- >           , [3,2]
--- >           , [2,1]
--- >           , [1]
--- >           , [1]
--- >           ]
--- 
-{-# LANGUAGE BangPatterns #-}
-module Math.Combinat.Partitions.Plane where
-
---------------------------------------------------------------------------------
-
-import Data.List
-import Data.Array
-
-import Math.Combinat.Partitions
-import Math.Combinat.Tableaux as Tableaux
-import Math.Combinat.Helper
-
---------------------------------------------------------------------------------
--- * the type of plane partitions
-
--- | A plane partition encoded as a tablaeu (the \"Z\" heights are the numbers)
-newtype PlanePart = PlanePart [[Int]] deriving (Eq,Ord,Show)
-
-fromPlanePart :: PlanePart -> [[Int]]
-fromPlanePart (PlanePart xs) = xs
-
-isValidPlanePart :: [[Int]] -> Bool
-isValidPlanePart pps = 
-  and [ table!(i,j) >= table!(i  ,j+1) &&
-        table!(i,j) >= table!(i+1,j  )
-      | i<-[0..y-1] , j<-[0..x-1] 
-      ]
-  where
-    table :: Array (Int,Int) Int
-    table = accumArray const 0 ((0,0),(y,x)) [ ((i,j),k) | (i,ps) <- zip [0..] pps , (j,k) <- zip [0..] ps ]
-    y = length pps
-    x = maximum (map length pps)
-
--- | Throws an exception if the input is not a plane partition
-toPlanePart :: [[Int]] -> PlanePart
-toPlanePart pps = if isValidPlanePart pps
-  then PlanePart $ filter (not . null) $ map (filter (>0)) $ pps
-  else error "toPlanePart: not a plane partition"
-
--- | The XY projected shape of a plane partition, as an integer partition
-planePartShape :: PlanePart -> Partition
-planePartShape = Tableaux.shape . fromPlanePart
-
--- | The Z height of a plane partition
-planePartZHeight :: PlanePart -> Int
-planePartZHeight (PlanePart xs) = 
-  case xs of
-    ((h:_):_) -> h
-    _         -> 0
-
-planePartWeight :: PlanePart -> Int
-planePartWeight (PlanePart xs) = sum' (map sum' xs)
-
---------------------------------------------------------------------------------
--- * constructing plane partitions
-
-singleLayer :: Partition -> PlanePart
-singleLayer = PlanePart . map (\k -> replicate k 1) . fromPartition 
-
--- |  Stacks layers of partitions into a plane partition.
--- Throws an exception if they do not form a plane partition.
-stackLayers :: [Partition] -> PlanePart
-stackLayers layers = if and [ isSubPartitionOf p q | (q,p) <- pairs layers ]
-  then unsafeStackLayers layers
-  else error "stackLayers: the layers do not form a plane partition"
-
--- | Stacks layers of partitions into a plane partition.
--- This is unsafe in the sense that we don't check that the partitions fit on the top of each other.
-unsafeStackLayers :: [Partition] -> PlanePart
-unsafeStackLayers []            = PlanePart []
-unsafeStackLayers (bottom:rest) = PlanePart $ foldl addLayer (fromPlanePart $ singleLayer bottom) rest where
-  addLayer :: [[Int]] -> Partition -> [[Int]]
-  addLayer xxs (Partition ps) = [ zipWith (+) xs (replicate p 1 ++ repeat 0) | (xs,p) <- zip xxs (ps ++ repeat 0) ] 
-
--- | The \"layers\" of a plane partition (in direction @Z@). We should have
---
--- > unsafeStackLayers (planePartLayers pp) == pp
--- 
-planePartLayers :: PlanePart -> [Partition]
-planePartLayers pp@(PlanePart xs) = [ layer h | h<-[1..planePartZHeight pp] ] where
-  layer h = Partition $ filter (>0) $ map sum' $ (map . map) (f h) xs
-  f h = \k -> if k>=h then 1 else 0
-
---------------------------------------------------------------------------------
--- * generating plane partitions
-
--- | Plane partitions of a given weight
-planePartitions :: Int -> [PlanePart]
-planePartitions d 
-  | d <  0     = []
-  | d == 0     = [PlanePart []]
-  | otherwise  = concat [ go (d-n) [p] | n<-[1..d] , p<-partitions n ]
-  where
-    go :: Int -> [Partition] -> [PlanePart]
-    go  0   acc       = [unsafeStackLayers (reverse acc)]
-    go !rem acc@(h:_) = concat [ go (rem-k) (this:acc) | k<-[1..rem] , this <- subPartitions k h ]
-
---------------------------------------------------------------------------------
++-- | Plane partitions. See eg. <http://en.wikipedia.org/wiki/Plane_partition>+--+-- Plane partitions are encoded as lists of lists of Z heights. For example the plane +-- partition in the picture+-- +-- <<svg/plane_partition.svg>>+--+-- is encoded as+--+-- > PlanePart [ [5,4,3,3,1]+-- >           , [4,4,2,1]+-- >           , [3,2]+-- >           , [2,1]+-- >           , [1]+-- >           , [1]+-- >           ]+-- +{-# LANGUAGE BangPatterns #-}+module Math.Combinat.Partitions.Plane where++--------------------------------------------------------------------------------++import Data.List+import Data.Array++import Math.Combinat.Partitions+import Math.Combinat.Tableaux as Tableaux+import Math.Combinat.Helper++--------------------------------------------------------------------------------+-- * the type of plane partitions++-- | A plane partition encoded as a tablaeu (the \"Z\" heights are the numbers)+newtype PlanePart = PlanePart [[Int]] deriving (Eq,Ord,Show)++fromPlanePart :: PlanePart -> [[Int]]+fromPlanePart (PlanePart xs) = xs++isValidPlanePart :: [[Int]] -> Bool+isValidPlanePart pps = +  and [ table!(i,j) >= table!(i  ,j+1) &&+        table!(i,j) >= table!(i+1,j  )+      | i<-[0..y-1] , j<-[0..x-1] +      ]+  where+    table :: Array (Int,Int) Int+    table = accumArray const 0 ((0,0),(y,x)) [ ((i,j),k) | (i,ps) <- zip [0..] pps , (j,k) <- zip [0..] ps ]+    y = length pps+    x = maximum (map length pps)++-- | Throws an exception if the input is not a plane partition+toPlanePart :: [[Int]] -> PlanePart+toPlanePart pps = if isValidPlanePart pps+  then PlanePart $ filter (not . null) $ map (filter (>0)) $ pps+  else error "toPlanePart: not a plane partition"++-- | The XY projected shape of a plane partition, as an integer partition+planePartShape :: PlanePart -> Partition+planePartShape = Tableaux.shape . fromPlanePart++-- | The Z height of a plane partition+planePartZHeight :: PlanePart -> Int+planePartZHeight (PlanePart xs) = +  case xs of+    ((h:_):_) -> h+    _         -> 0++planePartWeight :: PlanePart -> Int+planePartWeight (PlanePart xs) = sum' (map sum' xs)++--------------------------------------------------------------------------------+-- * constructing plane partitions++singleLayer :: Partition -> PlanePart+singleLayer = PlanePart . map (\k -> replicate k 1) . fromPartition ++-- |  Stacks layers of partitions into a plane partition.+-- Throws an exception if they do not form a plane partition.+stackLayers :: [Partition] -> PlanePart+stackLayers layers = if and [ isSubPartitionOf p q | (q,p) <- pairs layers ]+  then unsafeStackLayers layers+  else error "stackLayers: the layers do not form a plane partition"++-- | Stacks layers of partitions into a plane partition.+-- This is unsafe in the sense that we don't check that the partitions fit on the top of each other.+unsafeStackLayers :: [Partition] -> PlanePart+unsafeStackLayers []            = PlanePart []+unsafeStackLayers (bottom:rest) = PlanePart $ foldl addLayer (fromPlanePart $ singleLayer bottom) rest where+  addLayer :: [[Int]] -> Partition -> [[Int]]+  addLayer xxs (Partition ps) = [ zipWith (+) xs (replicate p 1 ++ repeat 0) | (xs,p) <- zip xxs (ps ++ repeat 0) ] ++-- | The \"layers\" of a plane partition (in direction @Z@). We should have+--+-- > unsafeStackLayers (planePartLayers pp) == pp+-- +planePartLayers :: PlanePart -> [Partition]+planePartLayers pp@(PlanePart xs) = [ layer h | h<-[1..planePartZHeight pp] ] where+  layer h = Partition $ filter (>0) $ map sum' $ (map . map) (f h) xs+  f h = \k -> if k>=h then 1 else 0++--------------------------------------------------------------------------------+-- * generating plane partitions++-- | Plane partitions of a given weight+planePartitions :: Int -> [PlanePart]+planePartitions d +  | d <  0     = []+  | d == 0     = [PlanePart []]+  | otherwise  = concat [ go (d-n) [p] | n<-[1..d] , p<-partitions n ]+  where+    go :: Int -> [Partition] -> [PlanePart]+    go  0   acc       = [unsafeStackLayers (reverse acc)]+    go !rem acc@(h:_) = concat [ go (rem-k) (this:acc) | k<-[1..rem] , this <- subPartitions k h ]++--------------------------------------------------------------------------------
Math/Combinat/Partitions/Set.hs view
@@ -1,99 +1,99 @@-
--- | Set partitions.
---
--- See eg. <http://en.wikipedia.org/wiki/Partition_of_a_set>
--- 
-
-{-# LANGUAGE BangPatterns #-}
-module Math.Combinat.Partitions.Set where
-
---------------------------------------------------------------------------------
-
-import Data.List
-import Data.Ord
-
-import System.Random
-
-import Math.Combinat.Sets
-import Math.Combinat.Numbers
-import Math.Combinat.Helper
-import Math.Combinat.Partitions ( HasNumberOfParts(..) )
-
---------------------------------------------------------------------------------
--- * The type of set partitions
-
--- | A partition of the set @[1..n]@ (in standard order)
-newtype SetPartition = SetPartition [[Int]] deriving (Eq,Ord,Show,Read)
-
-_standardizeSetPartition :: [[Int]] -> [[Int]]
-_standardizeSetPartition = sortBy (comparing myhead) . map sort where
-  myhead xs = case xs of
-    (x:xs) -> x
-    []     -> error "_standardizeSetPartition: empty subset"
-
-fromSetPartition :: SetPartition -> [[Int]]
-fromSetPartition (SetPartition zzs) = zzs
-
-toSetPartitionUnsafe :: [[Int]] -> SetPartition
-toSetPartitionUnsafe = SetPartition
-
-toSetPartition :: [[Int]] -> SetPartition
-toSetPartition zzs = if _isSetPartition zzs
-  then SetPartition (_standardizeSetPartition zzs)
-  else error "toSetPartition: not a set partition"
-
-_isSetPartition :: [[Int]] -> Bool
-_isSetPartition zzs = sort (concat zzs) == [1..n] where 
-  n = sum' (map length zzs)
-
-instance HasNumberOfParts SetPartition where
-  numberOfParts (SetPartition p) = length p
-
---------------------------------------------------------------------------------
--- * Generating set partitions
-
--- | Synonym for 'setPartitionsNaive'
-setPartitions :: Int -> [SetPartition]
-setPartitions = setPartitionsNaive
-
--- | Synonym for 'setPartitionsWithKPartsNaive'
---
--- > sort (setPartitionsWithKParts k n) == sort [ p | p <- setPartitions n , numberOfParts p == k ]
--- 
-setPartitionsWithKParts   
-  :: Int    -- ^ @k@ = number of parts
-  -> Int    -- ^ @n@ = size of the set
-  -> [SetPartition]
-setPartitionsWithKParts = setPartitionsWithKPartsNaive
-
--- | List all set partitions of @[1..n]@, naive algorithm
-setPartitionsNaive :: Int -> [SetPartition]
-setPartitionsNaive n = map (SetPartition . _standardizeSetPartition) $ go [1..n] where
-  go :: [Int] -> [[[Int]]]
-  go []     = [[]]
-  go (z:zs) = [ s : rest | k <- [1..n] , s0 <- choose (k-1) zs , let s = z:s0 , rest <- go (zs \\ s) ]
-
--- | Set partitions of the set @[1..n]@ into @k@ parts
-setPartitionsWithKPartsNaive 
-  :: Int    -- ^ @k@ = number of parts
-  -> Int    -- ^ @n@ = size of the set
-  -> [SetPartition]
-setPartitionsWithKPartsNaive k n = map (SetPartition . _standardizeSetPartition) $ go k [1..n] where
-  go :: Int -> [Int] -> [[[Int]]]
-  go !k []     = if k==0 then [[]] else []
-  go  1 zs     = [[zs]]
-  go !k (z:zs) = [ s : rest | l <- [1..n] , s0 <- choose (l-1) zs , let s = z:s0 , rest <- go (k-1) (zs \\ s) ]
-
-
--- | Set partitions are counted by the Bell numbers
-countSetPartitions :: Int -> Integer
-countSetPartitions = bellNumber 
-
--- | Set partitions of size @k@ are counted by the Stirling numbers of second kind
-countSetPartitionsWithKParts 
-  :: Int    -- ^ @k@ = number of parts
-  -> Int    -- ^ @n@ = size of the set
-  -> Integer
-countSetPartitionsWithKParts k n = stirling2nd n k
-
---------------------------------------------------------------------------------
++-- | Set partitions.+--+-- See eg. <http://en.wikipedia.org/wiki/Partition_of_a_set>+-- ++{-# LANGUAGE BangPatterns #-}+module Math.Combinat.Partitions.Set where++--------------------------------------------------------------------------------++import Data.List+import Data.Ord++import System.Random++import Math.Combinat.Sets+import Math.Combinat.Numbers+import Math.Combinat.Helper+import Math.Combinat.Partitions ( HasNumberOfParts(..) )++--------------------------------------------------------------------------------+-- * The type of set partitions++-- | A partition of the set @[1..n]@ (in standard order)+newtype SetPartition = SetPartition [[Int]] deriving (Eq,Ord,Show,Read)++_standardizeSetPartition :: [[Int]] -> [[Int]]+_standardizeSetPartition = sortBy (comparing myhead) . map sort where+  myhead xs = case xs of+    (x:xs) -> x+    []     -> error "_standardizeSetPartition: empty subset"++fromSetPartition :: SetPartition -> [[Int]]+fromSetPartition (SetPartition zzs) = zzs++toSetPartitionUnsafe :: [[Int]] -> SetPartition+toSetPartitionUnsafe = SetPartition++toSetPartition :: [[Int]] -> SetPartition+toSetPartition zzs = if _isSetPartition zzs+  then SetPartition (_standardizeSetPartition zzs)+  else error "toSetPartition: not a set partition"++_isSetPartition :: [[Int]] -> Bool+_isSetPartition zzs = sort (concat zzs) == [1..n] where +  n = sum' (map length zzs)++instance HasNumberOfParts SetPartition where+  numberOfParts (SetPartition p) = length p++--------------------------------------------------------------------------------+-- * Generating set partitions++-- | Synonym for 'setPartitionsNaive'+setPartitions :: Int -> [SetPartition]+setPartitions = setPartitionsNaive++-- | Synonym for 'setPartitionsWithKPartsNaive'+--+-- > sort (setPartitionsWithKParts k n) == sort [ p | p <- setPartitions n , numberOfParts p == k ]+-- +setPartitionsWithKParts   +  :: Int    -- ^ @k@ = number of parts+  -> Int    -- ^ @n@ = size of the set+  -> [SetPartition]+setPartitionsWithKParts = setPartitionsWithKPartsNaive++-- | List all set partitions of @[1..n]@, naive algorithm+setPartitionsNaive :: Int -> [SetPartition]+setPartitionsNaive n = map (SetPartition . _standardizeSetPartition) $ go [1..n] where+  go :: [Int] -> [[[Int]]]+  go []     = [[]]+  go (z:zs) = [ s : rest | k <- [1..n] , s0 <- choose (k-1) zs , let s = z:s0 , rest <- go (zs \\ s) ]++-- | Set partitions of the set @[1..n]@ into @k@ parts+setPartitionsWithKPartsNaive +  :: Int    -- ^ @k@ = number of parts+  -> Int    -- ^ @n@ = size of the set+  -> [SetPartition]+setPartitionsWithKPartsNaive k n = map (SetPartition . _standardizeSetPartition) $ go k [1..n] where+  go :: Int -> [Int] -> [[[Int]]]+  go !k []     = if k==0 then [[]] else []+  go  1 zs     = [[zs]]+  go !k (z:zs) = [ s : rest | l <- [1..n] , s0 <- choose (l-1) zs , let s = z:s0 , rest <- go (k-1) (zs \\ s) ]+++-- | Set partitions are counted by the Bell numbers+countSetPartitions :: Int -> Integer+countSetPartitions = bellNumber ++-- | Set partitions of size @k@ are counted by the Stirling numbers of second kind+countSetPartitionsWithKParts +  :: Int    -- ^ @k@ = number of parts+  -> Int    -- ^ @n@ = size of the set+  -> Integer+countSetPartitionsWithKParts k n = stirling2nd n k++--------------------------------------------------------------------------------
Math/Combinat/Partitions/Skew.hs view
@@ -1,67 +1,82 @@-
--- | Skew partitions.
---
--- Skew partitions are the difference of two integer partitions, denoted by @lambda/mu@.
---
-
-{-# LANGUAGE BangPatterns #-}
-module Math.Combinat.Partitions.Skew where
-
---------------------------------------------------------------------------------
-
-import Data.List
-
-import Math.Combinat.Partitions.Integer
-import Math.Combinat.ASCII
-
---------------------------------------------------------------------------------
-
--- | A skew partition @lambda/mu@ is 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@
-mkSkewPartition :: (Partition,Partition) -> SkewPartition
-mkSkewPartition ( lam@(Partition bs) , mu@(Partition as)) = if mu `isSubPartitionOf` lam 
-  then SkewPartition $ zipWith (\b a -> (a,b-a)) bs (as ++ repeat 0)
-  else error "mkSkewPartition: mu should be a subpartition of lambda!" 
-
--- | This function \"cuts off\" the \"uninteresting parts\" of a skew partition
-normalizeSkewPartition :: SkewPartition -> SkewPartition
-normalizeSkewPartition (SkewPartition abs) = SkewPartition abs' where
-  (as,bs) = unzip abs
-  a0 = minimum as
-  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
-fromSkewPartition :: SkewPartition -> (Partition,Partition)
-fromSkewPartition (SkewPartition list) = (toPartition (zipWith (+) as bs) , toPartition (filter (>0) as)) where
-  (as,bs) = unzip list
-
-outerPartition :: SkewPartition -> Partition  
-outerPartition = fst . fromSkewPartition 
-
-innerPartition :: SkewPartition -> Partition  
-innerPartition = snd . fromSkewPartition 
-
-dualSkewPartition :: SkewPartition -> SkewPartition
-dualSkewPartition = mkSkewPartition . f . fromSkewPartition where
-  f (lam,mu) = (dualPartition lam, dualPartition mu)
-
-asciiSkewFerrersDiagram :: SkewPartition -> ASCII
-asciiSkewFerrersDiagram = asciiSkewFerrersDiagram' ('@','.') EnglishNotation
-
-asciiSkewFerrersDiagram' 
-  :: (Char,Char)       
-  -> PartitionConvention -- Orientation
-  -> SkewPartition 
-  -> ASCII
-asciiSkewFerrersDiagram' (outer,inner) orient (SkewPartition abs) = asciiFromLines stuff where
-  stuff = case orient of
-    EnglishNotation    -> ls
-    EnglishNotationCCW -> reverse (transpose ls)
-    FrenchNotation     -> reverse ls
-  ls = [ replicate a inner ++ replicate b outer | (a,b) <- abs ]
-    
-  
---------------------------------------------------------------------------------
++-- | Skew partitions.+--+-- Skew partitions are the difference of two integer partitions, denoted by @lambda/mu@.+--++{-# LANGUAGE BangPatterns #-}+module Math.Combinat.Partitions.Skew where++--------------------------------------------------------------------------------++import Data.List++import Math.Combinat.Partitions.Integer+import Math.Combinat.ASCII++--------------------------------------------------------------------------------++-- | A skew partition @lambda/mu@ is 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@.+-- Throws an error if @mu@ is not a sub-partition of @lambda@.+mkSkewPartition :: (Partition,Partition) -> SkewPartition+mkSkewPartition ( lam@(Partition bs) , mu@(Partition as)) = if mu `isSubPartitionOf` lam +  then SkewPartition $ zipWith (\b a -> (a,b-a)) bs (as ++ repeat 0)+  else error "mkSkewPartition: mu should be a subpartition of lambda!" ++-- | Returns 'Nothing' if @mu@ is not a sub-partition of @lambda@.+safeSkewPartition :: (Partition,Partition) -> Maybe SkewPartition+safeSkewPartition ( lam@(Partition bs) , mu@(Partition as)) = if mu `isSubPartitionOf` lam +  then Just $ SkewPartition $ zipWith (\b a -> (a,b-a)) bs (as ++ repeat 0)+  else Nothing++skewPartitionWeight :: SkewPartition -> Int+skewPartitionWeight (SkewPartition abs) = foldl' (+) 0 (map snd abs)++-- | This function \"cuts off\" the \"uninteresting parts\" of a skew partition+normalizeSkewPartition :: SkewPartition -> SkewPartition+normalizeSkewPartition (SkewPartition abs) = SkewPartition abs' where+  (as,bs) = unzip abs+  a0 = minimum as+  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+fromSkewPartition :: SkewPartition -> (Partition,Partition)+fromSkewPartition (SkewPartition list) = (toPartition (zipWith (+) as bs) , toPartition (filter (>0) as)) where+  (as,bs) = unzip list++outerPartition :: SkewPartition -> Partition  +outerPartition = fst . fromSkewPartition ++innerPartition :: SkewPartition -> Partition  +innerPartition = snd . fromSkewPartition ++dualSkewPartition :: SkewPartition -> SkewPartition+dualSkewPartition = mkSkewPartition . f . fromSkewPartition where+  f (lam,mu) = (dualPartition lam, dualPartition mu)++--------------------------------------------------------------------------------++asciiSkewFerrersDiagram :: SkewPartition -> ASCII+asciiSkewFerrersDiagram = asciiSkewFerrersDiagram' ('@','.') EnglishNotation++asciiSkewFerrersDiagram' +  :: (Char,Char)       +  -> PartitionConvention -- Orientation+  -> SkewPartition +  -> ASCII+asciiSkewFerrersDiagram' (outer,inner) orient (SkewPartition abs) = asciiFromLines stuff where+  stuff = case orient of+    EnglishNotation    -> ls+    EnglishNotationCCW -> reverse (transpose ls)+    FrenchNotation     -> reverse ls+  ls = [ replicate a inner ++ replicate b outer | (a,b) <- abs ]++instance DrawASCII SkewPartition where+  ascii = asciiSkewFerrersDiagram     ++--------------------------------------------------------------------------------+
Math/Combinat/Partitions/Vector.hs view
@@ -1,82 +1,82 @@-
--- | Vector partitions. See:
---
---  * Donald E. Knuth: The Art of Computer Programming, vol 4, pre-fascicle 3B.
---
-
-{-# LANGUAGE BangPatterns #-}
-module Math.Combinat.Partitions.Vector where
-
---------------------------------------------------------------------------------
-
-import Data.Array.Unboxed
-import Data.List
-
---------------------------------------------------------------------------------
-
--- | Integer vectors. The indexing starts from 1.
-type IntVector = UArray Int Int
-
--- | Vector partitions. Basically a synonym for 'fasc3B_algorithm_M'.
-vectorPartitions :: IntVector -> [[IntVector]]
-vectorPartitions = fasc3B_algorithm_M . elems
-
-_vectorPartitions :: [Int] -> [[[Int]]]
-_vectorPartitions = map (map elems) . fasc3B_algorithm_M
-
--- | Generates all vector partitions 
---   (\"algorithm M\" in Knuth). 
---   The order is decreasing lexicographic.  
-fasc3B_algorithm_M :: [Int] -> [[IntVector]] 
-{- note to self: Knuth's descriptions of algorithms are still totally unreadable -}
-fasc3B_algorithm_M xs = worker [start] where
-
-  -- n = sum xs
-  m = length xs
-
-  start = [ (j,x,x) | (j,x) <- zip [1..] xs ]  
-  
-  worker stack@(last:_) = 
-    case decrease stack' of
-      Nothing -> [visited]
-      Just stack'' -> visited : worker stack''
-    where
-      stack'  = subtract_rec stack
-      visited = map to_vector stack'
-      
-  decrease (last:rest) = 
-    case span (\(_,_,v) -> v==0) (reverse last) of
-      ( _ , [(_,_,1)] ) -> case rest of
-        [] -> Nothing
-        _  -> decrease rest
-      ( second , (c,u,v):first ) -> Just (modified:rest) where 
-        modified =   
-          reverse first ++ 
-          (c,u,v-1) :  
-          [ (c,u,u) | (c,u,_) <- reverse second ] 
-      _ -> error "fasc3B_algorithm_M: should not happen"
-        
-  to_vector cuvs = 
-    accumArray (flip const) 0 (1,m)
-      [ (c,v) | (c,_,v) <- cuvs ] 
-
-  subtract_rec all@(last:_) = 
-    case subtract last of 
-      []  -> all
-      new -> subtract_rec (new:all) 
-
-  subtract [] = []
-  subtract full@((c,u,v):rest) = 
-    if w >= v 
-      then (c,w,v) : subtract   rest
-      else           subtract_b full
-    where w = u - v
-    
-  subtract_b [] = []
-  subtract_b ((c,u,v):rest) = 
-    if w /= 0 
-      then (c,w,w) : subtract_b rest
-      else           subtract_b rest
-    where w = u - v
-
---------------------------------------------------------------------------------
++-- | Vector partitions. See:+--+--  * Donald E. Knuth: The Art of Computer Programming, vol 4, pre-fascicle 3B.+--++{-# LANGUAGE BangPatterns #-}+module Math.Combinat.Partitions.Vector where++--------------------------------------------------------------------------------++import Data.Array.Unboxed+import Data.List++--------------------------------------------------------------------------------++-- | Integer vectors. The indexing starts from 1.+type IntVector = UArray Int Int++-- | Vector partitions. Basically a synonym for 'fasc3B_algorithm_M'.+vectorPartitions :: IntVector -> [[IntVector]]+vectorPartitions = fasc3B_algorithm_M . elems++_vectorPartitions :: [Int] -> [[[Int]]]+_vectorPartitions = map (map elems) . fasc3B_algorithm_M++-- | Generates all vector partitions +--   (\"algorithm M\" in Knuth). +--   The order is decreasing lexicographic.  +fasc3B_algorithm_M :: [Int] -> [[IntVector]] +{- note to self: Knuth's descriptions of algorithms are still totally unreadable -}+fasc3B_algorithm_M xs = worker [start] where++  -- n = sum xs+  m = length xs++  start = [ (j,x,x) | (j,x) <- zip [1..] xs ]  +  +  worker stack@(last:_) = +    case decrease stack' of+      Nothing -> [visited]+      Just stack'' -> visited : worker stack''+    where+      stack'  = subtract_rec stack+      visited = map to_vector stack'+      +  decrease (last:rest) = +    case span (\(_,_,v) -> v==0) (reverse last) of+      ( _ , [(_,_,1)] ) -> case rest of+        [] -> Nothing+        _  -> decrease rest+      ( second , (c,u,v):first ) -> Just (modified:rest) where +        modified =   +          reverse first ++ +          (c,u,v-1) :  +          [ (c,u,u) | (c,u,_) <- reverse second ] +      _ -> error "fasc3B_algorithm_M: should not happen"+        +  to_vector cuvs = +    accumArray (flip const) 0 (1,m)+      [ (c,v) | (c,_,v) <- cuvs ] ++  subtract_rec all@(last:_) = +    case subtract last of +      []  -> all+      new -> subtract_rec (new:all) ++  subtract [] = []+  subtract full@((c,u,v):rest) = +    if w >= v +      then (c,w,v) : subtract   rest+      else           subtract_b full+    where w = u - v+    +  subtract_b [] = []+  subtract_b ((c,u,v):rest) = +    if w /= 0 +      then (c,w,w) : subtract_b rest+      else           subtract_b rest+    where w = u - v++--------------------------------------------------------------------------------
Math/Combinat/Permutations.hs view
@@ -274,9 +274,9 @@     then Permutation result     else error "multiply: permutations of different sets"     where-	  (_,n) = bounds perm1-	  (_,m) = bounds perm2    -	  result = permute pi1 perm2    +    (_,n) = bounds perm1+    (_,m) = bounds perm2    +    result = permute pi1 perm2        infixr 7 `multiply`
Math/Combinat/Sign.hs view
@@ -1,48 +1,48 @@-
--- | Signs
-
-{-# LANGUAGE BangPatterns #-}
-module Math.Combinat.Sign where
-
---------------------------------------------------------------------------------
-
-import Data.Monoid
-
---------------------------------------------------------------------------------
-
-data Sign
-  = Plus
-  | Minus
-  deriving (Eq,Ord,Show,Read)
-
-instance Monoid Sign where
-  mempty  = Plus
-  mappend = mulSign
-  mconcat = productOfSigns
-
-signValue :: Num a => Sign -> a
-signValue s = case s of 
-  Plus  ->  1 
-  Minus -> -1 
-
-paritySign :: Integral a => a -> Sign
-paritySign x = if even x then Plus else Minus 
-
-oppositeSign :: Sign -> Sign
-oppositeSign s = case s of
-  Plus  -> Minus
-  Minus -> Plus
-
-mulSign :: Sign -> Sign -> Sign
-mulSign s1 s2 = case s1 of
-  Plus  -> s2
-  Minus -> oppositeSign s2
-
-productOfSigns :: [Sign] -> Sign
-productOfSigns = go Plus where
-  go !acc []     = acc
-  go !acc (x:xs) = case x of
-    Plus  -> go acc xs
-    Minus -> go (oppositeSign acc) xs
-
---------------------------------------------------------------------------------
++-- | Signs++{-# LANGUAGE BangPatterns #-}+module Math.Combinat.Sign where++--------------------------------------------------------------------------------++import Data.Monoid++--------------------------------------------------------------------------------++data Sign+  = Plus+  | Minus+  deriving (Eq,Ord,Show,Read)++instance Monoid Sign where+  mempty  = Plus+  mappend = mulSign+  mconcat = productOfSigns++signValue :: Num a => Sign -> a+signValue s = case s of +  Plus  ->  1 +  Minus -> -1 ++paritySign :: Integral a => a -> Sign+paritySign x = if even x then Plus else Minus ++oppositeSign :: Sign -> Sign+oppositeSign s = case s of+  Plus  -> Minus+  Minus -> Plus++mulSign :: Sign -> Sign -> Sign+mulSign s1 s2 = case s1 of+  Plus  -> s2+  Minus -> oppositeSign s2++productOfSigns :: [Sign] -> Sign+productOfSigns = go Plus where+  go !acc []     = acc+  go !acc (x:xs) = case x of+    Plus  -> go acc xs+    Minus -> go (oppositeSign acc) xs++--------------------------------------------------------------------------------
Math/Combinat/Tableaux.hs view
@@ -22,8 +22,11 @@ -- > ] -- +{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} module Math.Combinat.Tableaux where +--------------------------------------------------------------------------------+ import Data.List  import Math.Combinat.Helper@@ -31,6 +34,9 @@ import Math.Combinat.Partitions import Math.Combinat.ASCII +import Data.Map.Strict (Map)+import qualified Data.Map.Strict as Map+ -------------------------------------------------------------------------------- -- * Basic stuff @@ -41,6 +47,9 @@            $ (map . map) asciiShow            $ t +instance Show a => DrawASCII (Tableau a) where +  ascii = asciiTableau+ _shape :: Tableau a -> [Int] _shape t = map length t  @@ -93,6 +102,25 @@  columnWordToTableau :: Ord a => [a] -> Tableau a columnWordToTableau = transpose . rowWordToTableau++-- | Checks whether a sequence of positive integers is a /lattice word/, +-- which means that in every initial part of the sequence any number @i@+-- occurs at least as often as the number @i+1@+--+isLatticeWord :: [Int] -> Bool+isLatticeWord = go Map.empty where+  go :: Map Int Int -> [Int] -> Bool+  go _      []     = True+  go !table (i:is) =+    if check i+      then go table' is+      else False+    where+      table'  = Map.insertWith (+) i 1 table+      check j = j==1 || cnt (j-1) >= cnt j+      cnt j   = case Map.lookup j table' of+        Just k  -> k+        Nothing -> 0      -------------------------------------------------------------------------------- -- * Standard Young tableaux
Math/Combinat/Tableaux/GelfandTsetlin.hs view
@@ -1,343 +1,341 @@-
--- | Gelfand-Tsetlin patterns and Kostka numbers.
---
--- Gelfand-Tsetlin patterns (or tableaux) are triangular arrays like
---
--- > [ 3 ]
--- > [ 3 , 2 ]
--- > [ 3 , 1 , 0 ]
--- > [ 2 , 0 , 0 , 0 ]
---
--- with both rows and columns non-increasing non-negative integers.
--- Note: these are in bijection with the semi-standard Young tableaux.
---
--- If we add the further restriction that
--- the top diagonal reads @lambda@, 
--- and the diagonal sums are partial sums of @mu@, where @lambda@ and @mu@ are two
--- partitions (in this case @lambda=[3,2]@ and @mu=[2,1,1,1]@), 
--- then the number of the resulting patterns 
--- or tableaux is the Kostka number @K(lambda,mu)@.
--- Actually @mu@ doesn't even need to the be non-increasing.
---
-
-{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}
-module Math.Combinat.Tableaux.GelfandTsetlin where
-
---------------------------------------------------------------------------------
-
-import Data.List
-import Data.Maybe
-import Data.Monoid
-import Data.Ord
-
-import Control.Monad
-import Control.Monad.Trans.State
-
-import Data.Map (Map)
-import qualified Data.Map as Map
-
-import Math.Combinat.Partitions.Integer
-import Math.Combinat.Tableaux
-import Math.Combinat.Helper
-import Math.Combinat.ASCII
-
---------------------------------------------------------------------------------
--- * Kostka numbers
-
--- | Kostka numbers (via counting Gelfand-Tsetlin patterns). See for example <http://en.wikipedia.org/wiki/Kostka_number>
---
--- @K(lambda,mu)==0@ unless @lambda@ dominates @mu@:
---
--- > [ mu | mu <- partitions (weight lam) , kostkaNumber lam mu > 0 ] == dominatedPartitions lam
---
-kostkaNumber :: Partition -> Partition -> Int
-kostkaNumber = countKostkaGelfandTsetlinPatterns
-
--- | Very naive (and slow) implementation of Kostka numbers, for reference.
-kostkaNumberReferenceNaive :: Partition -> Partition -> Int
-kostkaNumberReferenceNaive plambda pmu@(Partition mu) = length stuff where
-  stuff = [ 1 | t <- semiStandardYoungTableaux k plambda , cond t ]
-  k = length mu
-  cond t = [ (head xs, length xs) | xs <- group (sort $ concat t) ] == zip [1..] mu 
-
---------------------------------------------------------------------------------
-
--- | Lists all (positive) Kostka numbers @K(lambda,mu)@ with the given @lambda@:
---
--- > kostkaNumbersWithGivenLambda lambda == Map.fromList [ (mu , kostkaNumber lambda mu) | mu <- dominatedPartitions lambda ]
---
--- It's much faster than computing the individual Kostka numbers, but not as fast
--- as it could be.
---
-{-# SPECIALIZE kostkaNumbersWithGivenLambda :: Partition -> Map Partition Int     #-}
-{-# SPECIALIZE kostkaNumbersWithGivenLambda :: Partition -> Map Partition Integer #-}
-kostkaNumbersWithGivenLambda :: forall coeff. Num coeff => Partition -> Map Partition coeff
-kostkaNumbersWithGivenLambda plambda@(Partition lam) = evalState (worker lam) Map.empty where
-
-  worker :: [Int] -> State (Map Partition (Map Partition coeff)) (Map Partition coeff)
-  worker unlam = case unlam of
-    [] -> return $ Map.singleton (Partition []) 1
-    _  -> do
-      cache <- get
-      case Map.lookup (Partition unlam) cache of
-        Just sol -> return sol
-        Nothing  -> do
-          let s = foldl' (+) 0 unlam
-          subsols <- forM (prevLambdas0 unlam) $ \p -> do
-            sub <- worker p 
-            let t = s - foldl' (+) 0 p              
-                f (Partition xs , c) = case xs of
-                  (y:_) -> if t >= y then Just (Partition (t:xs) , c) else Nothing
-                  []    -> if t >  0 then Just (Partition [t]    , c) else Nothing
-            if t > 0
-              then return $ Map.fromList $ mapMaybe f $ Map.toList sub
-              else return $ Map.empty
-
-          let sol = Map.unionsWith (+) subsols
-          put $! (Map.insert (Partition unlam) sol cache)
-          return sol
-
-  -- needs decreasing sequence
-  prevLambdas0 :: [Int] -> [[Int]]
-  prevLambdas0 (l:ls) = go l ls where
-    go b [a]    = [ [x]   | x <- [a..b] ] ++ [ [x,y] | x <- [a..b] , y<-[1..a] ]
-    go b (a:as) = [ x:xs  | x <- [a..b] , xs <- go a as ]
-    go b []     = [] : [ [j] | j <- [1..b] ]
-  prevLambdas0 []  = []
-
--- | Lists all (positive) Kostka numbers @K(lambda,mu)@ with the given @mu@:
---
--- > kostkaNumbersWithGivenMu mu == Map.fromList [ (lambda , kostkaNumber lambda mu) | lambda <- dominatingPartitions mu ]
---
--- This function uses the iterated Pieri rule, and is relatively fast.
---
-kostkaNumbersWithGivenMu :: Partition -> Map Partition Int
-kostkaNumbersWithGivenMu (Partition mu) = iteratedPieriRule (reverse mu)
-
---------------------------------------------------------------------------------
--- * Gelfand-Tsetlin patterns
-
--- | A Gelfand-Tstetlin tableau
-type GT = [[Int]]
-
-
-asciiGT :: GT -> ASCII
-asciiGT gt = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) 
-           $ (map . map) asciiShow
-           $ gt
-
-
-kostkaGelfandTsetlinPatterns :: Partition -> Partition -> [GT]
-kostkaGelfandTsetlinPatterns lambda (Partition mu) = kostkaGelfandTsetlinPatterns' lambda mu
-
--- | Generates all Kostka-Gelfand-Tsetlin tableau, that is, triangular arrays like
---
--- > [ 3 ]
--- > [ 3 , 2 ]
--- > [ 3 , 1 , 0 ]
--- > [ 2 , 0 , 0 , 0 ]
---
--- with both rows and column non-increasing such that
--- the top diagonal read lambda (in this case @lambda=[3,2]@) and the diagonal sums
--- are partial sums of mu (in this case @mu=[2,1,1,1]@)
---
--- The number of such GT tableaux is the Kostka
--- number K(lambda,mu).
---
-kostkaGelfandTsetlinPatterns' :: Partition -> [Int] -> [GT]
-kostkaGelfandTsetlinPatterns' plam@(Partition lambda0) mu0
-  | minimum mu0 < 0                       = []
-  | wlam == 0                             = if wmu == 0 then [ [] ] else []
-  | wmu  == wlam && plam `dominates` pmu  = list
-  | otherwise                             = []
-  where
-
-    pmu = mkPartition mu0
-
-    nlam = length lambda0
-    nmu  = length mu0
-
-    n = max nlam nmu
-
-    lambda = lambda0 ++ replicate (n - nlam) 0
-    mu     = mu0     ++ replicate (n - nmu ) 0
-
-    revlam = reverse lambda
-
-    wmu  = sum' mu
-    wlam = sum' lambda
-
-    list = worker 
-             revlam 
-             (scanl1 (+) mu) 
-             (replicate (n-1) 0) 
-             (replicate (n  ) 0) 
-             []
-
-    worker
-      :: [Int]       -- lambda_i in reverse order
-      -> [Int]       -- partial sums of mu
-      -> [Int]       -- sums of the tails of previous rows
-      -> [Int]       -- last row
-      -> [[Int]]     -- the lower part of GT tableau we accumulated so far (this is not needed if we only want to count)
-      -> [GT]   
-
-    worker (rl:rls) (smu:smus) (a:acc) (lastx0:lastrowt) table = stuff 
-      where
-        x0 = smu - a
-        stuff = concat 
-          [ worker rls smus (zipWith (+) acc (tail row)) (init row) (row:table)
-          | row <- boundedNonIncrSeqs' x0 (map (max rl) (max lastx0 x0 : lastrowt)) lambda
-          ]
-    worker [rl] _ _ _ table = [ [rl]:table ] 
-    worker []   _ _ _ _     = [ []         ]
-
-    boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]
-    boundedNonIncrSeqs' = go where
-      go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ]
-      go _  []     _      = [[]]
-      go _  _      []     = [[]]
-
---------------------------------------------------------------------------------
-
--- | This returns the corresponding Kostka number:
---
--- > countKostkaGelfandTsetlinPatterns lambda mu == length (kostkaGelfandTsetlinPatterns lambda mu)
--- 
-countKostkaGelfandTsetlinPatterns :: Partition -> Partition -> Int
-countKostkaGelfandTsetlinPatterns plam@(Partition lambda0) pmu@(Partition mu0) 
-  | wlam == 0                             = if wmu == 0 then 1 else 0
-  | wmu  == wlam && plam `dominates` pmu  = cnt
-  | otherwise                             = 0
-  where
-
-    nlam = length lambda0
-    nmu  = length mu0
-
-    n = max nlam nmu
-
-    lambda = lambda0 ++ replicate (n - nlam) 0
-    mu     = mu0     ++ replicate (n - nmu ) 0
-
-    revlam = reverse lambda
-
-    wmu  = sum' mu
-    wlam = sum' lambda
-
-    cnt = worker 
-            revlam 
-            (scanl1 (+) mu) 
-            (replicate (n-1) 0) 
-            (replicate (n  ) 0) 
-
-    worker
-      :: [Int]       -- lambda_i in reverse order
-      -> [Int]       -- partial sums of mu
-      -> [Int]       -- sums of the tails of previous rows
-      -> [Int]       -- last row
-      -> Int
-
-    worker (rl:rls) (smu:smus) (a:acc) (lastx0:lastrowt) = stuff 
-      where
-        x0 = smu - a
-        stuff = sum'
-          [ worker rls smus (zipWith (+) acc (tail row)) (init row) 
-          | row <- boundedNonIncrSeqs' x0 (map (max rl) (max lastx0 x0 : lastrowt)) lambda
-          ]
-    worker [rl] _ _ _ = 1 
-    worker []   _ _ _ = 1
-
-    boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]
-    boundedNonIncrSeqs' = go where
-      go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ]
-      go _  []     _      = [[]]
-      go _  _      []     = [[]]
-
---------------------------------------------------------------------------------
-
-{-
-
--- | All non-increasing sentences between a lower and an upper bound
-boundedNonIncrSeqs :: [Int] -> [Int] -> [[Int]]
-boundedNonIncrSeqs as bs = case bs of  
-  (h0:_) -> boundedNonIncrSeqs' h0 as bs
-  []     -> [[]]
-
--- | All non-increasing sentences between a lower and an upper bound, and also less-or-equal than the given number
-boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]
-boundedNonIncrSeqs' = go where
-  go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ]
-  go _  []     _      = [[]]
-  go _  _      []     = [[]]
-
--- | All non-decreasing sentences between a lower and an upper bound
-boundedNonDecrSeqs :: [Int] -> [Int] -> [[Int]]
-boundedNonDecrSeqs = boundedNonDecrSeqs' 0
-
--- | All non-decreasing sentences between a lower and an upper bound, and also greator-or-equal then the given number
-boundedNonDecrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]
-boundedNonDecrSeqs' h0 = go (max 0 h0) where
-  go h0 (a:as) (b:bs) = [ h:hs | h <- [(max h0 a)..b] , hs <- go h as bs ]
-  go _  []     _      = [[]]
-  go _  _      []     = [[]]
-
--}
-
---------------------------------------------------------------------------------
--- * The iterated Pieri rule 
-
--- | Computes the Schur expansion of @h[n1]*h[n2]*h[n3]*...*h[nk]@ via iterating the Pieri rule.
--- Note: the coefficients are actually the Kostka numbers; the following is true:
---
--- > Map.toList (iteratedPieriRule (fromPartition mu))  ==  [ (lam, kostkaNumber lam mu) | lam <- dominatingPartitions mu ]
--- 
--- This should be faster than individually computing all these Kostka numbers.
---
-iteratedPieriRule :: Num coeff => [Int] -> Map Partition coeff
-iteratedPieriRule = iteratedPieriRule' (Partition [])
-
--- | Iterating the Pieri rule, we can compute the Schur expansion of
--- @h[lambda]*h[n1]*h[n2]*h[n3]*...*h[nk]@
-iteratedPieriRule' :: Num coeff => Partition -> [Int] -> Map Partition coeff
-iteratedPieriRule' plambda ns = iteratedPieriRule'' (plambda,1) ns
-
-{-# SPECIALIZE iteratedPieriRule'' :: (Partition,Int    ) -> [Int] -> Map Partition Int     #-}
-{-# SPECIALIZE iteratedPieriRule'' :: (Partition,Integer) -> [Int] -> Map Partition Integer #-}
-iteratedPieriRule'' :: Num coeff => (Partition,coeff) -> [Int] -> Map Partition coeff
-iteratedPieriRule'' (plambda,coeff0) ns = worker (Map.singleton plambda coeff0) ns where
-  worker old []     = old
-  worker old (n:ns) = worker new ns where
-    stuff = [ (coeff, pieriRule lam n) | (lam,coeff) <- Map.toList old ] 
-    new   = foldl' f Map.empty stuff 
-    f t0 (c,ps) = foldl' (\t p -> Map.insertWith (+) p c t) t0 ps  
-
---------------------------------------------------------------------------------
-
--- | Computes the Schur expansion of @e[n1]*e[n2]*e[n3]*...*e[nk]@ via iterating the Pieri rule.
--- Note: the coefficients are actually the Kostka numbers; the following is true:
---
--- > Map.toList (iteratedDualPieriRule (fromPartition mu))  ==  
--- >   [ (dualPartition lam, kostkaNumber lam mu) | lam <- dominatingPartitions mu ]
--- 
--- This should be faster than individually computing all these Kostka numbers.
--- It is a tiny bit slower than 'iteratedPieriRule'.
---
-iteratedDualPieriRule :: Num coeff => [Int] -> Map Partition coeff
-iteratedDualPieriRule = iteratedDualPieriRule' (Partition [])
-
--- | Iterating the Pieri rule, we can compute the Schur expansion of
--- @e[lambda]*e[n1]*e[n2]*e[n3]*...*e[nk]@
-iteratedDualPieriRule' :: Num coeff => Partition -> [Int] -> Map Partition coeff
-iteratedDualPieriRule' plambda ns = iteratedDualPieriRule'' (plambda,1) ns
-
-{-# SPECIALIZE iteratedDualPieriRule'' :: (Partition,Int    ) -> [Int] -> Map Partition Int     #-}
-{-# SPECIALIZE iteratedDualPieriRule'' :: (Partition,Integer) -> [Int] -> Map Partition Integer #-}
-iteratedDualPieriRule'' :: Num coeff => (Partition,coeff) -> [Int] -> Map Partition coeff
-iteratedDualPieriRule'' (plambda,coeff0) ns = worker (Map.singleton plambda coeff0) ns where
-  worker old []     = old
-  worker old (n:ns) = worker new ns where
-    stuff = [ (coeff, dualPieriRule lam n) | (lam,coeff) <- Map.toList old ] 
-    new   = foldl' f Map.empty stuff 
-    f t0 (c,ps) = foldl' (\t p -> Map.insertWith (+) p c t) t0 ps  
-
---------------------------------------------------------------------------------
++-- | Gelfand-Tsetlin patterns and Kostka numbers.+--+-- Gelfand-Tsetlin patterns (or tableaux) are triangular arrays like+--+-- > [ 3 ]+-- > [ 3 , 2 ]+-- > [ 3 , 1 , 0 ]+-- > [ 2 , 0 , 0 , 0 ]+--+-- with both rows and columns non-increasing non-negative integers.+-- Note: these are in bijection with the semi-standard Young tableaux.+--+-- If we add the further restriction that+-- the top diagonal reads @lambda@, +-- and the diagonal sums are partial sums of @mu@, where @lambda@ and @mu@ are two+-- partitions (in this case @lambda=[3,2]@ and @mu=[2,1,1,1]@), +-- then the number of the resulting patterns +-- or tableaux is the Kostka number @K(lambda,mu)@.+-- Actually @mu@ doesn't even need to the be non-increasing.+--++{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}+module Math.Combinat.Tableaux.GelfandTsetlin where++--------------------------------------------------------------------------------++import Data.List+import Data.Maybe+import Data.Monoid+import Data.Ord++import Control.Monad+import Control.Monad.Trans.State++import Data.Map (Map)+import qualified Data.Map as Map++import Math.Combinat.Partitions.Integer+import Math.Combinat.Tableaux+import Math.Combinat.Helper+import Math.Combinat.ASCII++--------------------------------------------------------------------------------+-- * Kostka numbers++-- | Kostka numbers (via counting Gelfand-Tsetlin patterns). See for example <http://en.wikipedia.org/wiki/Kostka_number>+--+-- @K(lambda,mu)==0@ unless @lambda@ dominates @mu@:+--+-- > [ mu | mu <- partitions (weight lam) , kostkaNumber lam mu > 0 ] == dominatedPartitions lam+--+kostkaNumber :: Partition -> Partition -> Int+kostkaNumber = countKostkaGelfandTsetlinPatterns++-- | Very naive (and slow) implementation of Kostka numbers, for reference.+kostkaNumberReferenceNaive :: Partition -> Partition -> Int+kostkaNumberReferenceNaive plambda pmu@(Partition mu) = length stuff where+  stuff  = [ (1::Int) | t <- semiStandardYoungTableaux k plambda , cond t ]+  k      = length mu+  cond t = [ (head xs, length xs) | xs <- group (sort $ concat t) ] == zip [1..] mu ++--------------------------------------------------------------------------------++-- | Lists all (positive) Kostka numbers @K(lambda,mu)@ with the given @lambda@:+--+-- > kostkaNumbersWithGivenLambda lambda == Map.fromList [ (mu , kostkaNumber lambda mu) | mu <- dominatedPartitions lambda ]+--+-- It's much faster than computing the individual Kostka numbers, but not as fast+-- as it could be.+--+{-# SPECIALIZE kostkaNumbersWithGivenLambda :: Partition -> Map Partition Int     #-}+{-# SPECIALIZE kostkaNumbersWithGivenLambda :: Partition -> Map Partition Integer #-}+kostkaNumbersWithGivenLambda :: forall coeff. Num coeff => Partition -> Map Partition coeff+kostkaNumbersWithGivenLambda plambda@(Partition lam) = evalState (worker lam) Map.empty where++  worker :: [Int] -> State (Map Partition (Map Partition coeff)) (Map Partition coeff)+  worker unlam = case unlam of+    [] -> return $ Map.singleton (Partition []) 1+    _  -> do+      cache <- get+      case Map.lookup (Partition unlam) cache of+        Just sol -> return sol+        Nothing  -> do+          let s = foldl' (+) 0 unlam+          subsols <- forM (prevLambdas0 unlam) $ \p -> do+            sub <- worker p +            let t = s - foldl' (+) 0 p              +                f (Partition xs , c) = case xs of+                  (y:_) -> if t >= y then Just (Partition (t:xs) , c) else Nothing+                  []    -> if t >  0 then Just (Partition [t]    , c) else Nothing+            if t > 0+              then return $ Map.fromList $ mapMaybe f $ Map.toList sub+              else return $ Map.empty++          let sol = Map.unionsWith (+) subsols+          put $! (Map.insert (Partition unlam) sol cache)+          return sol++  -- needs decreasing sequence+  prevLambdas0 :: [Int] -> [[Int]]+  prevLambdas0 (l:ls) = go l ls where+    go b [a]    = [ [x]   | x <- [a..b] ] ++ [ [x,y] | x <- [a..b] , y<-[1..a] ]+    go b (a:as) = [ x:xs  | x <- [a..b] , xs <- go a as ]+    go b []     = [] : [ [j] | j <- [1..b] ]+  prevLambdas0 []  = []++-- | Lists all (positive) Kostka numbers @K(lambda,mu)@ with the given @mu@:+--+-- > kostkaNumbersWithGivenMu mu == Map.fromList [ (lambda , kostkaNumber lambda mu) | lambda <- dominatingPartitions mu ]+--+-- This function uses the iterated Pieri rule, and is relatively fast.+--+kostkaNumbersWithGivenMu :: Partition -> Map Partition Int+kostkaNumbersWithGivenMu (Partition mu) = iteratedPieriRule (reverse mu)++--------------------------------------------------------------------------------+-- * Gelfand-Tsetlin patterns++-- | A Gelfand-Tstetlin tableau+type GT = [[Int]]++asciiGT :: GT -> ASCII+asciiGT gt = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) +           $ (map . map) asciiShow+           $ gt++kostkaGelfandTsetlinPatterns :: Partition -> Partition -> [GT]+kostkaGelfandTsetlinPatterns lambda (Partition mu) = kostkaGelfandTsetlinPatterns' lambda mu++-- | Generates all Kostka-Gelfand-Tsetlin tableau, that is, triangular arrays like+--+-- > [ 3 ]+-- > [ 3 , 2 ]+-- > [ 3 , 1 , 0 ]+-- > [ 2 , 0 , 0 , 0 ]+--+-- with both rows and column non-increasing such that+-- the top diagonal read lambda (in this case @lambda=[3,2]@) and the diagonal sums+-- are partial sums of mu (in this case @mu=[2,1,1,1]@)+--+-- The number of such GT tableaux is the Kostka+-- number K(lambda,mu).+--+kostkaGelfandTsetlinPatterns' :: Partition -> [Int] -> [GT]+kostkaGelfandTsetlinPatterns' plam@(Partition lambda0) mu0+  | minimum mu0 < 0                       = []+  | wlam == 0                             = if wmu == 0 then [ [] ] else []+  | wmu  == wlam && plam `dominates` pmu  = list+  | otherwise                             = []+  where++    pmu = mkPartition mu0++    nlam = length lambda0+    nmu  = length mu0++    n = max nlam nmu++    lambda = lambda0 ++ replicate (n - nlam) 0+    mu     = mu0     ++ replicate (n - nmu ) 0++    revlam = reverse lambda++    wmu  = sum' mu+    wlam = sum' lambda++    list = worker +             revlam +             (scanl1 (+) mu) +             (replicate (n-1) 0) +             (replicate (n  ) 0) +             []++    worker+      :: [Int]       -- lambda_i in reverse order+      -> [Int]       -- partial sums of mu+      -> [Int]       -- sums of the tails of previous rows+      -> [Int]       -- last row+      -> [[Int]]     -- the lower part of GT tableau we accumulated so far (this is not needed if we only want to count)+      -> [GT]   ++    worker (rl:rls) (smu:smus) (a:acc) (lastx0:lastrowt) table = stuff +      where+        x0 = smu - a+        stuff = concat +          [ worker rls smus (zipWith (+) acc (tail row)) (init row) (row:table)+          | row <- boundedNonIncrSeqs' x0 (map (max rl) (max lastx0 x0 : lastrowt)) lambda+          ]+    worker [rl] _ _ _ table = [ [rl]:table ] +    worker []   _ _ _ _     = [ []         ]++    boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]+    boundedNonIncrSeqs' = go where+      go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ]+      go _  []     _      = [[]]+      go _  _      []     = [[]]++--------------------------------------------------------------------------------++-- | This returns the corresponding Kostka number:+--+-- > countKostkaGelfandTsetlinPatterns lambda mu == length (kostkaGelfandTsetlinPatterns lambda mu)+-- +countKostkaGelfandTsetlinPatterns :: Partition -> Partition -> Int+countKostkaGelfandTsetlinPatterns plam@(Partition lambda0) pmu@(Partition mu0) +  | wlam == 0                             = if wmu == 0 then 1 else 0+  | wmu  == wlam && plam `dominates` pmu  = cnt+  | otherwise                             = 0+  where++    nlam = length lambda0+    nmu  = length mu0++    n = max nlam nmu++    lambda = lambda0 ++ replicate (n - nlam) 0+    mu     = mu0     ++ replicate (n - nmu ) 0++    revlam = reverse lambda++    wmu  = sum' mu+    wlam = sum' lambda++    cnt = worker +            revlam +            (scanl1 (+) mu) +            (replicate (n-1) 0) +            (replicate (n  ) 0) ++    worker+      :: [Int]       -- lambda_i in reverse order+      -> [Int]       -- partial sums of mu+      -> [Int]       -- sums of the tails of previous rows+      -> [Int]       -- last row+      -> Int++    worker (rl:rls) (smu:smus) (a:acc) (lastx0:lastrowt) = stuff +      where+        x0 = smu - a+        stuff = sum'+          [ worker rls smus (zipWith (+) acc (tail row)) (init row) +          | row <- boundedNonIncrSeqs' x0 (map (max rl) (max lastx0 x0 : lastrowt)) lambda+          ]+    worker [rl] _ _ _ = 1 +    worker []   _ _ _ = 1++    boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]+    boundedNonIncrSeqs' = go where+      go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ]+      go _  []     _      = [[]]+      go _  _      []     = [[]]++--------------------------------------------------------------------------------++{-++-- | All non-increasing sentences between a lower and an upper bound+boundedNonIncrSeqs :: [Int] -> [Int] -> [[Int]]+boundedNonIncrSeqs as bs = case bs of  +  (h0:_) -> boundedNonIncrSeqs' h0 as bs+  []     -> [[]]++-- | All non-increasing sentences between a lower and an upper bound, and also less-or-equal than the given number+boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]+boundedNonIncrSeqs' = go where+  go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ]+  go _  []     _      = [[]]+  go _  _      []     = [[]]++-- | All non-decreasing sentences between a lower and an upper bound+boundedNonDecrSeqs :: [Int] -> [Int] -> [[Int]]+boundedNonDecrSeqs = boundedNonDecrSeqs' 0++-- | All non-decreasing sentences between a lower and an upper bound, and also greator-or-equal then the given number+boundedNonDecrSeqs' :: Int -> [Int] -> [Int] -> [[Int]]+boundedNonDecrSeqs' h0 = go (max 0 h0) where+  go h0 (a:as) (b:bs) = [ h:hs | h <- [(max h0 a)..b] , hs <- go h as bs ]+  go _  []     _      = [[]]+  go _  _      []     = [[]]++-}++--------------------------------------------------------------------------------+-- * The iterated Pieri rule ++-- | Computes the Schur expansion of @h[n1]*h[n2]*h[n3]*...*h[nk]@ via iterating the Pieri rule.+-- Note: the coefficients are actually the Kostka numbers; the following is true:+--+-- > Map.toList (iteratedPieriRule (fromPartition mu))  ==  [ (lam, kostkaNumber lam mu) | lam <- dominatingPartitions mu ]+-- +-- This should be faster than individually computing all these Kostka numbers.+--+iteratedPieriRule :: Num coeff => [Int] -> Map Partition coeff+iteratedPieriRule = iteratedPieriRule' (Partition [])++-- | Iterating the Pieri rule, we can compute the Schur expansion of+-- @h[lambda]*h[n1]*h[n2]*h[n3]*...*h[nk]@+iteratedPieriRule' :: Num coeff => Partition -> [Int] -> Map Partition coeff+iteratedPieriRule' plambda ns = iteratedPieriRule'' (plambda,1) ns++{-# SPECIALIZE iteratedPieriRule'' :: (Partition,Int    ) -> [Int] -> Map Partition Int     #-}+{-# SPECIALIZE iteratedPieriRule'' :: (Partition,Integer) -> [Int] -> Map Partition Integer #-}+iteratedPieriRule'' :: Num coeff => (Partition,coeff) -> [Int] -> Map Partition coeff+iteratedPieriRule'' (plambda,coeff0) ns = worker (Map.singleton plambda coeff0) ns where+  worker old []     = old+  worker old (n:ns) = worker new ns where+    stuff = [ (coeff, pieriRule lam n) | (lam,coeff) <- Map.toList old ] +    new   = foldl' f Map.empty stuff +    f t0 (c,ps) = foldl' (\t p -> Map.insertWith (+) p c t) t0 ps  ++--------------------------------------------------------------------------------++-- | Computes the Schur expansion of @e[n1]*e[n2]*e[n3]*...*e[nk]@ via iterating the Pieri rule.+-- Note: the coefficients are actually the Kostka numbers; the following is true:+--+-- > Map.toList (iteratedDualPieriRule (fromPartition mu))  ==  +-- >   [ (dualPartition lam, kostkaNumber lam mu) | lam <- dominatingPartitions mu ]+-- +-- This should be faster than individually computing all these Kostka numbers.+-- It is a tiny bit slower than 'iteratedPieriRule'.+--+iteratedDualPieriRule :: Num coeff => [Int] -> Map Partition coeff+iteratedDualPieriRule = iteratedDualPieriRule' (Partition [])++-- | Iterating the Pieri rule, we can compute the Schur expansion of+-- @e[lambda]*e[n1]*e[n2]*e[n3]*...*e[nk]@+iteratedDualPieriRule' :: Num coeff => Partition -> [Int] -> Map Partition coeff+iteratedDualPieriRule' plambda ns = iteratedDualPieriRule'' (plambda,1) ns++{-# SPECIALIZE iteratedDualPieriRule'' :: (Partition,Int    ) -> [Int] -> Map Partition Int     #-}+{-# SPECIALIZE iteratedDualPieriRule'' :: (Partition,Integer) -> [Int] -> Map Partition Integer #-}+iteratedDualPieriRule'' :: Num coeff => (Partition,coeff) -> [Int] -> Map Partition coeff+iteratedDualPieriRule'' (plambda,coeff0) ns = worker (Map.singleton plambda coeff0) ns where+  worker old []     = old+  worker old (n:ns) = worker new ns where+    stuff = [ (coeff, dualPieriRule lam n) | (lam,coeff) <- Map.toList old ] +    new   = foldl' f Map.empty stuff +    f t0 (c,ps) = foldl' (\t p -> Map.insertWith (+) p c t) t0 ps  ++--------------------------------------------------------------------------------
Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs view
@@ -47,6 +47,7 @@ -- to the dimension), which encode the combinatorics of Kostka numbers. -- +{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} module Math.Combinat.Tableaux.GelfandTsetlin.Cone   (      -- * Types@@ -115,12 +116,10 @@   range     (a,b) = map deIndex' [ index' a .. index' b ]    rangeSize (a,b) = index' b - index' a + 1  -{-# SPECIALIZE triangularArrayUnsafe :: Tableau Int -> TriangularArray Int #-} triangularArrayUnsafe :: Tableau a -> TriangularArray a triangularArrayUnsafe tableau = listArray (Tri (1,1),Tri (k,k)) (concat tableau)    where k = length tableau -{-# SPECIALIZE fromTriangularArray :: TriangularArray Int -> Tableau Int #-} fromTriangularArray :: TriangularArray a -> Tableau a fromTriangularArray arr = (map.map) snd $ groupBy (equating f) $ assocs arr   where f = fst . unTri . fst@@ -135,6 +134,12 @@                  $ (map . map) asciiShow                  $ xxs +instance Show a => DrawASCII (TriangularArray a) where+  ascii = asciiTriangularArray++-- instance Show a => DrawASCII (Tableau a) where+--   ascii = asciiTableau+ --------------------------------------------------------------------------------  -- "fractional fillings"@@ -149,7 +154,6 @@ nextHole :: Hole -> Hole nextHole (Hole k l) = Hole k (l+1) -{-# SPECIALIZE reverseTableau :: [[Int]] -> [[Int]] #-} reverseTableau :: [[a]] -> [[a]] reverseTableau = reverse . map reverse
Math/Combinat/Tableaux/LittlewoodRichardson.hs view
@@ -2,22 +2,38 @@ -- | The Littlewood-Richardson rule  module Math.Combinat.Tableaux.LittlewoodRichardson -  ( lrRule , _lrRule+  ( lrRule , _lrRule +  , lrRuleNaive   )    where  --------------------------------------------------------------------------------  import Data.List+import Data.Maybe  import Math.Combinat.Partitions.Integer import Math.Combinat.Partitions.Skew+import Math.Combinat.Tableaux+import Math.Combinat.Tableaux.Skew  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"+--+lrRuleNaive :: SkewPartition -> Map Partition Int+lrRuleNaive skew = final where+  n     = skewPartitionWeight skew+  ssst  = semiStandardSkewTableaux n skew +  final = foldl' f Map.empty $ catMaybes [ skewTableauRowContent skew | skew <- ssst  ]+  f old nu = Map.insertWith (+) nu 1 old++--------------------------------------------------------------------------------+ -- | @lrRule@ computes the expansion of a skew Schur function  -- @s[lambda/mu]@ via the Littlewood-Richardson rule. --@@ -35,7 +51,7 @@ {-# SPECIALIZE _lrRule :: Partition -> Partition -> Map Partition Integer #-} _lrRule :: Num coeff => Partition -> Partition -> Map Partition coeff _lrRule plam@(Partition lam) pmu@(Partition mu0) = -  if not (plam `dominates` pmu) +  if not (pmu `isSubPartitionOf` plam)      then Map.empty     else foldl' f Map.empty [ nu | (nu,_) <- fillings n diagram ]   where@@ -107,8 +123,19 @@   ub = if upper>0      then min (length shape) (lpart !! (upper-1))       else      length shape+   nlist = filter (>0) $ map f [lb+1..ub]    f j   = if j==1 || shape!!(j-2) > shape!!(j-1) then j else 0++{-+  -- another nlist implementation, but doesn't seem to be faster+  (h0:hs0) = drop lb (-666:shape)+  nlist = go h0 hs0 [lb+1..ub] where+    go !lasth (h:hs) (j:js) = if j==1 || lasth > h +      then j : go h hs js +      else     go h hs js+    go _      _      []     = []+-}    -- increments the i-th element (starting from 1)   incr :: Int -> [Int] -> [Int]
Math/Combinat/Tableaux/Skew.hs view
@@ -1,78 +1,180 @@-
--- | Skew tableaux are skew partitions filled with numbers.
-
-{-# LANGUAGE BangPatterns #-}
-
-module Math.Combinat.Tableaux.Skew where
-
---------------------------------------------------------------------------------
-
-import Data.List
-
-import Math.Combinat.Partitions.Integer
-import Math.Combinat.Partitions.Skew
-import Math.Combinat.Tableaux
-import Math.Combinat.ASCII
-
---------------------------------------------------------------------------------
-
--- | A skew tableau is represented by a list of offsets and entries
-newtype SkewTableau a = SkewTableau [(Int,[a])] deriving (Eq,Ord,Show)
-
-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 ]
-
--- | Semi-standard skew tableaux filled with numbers @[1..n]@
-semiStandardSkewTableaux :: Int -> SkewPartition -> [SkewTableau Int]
-semiStandardSkewTableaux n (SkewPartition abs) = map SkewTableau stuff where
-
-  stuff = worker as bs ds (repeat 1) 
-  (as,bs) = unzip abs
-  ds = diffSequence as
-  
-  -- | @worker inner outerMinusInner innerdiffs lowerbound
-  worker :: [Int] -> [Int] -> [Int] -> [Int] -> [[(Int,[Int])]]
-  worker (a:as) (b:bs) (d:ds) lb = [ (a,this):rest 
-                                   | this <- row b 1 lb 
-                                   , let lb' = (replicate d 1 ++ map (+1) this) 
-                                   , rest <- worker as bs ds lb' ] 
-  worker []     _      _      _  = [ [] ]
-
-  -- @row length minimum lowerbound@
-  row 0  _  _       = [[]]
-  row _  _  []      = []
-  row !k !m (!a:as) = [ x:xs | x <- [(max a m)..n] , xs <- row (k-1) x as ] 
-
-{-
--- | from a sequence @[a1,a2,..,an]@ computes the sequence of differences
--- @[a1-a2,a2-a3,...,an-0]@
-diffSequence :: [Int] -> [Int]
-diffSequence = go where
-  go (x:ys@(y:_)) = (x-y) : go ys 
-  go [x] = [x]
-  go []  = []
--}
-
---------------------------------------------------------------------------------
-
-asciiSkewTableau :: Show a => SkewTableau a -> ASCII
-asciiSkewTableau = asciiSkewTableau' "." EnglishNotation
-
-asciiSkewTableau' 
-  :: Show a
-  => 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
-  stuff = case orient of
-    EnglishNotation    -> es
-    EnglishNotationCCW -> reverse (transpose es)
-    FrenchNotation     -> reverse es
-  inner = asciiFromString innerstr
-  es = [ replicate a inner ++ map asciiShow xs | (a,xs) <- axs ]
-
---------------------------------------------------------------------------------
++-- | Skew tableaux are skew partitions filled with numbers.++{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}++module Math.Combinat.Tableaux.Skew where++--------------------------------------------------------------------------------++import Data.List++import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Skew+import Math.Combinat.Tableaux+import Math.Combinat.ASCII++import Data.Map.Strict (Map)+import qualified Data.Map.Strict as Map++--------------------------------------------------------------------------------++-- | A skew tableau is represented by a list of offsets and entries+newtype SkewTableau a = SkewTableau [(Int,[a])] deriving (Eq,Ord,Show)++-- unSkewTableau :: SkewTableau a -> [(Int,[a])]+-- unSkewTableau (SkewTableau a) = a++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 ]++--------------------------------------------------------------------------------++dualSkewTableau :: forall a. SkewTableau a -> SkewTableau a+dualSkewTableau (SkewTableau axs) = SkewTableau (go axs) where++  go []  = []  +  go axs = case sub 0 axs of+    (0,[]) -> []+    this   -> this : go (strip axs)++  strip :: [(Int,[a])] -> [(Int,[a])]+  strip []            = []+  strip ((a,xs):rest) = if a>0 +    then (a-1,xs) : strip rest+    else case xs of+      []     -> []+      (z:zs) -> case zs of+        []      -> []+        _       -> (0,zs) : strip rest++  sub :: Int -> [(Int,[a])] -> (Int,[a])+  sub !b [] = (b,[])+  sub !b ((a,this):rest) = if a>0 +    then sub (b+1) rest  +    else (b,ys) where      +      ys = map head $ takeWhile (not . null) (this : map snd rest)++{-+test_dualSkewTableau :: [SkewTableau Int]+test_dualSkewTableau = bad where +  ps = allPartitions 11+  bad = [ st +        | p<-ps , q<-ps +        , (q `isSubPartitionOf` p) +        , let sp = mkSkewPartition (p,q) +        , let st = fillSkewPartitionWithRowWord sp [1..] +        , dualSkewTableau (dualSkewTableau st) /= st+        ]+-}++--------------------------------------------------------------------------------++-- | Semi-standard skew tableaux filled with numbers @[1..n]@+semiStandardSkewTableaux :: Int -> SkewPartition -> [SkewTableau Int]+semiStandardSkewTableaux n (SkewPartition abs) = map SkewTableau stuff where++  stuff = worker as bs ds (repeat 1) +  (as,bs) = unzip abs+  ds = diffSequence as+  +  -- | @worker inner outerMinusInner innerdiffs lowerbound+  worker :: [Int] -> [Int] -> [Int] -> [Int] -> [[(Int,[Int])]]+  worker (a:as) (b:bs) (d:ds) lb = [ (a,this):rest +                                   | this <- row b 1 lb +                                   , let lb' = (replicate d 1 ++ map (+1) this) +                                   , rest <- worker as bs ds lb' ] +  worker []     _      _      _  = [ [] ]++  -- @row length minimum lowerbound@+  row 0  _  _       = [[]]+  row _  _  []      = []+  row !k !m (!a:as) = [ x:xs | x <- [(max a m)..n] , xs <- row (k-1) x as ] ++{-+-- | from a sequence @[a1,a2,..,an]@ computes the sequence of differences+-- @[a1-a2,a2-a3,...,an-0]@+diffSequence :: [Int] -> [Int]+diffSequence = go where+  go (x:ys@(y:_)) = (x-y) : go ys +  go [x] = [x]+  go []  = []+-}++--------------------------------------------------------------------------------++asciiSkewTableau :: Show a => SkewTableau a -> ASCII+asciiSkewTableau = asciiSkewTableau' "." EnglishNotation++asciiSkewTableau' +  :: Show a+  => 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+  stuff = case orient of+    EnglishNotation    -> es+    EnglishNotationCCW -> reverse (transpose es)+    FrenchNotation     -> reverse es+  inner = asciiFromString innerstr+  es = [ replicate a inner ++ map asciiShow xs | (a,xs) <- axs ]++instance Show a => DrawASCII (SkewTableau a) where+  ascii = asciiSkewTableau++--------------------------------------------------------------------------------++-- | The reversed rows, concatenated+skewTableauRowWord :: SkewTableau a -> [a]+skewTableauRowWord (SkewTableau axs) = concatMap (reverse . snd) axs++-- | The reversed rows, concatenated+skewTableauColumnWord :: SkewTableau a -> [a]+skewTableauColumnWord = skewTableauRowWord . dualSkewTableau++-- | Fills a skew partition with content, in row word order +fillSkewPartitionWithRowWord :: SkewPartition -> [a] -> SkewTableau a+fillSkewPartitionWithRowWord (SkewPartition abs) xs = SkewTableau $ go abs xs where+  go ((b,a):rest) xs = let (ys,zs) = splitAt a xs in (b,reverse ys) : go rest zs+  go []           xs = []++-- | Fills a skew partition with content, in column word order +fillSkewPartitionWithColumnWord :: SkewPartition -> [a] -> SkewTableau a+fillSkewPartitionWithColumnWord shape content +  = dualSkewTableau +  $ fillSkewPartitionWithRowWord (dualSkewPartition shape) content++--------------------------------------------------------------------------------++-- | If the skew tableau's row word is a lattice word, we can make a partition from its content+skewTableauRowContent :: SkewTableau Int -> Maybe Partition+skewTableauRowContent (SkewTableau axs) = go Map.empty rowword where++  rowword = concatMap (reverse . snd) axs++  finish table = Partition (f 1) where+    f !i = case lkp i of+      0 -> []+      y -> y : f (i+1) +    lkp j = case Map.lookup j table of+      Just k  -> k+      Nothing -> 0++  go :: Map Int Int -> [Int] -> Maybe Partition+  go !table []     = Just (finish table)+  go !table (i:is) =+    if check i+      then go table' is+      else Nothing+    where+      table'  = Map.insertWith (+) i 1 table+      check j = j==1 || cnt (j-1) >= cnt j+      cnt j   = case Map.lookup j table' of+        Just k  -> k+        Nothing -> 0++--------------------------------------------------------------------------------
Math/Combinat/Trees/Binary.hs view
@@ -7,6 +7,7 @@ -- <<svg/bintrees.svg>> -- +{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} module Math.Combinat.Trees.Binary    ( -- * Types     BinTree(..)@@ -267,10 +268,10 @@     new =        {- debug (reverse ls,l,r,rs) $ -}        case l of -	      RightParen -> Just ( ls , LeftParen:RightParen:rs )-	      LeftParen  -> -	        {- debug ("---",reverse ls,l,r,rs) $ -}-	        findj ( lls , [] ) ( reverse (RightParen:rs) , [] ) +        RightParen -> Just ( ls , LeftParen:RightParen:rs )+        LeftParen  -> +          {- debug ("---",reverse ls,l,r,rs) $ -}+          findj ( lls , [] ) ( reverse (RightParen:rs) , [] )    next _ = error "fasc4A_algorithm_P: fatal error shouldn't happen"    findj :: ([Paren],[Paren]) -> ([Paren],[Paren]) -> Maybe ([Paren],[Paren])@@ -278,10 +279,10 @@   findj ( lls@(l:ls) , rs) ( xs , ys ) =      {- debug ((reverse ls,l,rs),(reverse xs,ys)) $ -}     case l of-	    LeftParen  -> case xs of-	      (a:_:as) -> findj ( ls, RightParen:rs ) ( as , LeftParen:a:ys )-	      _ -> findj ( lls, [] ) ( reverse rs ++ xs , ys) -	    RightParen -> Just ( reverse ys ++ xs ++ reverse (LeftParen:rs) ++ ls , [] )+      LeftParen  -> case xs of+        (a:_:as) -> findj ( ls, RightParen:rs ) ( as , LeftParen:a:ys )+        _ -> findj ( lls, [] ) ( reverse rs ++ xs , ys) +      RightParen -> Just ( reverse ys ++ xs ++ reverse (LeftParen:rs) ++ ls , [] )   findj _ _ = error "fasc4A_algorithm_P: fatal error shouldn't happen"      -- | Generates a uniformly random sequence of nested parentheses of length 2n.    @@ -409,5 +410,8 @@   blockWidth ls = case ls of     (l:_) -> length l     []    -> 0++instance DrawASCII (BinTree ()) where+  ascii = asciiBinaryTree_   --------------------------------------------------------------------------------
Math/Combinat/Trees/Nary.hs view
@@ -1,6 +1,7 @@  -- | N-ary trees. +{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} module Math.Combinat.Trees.Nary    (           -- * Regular trees @@ -185,6 +186,9 @@                                   else if bf then "@-"                                              else "+-"                    in  (branch++l) : map (indent++) ls ++ gap++instance DrawASCII (Tree ()) where+  ascii = asciiTreeVertical_  -- | Prints all labels. Example: --
combinat.cabal view
@@ -1,5 +1,5 @@ Name:                combinat-Version:             0.2.7.1+Version:             0.2.7.2 Synopsis:            Generate and manipulate various combinatorial objects. Description:         A collection of functions to generate, count and manipulate                      all kinds of combinatorial objects like partitions, @@ -83,5 +83,5 @@   if flag(withQuickCheck)     cpp-options:         -DQUICKCHECK -  ghc-options:         -Wall -fno-warn-unused-matches+  ghc-options:         -Wall -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports
svg/src/gen_figures.hs view
@@ -1,66 +1,66 @@-
--- | A script to generate the SVG figures in the documentation.
--- We use the @combinat-diagrams@ library for that.
-
-module Main where
-
---------------------------------------------------------------------------------
-
-import Math.Combinat.Partitions.Integer
-import Math.Combinat.Partitions.Plane
-import Math.Combinat.Partitions.NonCrossing
-import Math.Combinat.Tableaux
-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.Tableaux
-import Math.Combinat.Diagrams.LatticePaths
-import Math.Combinat.Diagrams.Trees.Binary
-
-import Diagrams.Core
-import Diagrams.Prelude
-import Diagrams.Backend.SVG
-
---------------------------------------------------------------------------------
-
-export fpath size what = renderSVG fpath size $ pad 1.10 what
-
-vcatSep = vcat' (with & sep .~ 1) 
-hcatSep = hcat' (with & sep .~ 1) 
-
-boxSep m xs = pad 1.05 $ vcatSep $ map hcatSep $ yys where
-  yys = go xs where
-    go [] = []
-    go zs = take m zs : go (drop m zs) 
-
---------------------------------------------------------------------------------
-
-main = do 
-
-  export "plane_partition.svg" (Width 320) $ 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 $
-    NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]
-
-  export "young_tableau.svg" (Width 256) $ drawTableau $ 
-    [ [ 1 , 3 , 4 , 6 , 7 ]
-    , [ 2 , 5 , 8 ,10 ]
-    , [ 9 ]
-    ]
-
-  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
-  -- print (pathHeight path, pathNumberOfZeroTouches path, pathNumberOfPeaks path)
-
-  export "ferrers.svg" (Width 256) $ drawFerrersDiagram' EnglishNotation red True $
-    Partition [8,6,3,3,1]
-
-  export "bintrees.svg" (Width 750) $ boxSep 7 $ map drawBinTree_ (binaryTrees 4)
-
---------------------------------------------------------------------------------
++-- | A script to generate the SVG figures in the documentation.+-- We use the @combinat-diagrams@ library for that.++module Main where++--------------------------------------------------------------------------------++import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Plane+import Math.Combinat.Partitions.NonCrossing+import Math.Combinat.Tableaux+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.Tableaux+import Math.Combinat.Diagrams.LatticePaths+import Math.Combinat.Diagrams.Trees.Binary++import Diagrams.Core+import Diagrams.Prelude+import Diagrams.Backend.SVG++--------------------------------------------------------------------------------++export fpath size what = renderSVG fpath size $ pad 1.10 what++vcatSep = vcat' (with & sep .~ 1) +hcatSep = hcat' (with & sep .~ 1) ++boxSep m xs = pad 1.05 $ vcatSep $ map hcatSep $ yys where+  yys = go xs where+    go [] = []+    go zs = take m zs : go (drop m zs) ++--------------------------------------------------------------------------------++main = do ++  export "plane_partition.svg" (Width 320) $ 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 $+    NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]++  export "young_tableau.svg" (Width 256) $ drawTableau $ +    [ [ 1 , 3 , 4 , 6 , 7 ]+    , [ 2 , 5 , 8 ,10 ]+    , [ 9 ]+    ]++  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+  -- print (pathHeight path, pathNumberOfZeroTouches path, pathNumberOfPeaks path)++  export "ferrers.svg" (Width 256) $ drawFerrersDiagram' EnglishNotation red True $+    Partition [8,6,3,3,1]++  export "bintrees.svg" (Width 750) $ boxSep 7 $ map drawBinTree_ (binaryTrees 4)++--------------------------------------------------------------------------------