-- | The Littlewood-Richardson rule

module Math.Combinat.Tableaux.LittlewoodRichardson 
  ( lrCoeff , lrCoeff'
  , lrMult
  , lrRule  , _lrRule , lrRuleNaive
  , lrScalar , _lrScalar
  ) 
  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 Math.Combinat.Helper

import Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map

--------------------------------------------------------------------------------

-- | Naive (very slow) reference implementation of the Littlewood-Richardson rule, based 
-- on the definition "count the semistandard skew tableaux whose row content is a lattice word"
--
lrRuleNaive :: SkewPartition -> Map Partition Int
lrRuleNaive :: SkewPartition -> Map Partition Int
lrRuleNaive SkewPartition
skew = Map Partition Int
final where
  n :: Int
n     = SkewPartition -> Int
skewPartitionWeight SkewPartition
skew
  ssst :: [SkewTableau Int]
ssst  = Int -> SkewPartition -> [SkewTableau Int]
semiStandardSkewTableaux Int
n SkewPartition
skew 
  final :: Map Partition Int
final = (Map Partition Int -> Partition -> Map Partition Int)
-> Map Partition Int -> [Partition] -> Map Partition Int
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Map Partition Int -> Partition -> Map Partition Int
forall k a. (Ord k, Num a) => Map k a -> k -> Map k a
f Map Partition Int
forall k a. Map k a
Map.empty ([Partition] -> Map Partition Int)
-> [Partition] -> Map Partition Int
forall a b. (a -> b) -> a -> b
$ [Maybe Partition] -> [Partition]
forall a. [Maybe a] -> [a]
catMaybes [ SkewTableau Int -> Maybe Partition
skewTableauRowContent SkewTableau Int
skew | SkewTableau Int
skew <- [SkewTableau Int]
ssst  ]
  f :: Map k a -> k -> Map k a
f Map k a
old k
nu = (a -> a -> a) -> k -> a -> Map k a -> Map k a
forall k a. Ord k => (a -> a -> a) -> k -> a -> Map k a -> Map k a
Map.insertWith a -> a -> a
forall a. Num a => a -> a -> a
(+) k
nu a
1 Map k a
old

--------------------------------------------------------------------------------
-- SKEW EXPANSION

-- | @lrRule@ computes the expansion of a skew Schur function 
-- @s[lambda\/mu]@ via the Littlewood-Richardson rule.
--
-- Adapted from John Stembridge's Maple code: 
-- <http://www.math.lsa.umich.edu/~jrs/software/SFexamples/LR_rule>
--
-- > lrRule (mkSkewPartition (lambda,nu)) == Map.fromList list where
-- >   muw  = weight lambda - weight nu
-- >   list = [ (mu, coeff) 
-- >          | mu <- partitions muw 
-- >          , let coeff = lrCoeff lambda (mu,nu)
-- >          , coeff /= 0
-- >          ] 
--
lrRule :: SkewPartition -> Map Partition Int
lrRule :: SkewPartition -> Map Partition Int
lrRule SkewPartition
skew = Partition -> Partition -> Map Partition Int
_lrRule Partition
lam Partition
mu where
  (Partition
lam,Partition
mu) = SkewPartition -> (Partition, Partition)
fromSkewPartition SkewPartition
skew

-- | @_lrRule lambda mu@ computes the expansion of the skew
-- Schur function @s[lambda\/mu]@ via the Littlewood-Richardson rule.
--
_lrRule :: Partition -> Partition -> Map Partition Int
_lrRule :: Partition -> Partition -> Map Partition Int
_lrRule plam :: Partition
plam@(Partition [Int]
lam) pmu :: Partition
pmu@(Partition [Int]
mu0) = 
  if Bool -> Bool
not (Partition
pmu Partition -> Partition -> Bool
`isSubPartitionOf` Partition
plam) 
    then Map Partition Int
forall k a. Map k a
Map.empty
    else (Map Partition Int -> [Int] -> Map Partition Int)
-> Map Partition Int -> [[Int]] -> Map Partition Int
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Map Partition Int -> [Int] -> Map Partition Int
forall a. Num a => Map Partition a -> [Int] -> Map Partition a
f Map Partition Int
forall k a. Map k a
Map.empty [ [Int]
nu | ([Int]
nu,[Int]
_) <- Int -> Diagram -> [([Int], [Int])]
fillings Int
n Diagram
diagram ]
  where
    f :: Map Partition a -> [Int] -> Map Partition a
f Map Partition a
old [Int]
nu = (a -> a -> a)
-> Partition -> a -> Map Partition a -> Map Partition a
forall k a. Ord k => (a -> a -> a) -> k -> a -> Map k a -> Map k a
Map.insertWith a -> a -> a
forall a. Num a => a -> a -> a
(+) ([Int] -> Partition
Partition [Int]
nu) a
1 Map Partition a
old
    diagram :: Diagram
diagram  = [ (Int
i,Int
j) | (Int
i,Int
a,Int
b) <- [(Int, Int, Int)] -> [(Int, Int, Int)]
forall a. [a] -> [a]
reverse ([Int] -> [Int] -> [Int] -> [(Int, Int, Int)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [Int
1..] [Int]
lam [Int]
mu) , Int
j <- [Int
bInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
a] ]    
    mu :: [Int]
mu       = [Int]
mu0 [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ Int -> [Int]
forall a. a -> [a]
repeat Int
0
    n :: Int
n        = [Int] -> Int
forall a. Num a => [a] -> a
sum' ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ (Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Int]
lam [Int]
mu    -- n == length diagram

{-
LR_rule:=proc(lambda) local l,mu,alpha,beta,i,j,dgrm;
  if not `LR_rule/fit`(lambda,args[2]) then RETURN(0) fi;
  l:=nops(lambda); mu:=[op(args[2]),0$l];
  dgrm:=[seq(seq([i,-j],j=-lambda[i]..-1-mu[i]),i=1..l)];
  if nargs>2 then alpha:=args[3];
    if nargs>3 then beta:=args[4] else beta:=[] fi;
    if not `LR_rule/fit`(alpha,beta) then RETURN(0) fi;
    l:=convert([op(lambda),op(beta)],`+`);
    if l<>convert([op(alpha),op(mu)],`+`) then RETURN(0) fi;
    nops(LR_fillings(dgrm,[alpha,beta]))
  else
    convert([seq(s[op(i[1])],i=LR_fillings(dgrm))],`+`)
  fi
end;
-}

--------------------------------------------------------------------------------

-- | A filling is a pair consisting a shape @nu@ and a lattice permutation @lp@
type Filling = ( [Int] , [Int] )

-- | A diagram is a set of boxes in a skew shape (in the right order)
type Diagram = [ (Int,Int) ]

-- | Note: we use reverse ordering in Diagrams compared the Stembridge's code.
-- Also, for performance reasons, we need the length of the diagram
fillings :: Int -> Diagram -> [Filling]
fillings :: Int -> Diagram -> [([Int], [Int])]
fillings Int
_ []                   = [ ([],[]) ]
fillings Int
n diagram :: Diagram
diagram@((Int
x,Int
y):Diagram
rest) = (([Int], [Int]) -> [([Int], [Int])])
-> [([Int], [Int])] -> [([Int], [Int])]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Int -> Int -> ([Int], [Int]) -> [([Int], [Int])]
nextLetter Int
lower Int
upper) (Int -> Diagram -> [([Int], [Int])]
fillings (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Diagram
rest) where
  upper :: Int
upper = case ((Int, Int) -> Bool) -> Diagram -> Maybe Int
forall a. (a -> Bool) -> [a] -> Maybe Int
findIndex ((Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
==(Int
x  ,Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Diagram
diagram of { Just Int
j -> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j ; Maybe Int
Nothing -> Int
0 }
  lower :: Int
lower = case ((Int, Int) -> Bool) -> Diagram -> Maybe Int
forall a. (a -> Bool) -> [a] -> Maybe Int
findIndex ((Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
==(Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
y  )) Diagram
diagram of { Just Int
j -> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j ; Maybe Int
Nothing -> Int
0 }

{-
LR_fillings:=proc(dgrm) local n,x,upper,lower;
  if dgrm=[] then
    if nargs=1 then x:=[] else x:=args[2][2] fi;
    RETURN([[x,[]]])
  fi;
  n:=nops(dgrm); x:=dgrm[n];
  if not member([x[1],x[2]+1],dgrm,'upper') then upper:=0 fi;
  if not member([x[1]-1,x[2]],dgrm,'lower') then lower:=0 fi;
  if nargs=1 then
    map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)]),lower,upper)
  else
    map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)],args[2]),
      lower,upper,args[2][1])
  fi;
end:
-}

--------------------------------------------------------------------------------

nextLetter :: Int -> Int -> Filling -> [Filling]
nextLetter :: Int -> Int -> ([Int], [Int]) -> [([Int], [Int])]
nextLetter Int
lower Int
upper ([Int]
nu,[Int]
lpart) = [([Int], [Int])]
stuff where
  stuff :: [([Int], [Int])]
stuff = [ ( Int -> [Int] -> [Int]
incr Int
i [Int]
shape , [Int]
lpart[Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++[Int
i] ) | Int
i<-[Int]
nlist ] 
  shape :: [Int]
shape = [Int]
nu [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
0] 
  lb :: Int
lb = if Int
lowerInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0
    then [Int]
lpart [Int] -> Int -> Int
forall a. [a] -> Int -> a
!! (Int
lowerInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
    else Int
0
  ub :: Int
ub = if Int
upperInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0 
    then Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
shape) ([Int]
lpart [Int] -> Int -> Int
forall a. [a] -> Int -> a
!! (Int
upperInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))  
    else      [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
shape

  nlist :: [Int]
nlist = (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0) ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Int
f [Int
lbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
ub] 
  f :: Int -> Int
f Int
j   = if Int
jInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
1 Bool -> Bool -> Bool
|| [Int]
shape[Int] -> Int -> Int
forall a. [a] -> Int -> a
!!(Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> [Int]
shape[Int] -> Int -> Int
forall a. [a] -> Int -> a
!!(Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) then Int
j else Int
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]
  incr :: Int -> [Int] -> [Int]
incr Int
i (Int
x:[Int]
xs) = case Int
i of
    Int
0 ->         [Int] -> [Int]
finish (Int
xInt -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int]
xs)
    Int
1 -> (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int] -> [Int]
finish [Int]
xs
    Int
_ -> Int
x     Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: Int -> [Int] -> [Int]
incr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) [Int]
xs
  incr Int
_ []     = []

  -- removes tailing zeros
  finish :: [Int] -> [Int]
  finish :: [Int] -> [Int]
finish (Int
x:[Int]
xs) = if Int
xInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0 then Int
x Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int] -> [Int]
finish [Int]
xs else []    
  finish []     = [] 

{-
`LR/nextletter`:=proc(T) local shape,lp,lb,ub,i,nl;
  shape:=[op(T[1]),0]; lp:=T[2]; ub:=nops(shape);
  if nargs>3 then ub:=min(ub,nops(args[4])) fi;
  if args[2]=0 then lb:=0 else lb:=lp[args[2]] fi;
  if args[3]>0 then ub:=min(lp[args[3]],ub) fi;
  if nargs<4 then
    nl:=map(proc(x,y) if x=1 or y[x-1]>y[x] then x fi end,[$lb+1..ub],shape)
  else
    nl:=map(proc(x,y,z) if y[x]<z[x] and (x=1 or y[x-1]>y[x]) then x fi end,
      [$lb+1..ub],shape,args[4])
  fi;
  nl:=[seq([subsop(i=shape[i]+1,shape),[op(lp),i]],i=nl)];
  op(subs(0=NULL,nl))
end:
-}

--------------------------------------------------------------------------------
-- COEFF

-- | @lrCoeff lam (mu,nu)@ computes the coressponding Littlewood-Richardson coefficients.
-- This is also the coefficient of @s[lambda]@ in the product @s[mu]*s[nu]@
--
-- Note: This is much slower than using 'lrRule' or 'lrMult' to compute several coefficients
-- at the same time!
lrCoeff :: Partition -> (Partition,Partition) -> Int
lrCoeff :: Partition -> (Partition, Partition) -> Int
lrCoeff Partition
lam (Partition
mu,Partition
nu) = 
  if Partition
nu Partition -> Partition -> Bool
`isSubPartitionOf` Partition
lam
    then SkewPartition -> SkewPartition -> Int
lrScalar ((Partition, Partition) -> SkewPartition
mkSkewPartition (Partition
lam,Partition
nu)) ((Partition, Partition) -> SkewPartition
mkSkewPartition (Partition
mu,Partition
emptyPartition))
    else Int
0

-- | @lrCoeff (lam\/nu) mu@ computes the coressponding Littlewood-Richardson coefficients.
-- This is also the coefficient of @s[mu]@ in the product @s[lam\/nu]@
--
-- Note: This is much slower than using 'lrRule' or 'lrMult' to compute several coefficients
-- at the same time!
lrCoeff' :: SkewPartition -> Partition -> Int
lrCoeff' :: SkewPartition -> Partition -> Int
lrCoeff' SkewPartition
skew Partition
p = SkewPartition -> SkewPartition -> Int
lrScalar SkewPartition
skew ((Partition, Partition) -> SkewPartition
mkSkewPartition (Partition
p,Partition
emptyPartition))
  
--------------------------------------------------------------------------------
-- SCALAR PRODUCT

-- | @lrScalar (lambda\/mu) (alpha\/beta)@ computes the scalar product of the two skew
-- Schur functions @s[lambda\/mu]@ and @s[alpha\/beta]@ via the Littlewood-Richardson rule.
--
-- Adapted from John Stembridge Maple code: 
-- <http://www.math.lsa.umich.edu/~jrs/software/SFexamples/LR_rule>
--
lrScalar :: SkewPartition -> SkewPartition -> Int
lrScalar :: SkewPartition -> SkewPartition -> Int
lrScalar SkewPartition
lambdaMu SkewPartition
alphaBeta = (Partition, Partition) -> (Partition, Partition) -> Int
_lrScalar (SkewPartition -> (Partition, Partition)
fromSkewPartition SkewPartition
lambdaMu) (SkewPartition -> (Partition, Partition)
fromSkewPartition SkewPartition
alphaBeta)

_lrScalar :: (Partition,Partition) -> (Partition,Partition) -> Int
_lrScalar :: (Partition, Partition) -> (Partition, Partition) -> Int
_lrScalar ( plam :: Partition
plam@(  Partition [Int]
lam  ) , pmu :: Partition
pmu@(  Partition [Int]
mu0 ) ) 
          ( palpha :: Partition
palpha@(Partition [Int]
alpha) , pbeta :: Partition
pbeta@(Partition [Int]
beta) ) = 
  if    Bool -> Bool
not (Partition
pmu   Partition -> Partition -> Bool
`isSubPartitionOf` Partition
plam  ) 
     Bool -> Bool -> Bool
|| Bool -> Bool
not (Partition
pbeta Partition -> Partition -> Bool
`isSubPartitionOf` Partition
palpha) 
     Bool -> Bool -> Bool
|| ([Int] -> Int
forall a. Num a => [a] -> a
sum' [Int]
lam Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Int] -> Int
forall a. Num a => [a] -> a
sum' [Int]
beta) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= ([Int] -> Int
forall a. Num a => [a] -> a
sum' [Int]
alpha Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Int] -> Int
forall a. Num a => [a] -> a
sum' [Int]
mu0)     -- equivalent to (lambda-mu) /= (alpha-beta)
    then Int
0
    else [([Int], [Int])] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([([Int], [Int])] -> Int) -> [([Int], [Int])] -> Int
forall a b. (a -> b) -> a -> b
$ Int -> Diagram -> ([Int], [Int]) -> [([Int], [Int])]
fillings' Int
n Diagram
diagram ([Int]
alpha,[Int]
beta) 
  where
    f :: Map Partition a -> [Int] -> Map Partition a
f Map Partition a
old [Int]
nu = (a -> a -> a)
-> Partition -> a -> Map Partition a -> Map Partition a
forall k a. Ord k => (a -> a -> a) -> k -> a -> Map k a -> Map k a
Map.insertWith a -> a -> a
forall a. Num a => a -> a -> a
(+) ([Int] -> Partition
Partition [Int]
nu) a
1 Map Partition a
old
    diagram :: Diagram
diagram  = [ (Int
i,Int
j) | (Int
i,Int
a,Int
b) <- [(Int, Int, Int)] -> [(Int, Int, Int)]
forall a. [a] -> [a]
reverse ([Int] -> [Int] -> [Int] -> [(Int, Int, Int)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [Int
1..] [Int]
lam [Int]
mu) , Int
j <- [Int
bInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
a] ]    
    mu :: [Int]
mu       = [Int]
mu0 [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ Int -> [Int]
forall a. a -> [a]
repeat Int
0
    n :: Int
n        = [Int] -> Int
forall a. Num a => [a] -> a
sum' ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ (Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Int]
lam [Int]
mu    -- n == length diagram

{-
LR_rule:=proc(lambda) local l,mu,alpha,beta,i,j,dgrm;
  if not `LR_rule/fit`(lambda,args[2]) then RETURN(0) fi;
  l:=nops(lambda); mu:=[op(args[2]),0$l];
  dgrm:=[seq(seq([i,-j],j=-lambda[i]..-1-mu[i]),i=1..l)];
  if nargs>2 then alpha:=args[3];
    if nargs>3 then beta:=args[4] else beta:=[] fi;
    if not `LR_rule/fit`(alpha,beta) then RETURN(0) fi;
    l:=convert([op(lambda),op(beta)],`+`);
    if l<>convert([op(alpha),op(mu)],`+`) then RETURN(0) fi;
    nops(LR_fillings(dgrm,[alpha,beta]))
  else
    convert([seq(s[op(i[1])],i=LR_fillings(dgrm))],`+`)
  fi
end;
-}

--------------------------------------------------------------------------------

-- | Note: we use reverse ordering in Diagrams compared the Stembridge's code.
-- Also, for performance reasons, we need the length of the diagram
fillings' :: Int -> Diagram -> ([Int],[Int]) -> [Filling]
fillings' :: Int -> Diagram -> ([Int], [Int]) -> [([Int], [Int])]
fillings' Int
_         []                     ([Int]
alpha,[Int]
beta) = [ ([Int]
beta,[]) ]
fillings' Int
n diagram :: Diagram
diagram@((Int
x,Int
y):Diagram
rest) alphaBeta :: ([Int], [Int])
alphaBeta@([Int]
alpha,[Int]
beta) = [([Int], [Int])]
stuff where
  stuff :: [([Int], [Int])]
stuff = (([Int], [Int]) -> [([Int], [Int])])
-> [([Int], [Int])] -> [([Int], [Int])]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Int -> Int -> [Int] -> ([Int], [Int]) -> [([Int], [Int])]
nextLetter' Int
lower Int
upper [Int]
alpha) (Int -> Diagram -> ([Int], [Int]) -> [([Int], [Int])]
fillings' (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Diagram
rest ([Int], [Int])
alphaBeta) 
  upper :: Int
upper = case ((Int, Int) -> Bool) -> Diagram -> Maybe Int
forall a. (a -> Bool) -> [a] -> Maybe Int
findIndex ((Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
==(Int
x  ,Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Diagram
diagram of { Just Int
j -> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j ; Maybe Int
Nothing -> Int
0 }
  lower :: Int
lower = case ((Int, Int) -> Bool) -> Diagram -> Maybe Int
forall a. (a -> Bool) -> [a] -> Maybe Int
findIndex ((Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
==(Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
y  )) Diagram
diagram of { Just Int
j -> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j ; Maybe Int
Nothing -> Int
0 }

{-
LR_fillings:=proc(dgrm) local n,x,upper,lower;
  if dgrm=[] then
    if nargs=1 then x:=[] else x:=args[2][2] fi;
    RETURN([[x,[]]])
  fi;
  n:=nops(dgrm); x:=dgrm[n];
  if not member([x[1],x[2]+1],dgrm,'upper') then upper:=0 fi;
  if not member([x[1]-1,x[2]],dgrm,'lower') then lower:=0 fi;
  if nargs=1 then
    map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)]),lower,upper)
  else
    map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)],args[2]),
      lower,upper,args[2][1])
  fi;
end:
-}

--------------------------------------------------------------------------------

nextLetter' :: Int -> Int -> [Int] -> Filling -> [Filling]
nextLetter' :: Int -> Int -> [Int] -> ([Int], [Int]) -> [([Int], [Int])]
nextLetter' Int
lower Int
upper [Int]
alpha ([Int]
nu,[Int]
lpart) = [([Int], [Int])]
stuff where
  stuff :: [([Int], [Int])]
stuff = [ ( Int -> [Int] -> [Int]
incr Int
i [Int]
shape , [Int]
lpart[Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++[Int
i] ) | Int
i<-[Int]
nlist ] 
  shape :: [Int]
shape = [Int]
nu [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
0] 
  lb :: Int
lb = if Int
lowerInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0
    then [Int]
lpart [Int] -> Int -> Int
forall a. [a] -> Int -> a
!! (Int
lowerInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
    else Int
0
  ub1 :: Int
ub1 = if Int
upperInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0 
    then Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
shape) ([Int]
lpart [Int] -> Int -> Int
forall a. [a] -> Int -> a
!! (Int
upperInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))  
    else      [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
shape
  ub :: Int
ub = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
alpha) Int
ub1
  nlist :: [Int]
nlist = (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0) ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Int
f [Int
lbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
ub] 
  f :: Int -> Int
f Int
j   = if (        [Int]
shape[Int] -> Int -> Int
forall a. [a] -> Int -> a
!!(Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [Int]
alpha[Int] -> Int -> Int
forall a. [a] -> Int -> a
!!(Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) Bool -> Bool -> Bool
&&
             (Int
jInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
1 Bool -> Bool -> Bool
|| [Int]
shape[Int] -> Int -> Int
forall a. [a] -> Int -> a
!!(Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> [Int]
shape[Int] -> Int -> Int
forall a. [a] -> Int -> a
!!(Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) 
          then Int
j 
          else Int
0

  -- increments the i-th element (starting from 1)
  incr :: Int -> [Int] -> [Int]
  incr :: Int -> [Int] -> [Int]
incr Int
i (Int
x:[Int]
xs) = case Int
i of
    Int
0 ->         [Int] -> [Int]
finish (Int
xInt -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int]
xs)
    Int
1 -> (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int] -> [Int]
finish [Int]
xs
    Int
_ -> Int
x     Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: Int -> [Int] -> [Int]
incr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) [Int]
xs
  incr Int
_ []     = []

  -- removes tailing zeros
  finish :: [Int] -> [Int]
  finish :: [Int] -> [Int]
finish (Int
x:[Int]
xs) = if Int
xInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0 then Int
x Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int] -> [Int]
finish [Int]
xs else []    
  finish []     = [] 

{-
`LR/nextletter`:=proc(T) local shape,lp,lb,ub,i,nl;
  shape:=[op(T[1]),0]; lp:=T[2]; ub:=nops(shape);
  if nargs>3 then ub:=min(ub,nops(args[4])) fi;
  if args[2]=0 then lb:=0 else lb:=lp[args[2]] fi;
  if args[3]>0 then ub:=min(lp[args[3]],ub) fi;
  if nargs<4 then
    nl:=map(proc(x,y) if x=1 or y[x-1]>y[x] then x fi end,[$lb+1..ub],shape)
  else
    nl:=map(proc(x,y,z) if y[x]<z[x] and (x=1 or y[x-1]>y[x]) then x fi end,
      [$lb+1..ub],shape,args[4])
  fi;
  nl:=[seq([subsop(i=shape[i]+1,shape),[op(lp),i]],i=nl)];
  op(subs(0=NULL,nl))
end:
-}

--------------------------------------------------------------------------------
-- MULTIPLICATION

type Part = [Int]

-- | Computes the expansion of the product of Schur polynomials @s[mu]*s[nu]@ using the
-- Littlewood-Richardson rule. Note: this is symmetric in the two arguments.
--
-- Based on the wikipedia article <https://en.wikipedia.org/wiki/Littlewood-Richardson_rule>
--
-- > lrMult mu nu == Map.fromList list where
-- >   lamw = weight nu + weight mu
-- >   list = [ (lambda, coeff) 
-- >          | lambda <- partitions lamw 
-- >          , let coeff = lrCoeff lambda (mu,nu)
-- >          , coeff /= 0
-- >          ] 
--
lrMult :: Partition -> Partition -> Map Partition Int
lrMult :: Partition -> Partition -> Map Partition Int
lrMult pmu :: Partition
pmu@(Partition [Int]
mu) pnu :: Partition
pnu@(Partition [Int]
nu) = Map Partition Int
result where
  result :: Map Partition Int
result = (Map Partition Int -> [Int] -> Map Partition Int)
-> Map Partition Int -> [[Int]] -> Map Partition Int
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Map Partition Int -> [Int] -> Map Partition Int
forall a. Num a => Map Partition a -> [Int] -> Map Partition a
add Map Partition Int
forall k a. Map k a
Map.empty ([Int] -> [Int] -> [[Int]]
addMu [Int]
mu [Int]
nu) where
  add :: Map Partition a -> [Int] -> Map Partition a
add !Map Partition a
old [Int]
lambda = (a -> a -> a)
-> Partition -> a -> Map Partition a -> Map Partition a
forall k a. Ord k => (a -> a -> a) -> k -> a -> Map k a -> Map k a
Map.insertWith a -> a -> a
forall a. Num a => a -> a -> a
(+) ([Int] -> Partition
Partition [Int]
lambda) a
1 Map Partition a
old

-- | This basically lists all the outer shapes (with multiplicities) which can be result from the LR rule
addMu :: Part -> Part -> [Part]
addMu :: [Int] -> [Int] -> [[Int]]
addMu [Int]
mu [Int]
part = [Int] -> [Int] -> [Int] -> [Int] -> [[Int]]
forall a. [Int] -> [a] -> [Int] -> [Int] -> [[Int]]
go [Int]
ubs0 [Int]
mu [Int]
dmu [Int]
part where

  go :: [Int] -> [a] -> [Int] -> [Int] -> [[Int]]
go [Int]
_   []     [Int]
_      [Int]
part = [[Int]
part]
  go [Int]
ubs (a
m:[a]
ms) (Int
d:[Int]
ds) [Int]
part = [[[Int]]] -> [[Int]]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [ [Int] -> [a] -> [Int] -> [Int] -> [[Int]]
go (Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
drop Int
d [Int]
ubs') [a]
ms [Int]
ds [Int]
part' | ([Int]
ubs',[Int]
part') <- [Int] -> [Int] -> [([Int], [Int])]
addRowOf [Int]
ubs [Int]
part ]

  ubs0 :: [Int]
ubs0 = Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
take ([Int] -> Int
headOrZero [Int]
mu) [[Int] -> Int
headOrZero [Int]
part Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1..]
  dmu :: [Int]
dmu  = [Int] -> [Int]
diffSeq [Int]
mu


-- | Adds a full row of @(length pcols)@ boxes containing to a partition, with
-- pcols being the upper bounds of the columns, respectively. We also return the
-- newly added columns
addRowOf :: [Int] -> Part -> [([Int],Part)]
addRowOf :: [Int] -> [Int] -> [([Int], [Int])]
addRowOf [Int]
pcols [Int]
part = Int -> [Int] -> [Int] -> [Int] -> [([Int], [Int])]
go Int
0 [Int]
pcols [Int]
part [] where
  go :: Int -> [Int] -> [Int] -> [Int] -> [([Int], [Int])]
go !Int
lb []        [Int]
p [Int]
ncols = [([Int] -> [Int]
forall a. [a] -> [a]
reverse [Int]
ncols , [Int]
p)]
  go !Int
lb (!Int
ub:[Int]
ubs) [Int]
p [Int]
ncols = [[([Int], [Int])]] -> [([Int], [Int])]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [ Int -> [Int] -> [Int] -> [Int] -> [([Int], [Int])]
go Int
col [Int]
ubs ((Int, Int) -> [Int] -> [Int]
addBox (Int, Int)
ij [Int]
p) (Int
colInt -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int]
ncols) | ij :: (Int, Int)
ij@(Int
row,Int
col) <- Int -> Int -> [Int] -> Diagram
newBoxes (Int
lbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
ub [Int]
p ]

-- | Returns the (row,column) pairs of the new boxes which 
-- can be added to the given partition with the given column bounds
-- and the 1-Rieri rule 
newBoxes :: Int -> Int -> Part -> [(Int,Int)]
newBoxes :: Int -> Int -> [Int] -> Diagram
newBoxes Int
lb Int
ub [Int]
part = Diagram -> Diagram
forall a. [a] -> [a]
reverse (Diagram -> Diagram) -> Diagram -> Diagram
forall a b. (a -> b) -> a -> b
$ [Int] -> [Int] -> Int -> Diagram
forall a. [a] -> [Int] -> Int -> [(a, Int)]
go [Int
1..] [Int]
part ([Int] -> Int
headOrZero [Int]
part Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) where
  go :: [a] -> [Int] -> Int -> [(a, Int)]
go (!a
i:[a]
_ ) []      !Int
lp
    | Int
lb Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1 Bool -> Bool -> Bool
&& Int
1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
ub Bool -> Bool -> Bool
&& Int
lp Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0  =  [(a
i,Int
1)]
    | Bool
otherwise                     =  []
  go (!a
i:[a]
is) (!Int
j:[Int]
js) !Int
lp 
    | Int
j1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  Int
lb            =  []
    | Int
j1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
ub Bool -> Bool -> Bool
&& Int
lp Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
j  =  (a
i,Int
j1) (a, Int) -> [(a, Int)] -> [(a, Int)]
forall a. a -> [a] -> [a]
: [a] -> [Int] -> Int -> [(a, Int)]
go [a]
is [Int]
js Int
j       
    | Bool
otherwise           =           [a] -> [Int] -> Int -> [(a, Int)]
go [a]
is [Int]
js Int
j
    where 
      j1 :: Int
j1 = Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1

-- | Adds a box to a partition
addBox :: (Int,Int) -> Part -> Part
addBox :: (Int, Int) -> [Int] -> [Int]
addBox (Int
k,Int
_) [Int]
part = Int -> [Int] -> [Int]
forall a. Num a => Int -> [a] -> [a]
go Int
1 [Int]
part where
  go :: Int -> [a] -> [a]
go !Int
i (a
p:[a]
ps) = if Int
iInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
k then (a
pa -> a -> a
forall a. Num a => a -> a -> a
+a
1)a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
ps else a
p a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> [a] -> [a]
go (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) [a]
ps
  go !Int
i []     = if Int
iInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
k then [a
1] else [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"addBox: shouldn't happen"

-- | Safe head defaulting to zero           
headOrZero :: [Int] -> Int
headOrZero :: [Int] -> Int
headOrZero [Int]
xs = case [Int]
xs of 
  (!Int
x:[Int]
_) -> Int
x
  []     -> Int
0

-- | Computes the sequence of differences from a partition (including the last difference to zero)
diffSeq :: Part -> [Int]
diffSeq :: [Int] -> [Int]
diffSeq = [Int] -> [Int]
forall a. Num a => [a] -> [a]
go where
  go :: [a] -> [a]
go (a
p:ps :: [a]
ps@(a
q:[a]
_)) = (a
pa -> a -> a
forall a. Num a => a -> a -> a
-a
q) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
go [a]
ps
  go [a
p]          = [a
p]
  go []           = []

--------------------------------------------------------------------------------