diff --git a/Math/Combinat/ASCII.hs b/Math/Combinat/ASCII.hs --- a/Math/Combinat/ASCII.hs +++ b/Math/Combinat/ASCII.hs @@ -1,296 +1,300 @@ - --- | A mini-DSL for ASCII drawing of structures. --- --- --- From some structures there is also Graphviz and\/or @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@ +-- () 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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Compositions.hs b/Math/Combinat/Compositions.hs --- a/Math/Combinat/Compositions.hs +++ b/Math/Combinat/Compositions.hs @@ -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) diff --git a/Math/Combinat/FreeGroups.hs b/Math/Combinat/FreeGroups.hs --- a/Math/Combinat/FreeGroups.hs +++ b/Math/Combinat/FreeGroups.hs @@ -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 ) diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs --- a/Math/Combinat/Helper.hs +++ b/Math/Combinat/Helper.hs @@ -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) diff --git a/Math/Combinat/LatticePaths.hs b/Math/Combinat/LatticePaths.hs --- a/Math/Combinat/LatticePaths.hs +++ b/Math/Combinat/LatticePaths.hs @@ -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: --- --- <> --- - -{-# 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: +-- +-- <> +-- + +{-# 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 + +-------------------------------------------------------------------------------- + diff --git a/Math/Combinat/Numbers/Series.hs b/Math/Combinat/Numbers/Series.hs --- a/Math/Combinat/Numbers/Series.hs +++ b/Math/Combinat/Numbers/Series.hs @@ -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 diff --git a/Math/Combinat/Partitions/Integer.hs b/Math/Combinat/Partitions/Integer.hs --- a/Math/Combinat/Partitions/Integer.hs +++ b/Math/Combinat/Partitions/Integer.hs @@ -526,6 +526,9 @@ EnglishNotationCCW -> reverse $ fromPartition $ dualPartition part FrenchNotation -> reverse $ fromPartition $ part +instance DrawASCII Partition where + ascii = asciiFerrersDiagram + -------------------------------------------------------------------------------- {- diff --git a/Math/Combinat/Partitions/Multiset.hs b/Math/Combinat/Partitions/Multiset.hs --- a/Math/Combinat/Partitions/Multiset.hs +++ b/Math/Combinat/Partitions/Multiset.hs @@ -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 + -------------------------------------------------------------------------------- diff --git a/Math/Combinat/Partitions/NonCrossing.hs b/Math/Combinat/Partitions/NonCrossing.hs --- a/Math/Combinat/Partitions/NonCrossing.hs +++ b/Math/Combinat/Partitions/NonCrossing.hs @@ -1,205 +1,205 @@ - --- | Non-crossing partitions. --- --- See eg. --- --- 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 --- --- <> --- --- 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. +-- +-- 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 +-- +-- <> +-- +-- 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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Partitions/Plane.hs b/Math/Combinat/Partitions/Plane.hs --- a/Math/Combinat/Partitions/Plane.hs +++ b/Math/Combinat/Partitions/Plane.hs @@ -1,116 +1,116 @@ - --- | Plane partitions. See eg. --- --- Plane partitions are encoded as lists of lists of Z heights. For example the plane --- partition in the picture --- --- <> --- --- 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. +-- +-- Plane partitions are encoded as lists of lists of Z heights. For example the plane +-- partition in the picture +-- +-- <> +-- +-- 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 ] + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Partitions/Set.hs b/Math/Combinat/Partitions/Set.hs --- a/Math/Combinat/Partitions/Set.hs +++ b/Math/Combinat/Partitions/Set.hs @@ -1,99 +1,99 @@ - --- | Set partitions. --- --- See eg. --- - -{-# 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. +-- + +{-# 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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Partitions/Skew.hs b/Math/Combinat/Partitions/Skew.hs --- a/Math/Combinat/Partitions/Skew.hs +++ b/Math/Combinat/Partitions/Skew.hs @@ -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 + +-------------------------------------------------------------------------------- + diff --git a/Math/Combinat/Partitions/Vector.hs b/Math/Combinat/Partitions/Vector.hs --- a/Math/Combinat/Partitions/Vector.hs +++ b/Math/Combinat/Partitions/Vector.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs --- a/Math/Combinat/Permutations.hs +++ b/Math/Combinat/Permutations.hs @@ -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` diff --git a/Math/Combinat/Sign.hs b/Math/Combinat/Sign.hs --- a/Math/Combinat/Sign.hs +++ b/Math/Combinat/Sign.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Tableaux.hs b/Math/Combinat/Tableaux.hs --- a/Math/Combinat/Tableaux.hs +++ b/Math/Combinat/Tableaux.hs @@ -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 diff --git a/Math/Combinat/Tableaux/GelfandTsetlin.hs b/Math/Combinat/Tableaux/GelfandTsetlin.hs --- a/Math/Combinat/Tableaux/GelfandTsetlin.hs +++ b/Math/Combinat/Tableaux/GelfandTsetlin.hs @@ -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 --- --- @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 +-- +-- @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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs b/Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs --- a/Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs +++ b/Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs @@ -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 diff --git a/Math/Combinat/Tableaux/LittlewoodRichardson.hs b/Math/Combinat/Tableaux/LittlewoodRichardson.hs --- a/Math/Combinat/Tableaux/LittlewoodRichardson.hs +++ b/Math/Combinat/Tableaux/LittlewoodRichardson.hs @@ -2,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] diff --git a/Math/Combinat/Tableaux/Skew.hs b/Math/Combinat/Tableaux/Skew.hs --- a/Math/Combinat/Tableaux/Skew.hs +++ b/Math/Combinat/Tableaux/Skew.hs @@ -1,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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Trees/Binary.hs b/Math/Combinat/Trees/Binary.hs --- a/Math/Combinat/Trees/Binary.hs +++ b/Math/Combinat/Trees/Binary.hs @@ -7,6 +7,7 @@ -- <> -- +{-# 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_ -------------------------------------------------------------------------------- diff --git a/Math/Combinat/Trees/Nary.hs b/Math/Combinat/Trees/Nary.hs --- a/Math/Combinat/Trees/Nary.hs +++ b/Math/Combinat/Trees/Nary.hs @@ -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: -- diff --git a/combinat.cabal b/combinat.cabal --- a/combinat.cabal +++ b/combinat.cabal @@ -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 diff --git a/svg/src/gen_figures.hs b/svg/src/gen_figures.hs --- a/svg/src/gen_figures.hs +++ b/svg/src/gen_figures.hs @@ -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) + +--------------------------------------------------------------------------------