######################################################################
## StanleySolitaire.txt Save this file as   StanleySolitaire.txt     #
#to use it,                                                          # 
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `StanleySolitaire.txt` <Enter>               #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################



print(`First Written: July 2023: tested for Maple 2020 `):
print(`Version  : June 2024  `):
print():
print(`This is StanleySolitaire.txt, a  Maple package that`):
print(`accompanies Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` In How Many Ways Can You Play Stanley Solitaire?`):

print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):

print(`--------------------------`):
print(`For a list of the supporting functions type: ezra1();`):
print(`For help with a specific procedure type:`):
print(`ezra(ProcedureName);`):
print(`--------------------------`):

print():
print(`--------------------------`):



print(`--------------------------`):
print():
print(`For a list of the  231 functions type: ezra231();`):
print(`For help with a specific procedure type:`):
print(`ezra(ProcedureName);`):

print():
print(`--------------------------`):



print(`--------------------------`):
print():
print(`For a list of  Story functions type: ezraSt();`):
print(`For help with a specific procedure type:`):
print(`ezra(ProcedureName);`):

print():
print(`--------------------------`):


print(`--------------------------`):
print():
print(`For a list of the Operator type functions type: ezraOp();`):
print(`For help with a specific procedure type:`):
print(`ezra(ProcedureName);`):

print():
print(`--------------------------`):


print(`--------------------------`):
print():
print(`For a list of the General tableaux functions type: ezraG();`):
print(`For help with a specific procedure type:`):
print(`ezra(ProcedureName);`):

print():
print(`--------------------------`):





ezraSt:=proc()
if args=NULL then

print(`The Story procedures are: `):
print(` SakajStory, SbkaStory, S231Story `):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: `):
print(` AllBTslow, Arm, ChaBT, Conj, CT, GenBT1, GuessPol, GuessRat,GuessRatM,   IsGood, IsGoodS, IT, IT0, ITr, iITr, IsBa, Kids, Kids1, Lam, Leg, Mu, permsCl, PermsDIV, PosToPer, PrT, RAN, SSnuN,YF `):

else
ezra(args):

fi:


end:



ezra231:=proc()
if args=NULL then

print(`The 231 procedures are: `):
print( ` Check231D1,   Check231D2, S231`):

else
ezra(args):

fi:


end:


ezraOp:=proc()
if args=NULL then

print(`The OPERATORS procedures are: `):
print(` CheckSbka, Hafel2, Ker2, Opek, Opek1j, SakaC, SakajS, Sakaj, Sbka, SbkaS, SbkaC, `):
#print(` Opek10, Opek11, Opek12, Saka, SakaS,  Saka1, Saka1S, Saka2S,  Saka2`):
else
ezra(args):

fi:


end:



ezraG:=proc()
if args=NULL then

print(`The General Tableaux procedures are: `):
print(` HookT, RandTab  `):

else
ezra(args):

fi:


end:


ezraBT:=proc()
if args=NULL then

print(`The Balanced Tableaux procedures are: `):
print(` AllBT, BTnu   `):

else
ezra(args):

fi:


end:



ezra:=proc()
if args=NULL then

print(` StanleySoliare.txt: A Maple package for studying Stanley Solitaire `):

print(`The MAIN procedures are:  `):
print(` GenSS, GoodPerms, GoodPermsS, KidsSS, NRD, S231, SScha, SSnu, SSnuCF3 `):
print(` `):




elif nargs=1 and args[1]=AllBTslow then
print(`AllBTslow(L): All the balanced tableax of shape L. Try:`):
print(`DONE Slowly. Only for checking puropses. Try:`):
print(`AllBTslow([3,2,1]);`):


elif nargs=1 and args[1]=AllBT then
print(`AllBT(L): All the balanced tableax of shape L, via dynamical programming. Try:`):
print(`AllBT([3,2,1]);`):


elif nargs=1 and args[1]=AllBTc then
print(`AllBTc(L): All the balanced tableax of shape L, represented as chains of states. Try:`):
print(`AllBTc([3,2,1]);`):

elif nargs=1 and args[1]=Arm then
print(`Arm(L,c): The arm-length of the cell c in the shape L. Try: `):
print(`Arm([4,3,2],[2,2]);`):


elif nargs=1 and args[1]=BTnu then
print(`BTnu(L): the number of balanced tableax of shape L, via dynamical programming. Try:`):
print(`BTnu([3,2,1]);`):

elif nargs=1 and args[1]=ChaBT then
print(`ChaBT(L,T): All the balanced tableaux chains  chains starting at the state T. Try:`):
print(`ChaBT([3,2,1],[[0,0,0],[0,0],[0]]);`):


elif nargs=1 and args[1]=Check231D1 then
print(`Check231D1(K1,K2): checks the explicit formula for  SSnu([a1+1,a1,a3])/YF([a3,a1+2,a1])  for a1=1..K1, and a3=a1+2..a1+K2`):
print(`Check231D1(5,10);`):

elif nargs=1 and args[1]=Check213D2 then
print(`Check231D2(K): checks the explcit formula for SSnu([a1+2,a1,a3])/YF([a3,a1+2,a1]) for a1=1..K1 and a3=a1+3..a1+K2`):
print(`Check231D2(5,10);`):

elif nargs=1 and args[1]=CheckSbka then
print(`CheckSbka(k,K): Checks that the formula gotten from SbkaS(b,k,a) is indeed equal to SSnu([b,0$k,a]) for 1<=b<=a<=K`):
print(`CheckSbka(4,10);`):

elif nargs=1 and args[1]=Conj then
print(`Conj(L): The conjugate of the partition L. Try:`):
print(`Conj([5,4,3]);`):


elif nargs=1 and args[1]=CT then
print(`CT(F,x): Given a rational function F in the list of variables x, finds the constant term. Try:`):
print(`CT((1-x1/x2)/(x1^7*x2^6)/(1-x1-x2),[x1,x2]);`):


elif nargs=1 and args[1]=Hafel2 then
print(`Hafel2(ope,F,a,b,A,B): applies the operator ope (in terms of the shift operators A, B of a and b respectively)  to the expression F in a, b. Try:`):
print(`Hafel2(A+B,(a+b)!/a!/b!,a,b,A,B);`):

elif nargs=1 and args[1]=HookT then
print(`HookT(T): Given a tableau T outputs the Rank table. Try:`):
print(`HookT([[3,5,6],[2,4],[1]]);`):


elif nargs=1 and args[1]=GenBT1 then
print(`GenBT1(L,k): all the descendants of the balanced tableux state L, Try:`):
print(`GenBT([[0,0,0],[0,0]],2);`):

elif nargs=1 and args[1]=GenSS then
print(`GenSS(L,k): all the descendants of the state L in Stanley Solitaire, L at the k-th generation, Try:`):
print(`GenSS([3,2,0,0,0],2);`):

elif nargs=1 and args[1]=GoodPerms then
print(`GoodPerms(k,K): The set of permutations such that  for any increasing sequence of pos. integers of length k<=K v, SSnu(pi[v])=YF(Rev(v)). Try:`):
print(`GoodPerms(3,6);`):


elif nargs=1 and args[1]=GuessPol then
print(`GuessPol(L,n): inputs a list of pairs [input,output] and tries to guess`):
print(`a polynomial of degree <=nops(L)-3. Try:`):
print(`GuessPol([seq([i,i^3],i=5..15)],x);`):

elif nargs=1 and args[1]=GuessRat then
print(`GuessRat(L,n): inputs a list of pairs [[input,output], ...] a variable name n`):
print(`tries to guess a rational function in the variable n, for it. try:`):
print(`GuessRat([seq([i,(i+1)/(i+2)],i=1..10)]):`):

elif nargs=1 and args[1]=GuessRatM then
print(`GuessRatM(DS,x): Given a  a data set DS, and a list of variables x, tries to fit a rational function to it. Try:`):
print(`GuessRatM({seq(seq([[i1,i2],(i1+i2)/(i1+i2+1)],i2=1..i1),i1=10..20)},[x,y]); `):




elif nargs=1 and args[1]=iITr then
print(`iITr(w): Given an inversion table, finds the permutation. Try:`):
print(`iITr([4,3,2,1]);`):

elif nargs=1 and args[1]=IT then
print(`IT(w): Given a word w of integers, finds the inversion table. Try:`):
print(`IT([4,3,2,1]);`):


elif nargs=1 and args[1]=ITr then
print(`ITr(w): Given a permutation pi, finds the inversion tableTry:`):
print(`ITr([4,3,2,1]);`):


elif nargs=1 and args[1]=IT0 then
print(`IT0(w): Given a word w of integers, finds the reduced inversion table (with 0 removed from the front and back. Try:`):
print(`IT0([4,3,2,1]);`):

elif nargs=1 and args[1]=IsBa then
print(`IsBa(T): Is the tableu T balanced?. try:`):
print(` IsBa([[3,5,6],[2,4],[1]]);`):

elif nargs=1 and args[1]=IsGood then
print(`IsGood(pi): Is the permutation pi good in the sense that SSnu([a[pi[1]],..., a[pi[k]]])=YF([a[k],..., a[1]]) for all decreasing sequences?. Try:`):
print(`IsGood([2,1,3]);`):

elif nargs=1 and args[1]=IsGoodS then
print(`IsGoodS(pi): Is the permutation pi good in Stanley's sense? Try:`):
print(`IsGoodS([2,1,3]);`):

elif nargs=1 and args[1]=Ker2 then
print(`Ker2(x1,x2,z): The constant-termand for SbkaC(b,k,a) (q.v.). Try:`):
print(`Ker2(x1,x2,z)`):

elif nargs=1 and args[1]=Kids then
print(`Kids(T): given a partially filled balanced tableu finds all its continuations. Try:`):
print(`Kids([[0,5,6],[0,0],[0]]);`):

elif nargs=1 and args[1]=Kids1 then
print(`Kids1(T): given a general shape consisting of 0 and n's finds all its legal children`):
print(`Kids1([[0,6,6],[0,0],[0]]);`):

elif nargs=1 and args[1]=KidsSS then
print(`KidsSS(L) All the children of the state L, try:`):
print(`KidsSS([3,2]);`):


elif nargs=1 and args[1]=Lam then
print(`Lam(pi): Given a permutation pi, finds the quantity mu(w) defined on top of p. 367 of Richard Stanley's article "On the Number of Reduced Decompositions of Elements in the Coxeter Group"`):
print(`Eur. J. Comb. 5 (1984). 359-372. Try:`):
print(`Lam([1,4,3,2]);`):

elif nargs=1 and args[1]=Leg then
print(`Leg(L,c): The leg-length of the cell c in the shape L. Try: `):
print(`Leg([4,3,2],[2,2]);`):

elif nargs=1 and args[1]=Mu then
print(`Mu(pi): Given a permutation pi, finds the quantity mu(w) defined on top of p. 367 of Richard Stanley's article "On the Number of Reduced Decompositions of Elements in the Coxeter Group"`):
print(`Eur. J. Comb. 5 (1984). 359-372. Try:`):
print(`Mu([1,4,3,2]);`):


elif nargs=1 and args[1]=NRD then
print(`NRD(pi): The number of reduced decompositions in S_n of the permutation pi (where n=nops(pi)). Try:`):
print(`NRD([5,4,3,2,1]);`):


elif nargs=1 and args[1]=Opek then
print(`Opek(k,A,B): The operator with constant coefficients in the shift operators A and B such that`):
print(`if a>=b`):
print(`SSnu([b,0^k,a])=Opek(k,A,B)YF([a,b]). Try:`):
print(`Opek(4,A,B);`):

elif nargs=1 and args[1]=Opek10 then
print(`Opek10(k,A,B): The operator such that SSnu([a1,0$k,a1])=Opek10(k,A,B)YF([a1,a1]). Try:`):
print(`Opek10(3,A,B):`):


elif nargs=1 and args[1]=Opek11 then
print(`Opek11(k,A,B): The operator such that SSnu([a1,0$k,a1-1])=Opek11(k,A,B)YF([a1,a1-1]). Try:`):
print(`Opek11(3,A,B):`):


elif nargs=1 and args[1]=Opek12 then
print(`Opek12(k,A,B): The operator such that SSnu([a1,0$k,a1-2])=Opek12(k,A,B)YF([a1,a1-2]). Try:`):
print(`Opek12(3,A,B):`):

elif nargs=1 and args[1]=Opek1j then
print(`Opek1j(j,k,A,B): The operator such that SSnu([a1,0$k,a1-j])=Opek1j(j,k,A,B)YF([a1,a1-j]). Try:`):
print(`Opek1j(2,3,A,B):`):


elif nargs=1 and args[1]=PermsCl then
print(`PermsCl(n): inputs a pos. integer n, and outputs a list of sets of length n-1 where the i-th entry is the set of permutations whose inversion table has i items. Try:`):
print(`PermsCl(5);`):

elif nargs=1 and args[1]=PermsDIV then
print(`PermsDIV(n): all the permutations of n that have a decreasing reduced inversion table. Try:, together with it`):
print(`PermsDIV(5);`):

elif nargs=1 and args[1]=PosToPer then
print(`PosToPer(L): Given a starting position in Stanley Solitaire,L, finds the permutation of smallest length such that its inversion table is L followed by 0's. Try:`):
print(`PosToPer([20,10,30]);`):

elif nargs=1 and args[1]=PrT then
print(`PrT(T): prints the tablea T nicely. Try:`):
print(`PrT([[1,2,5],[3,4]]);`):


elif nargs=1 and args[1]=RAN then
print(`RAN(L,a): The rank of the item a in the list L. Try:`):
print(`RAN([5,1,3,2],2);`):

elif nargs=1 and args[1]=RandTab then
print(`RandTab(L): R random tableau of shape L (no restrictions). Try:`):
print(`RandTab([3,3,3]);`):

elif nargs=1 and args[1]=S231 then
print(`S231(a2,a3,a1): The explicit formula for SSnu([a2,a3,a1]) if a1>=a2>=a3. Try:`):
print(`S231(a2,a3,a1); `):
print(`S231(30,20,33); `):

elif nargs=1 and args[1]=S231Story then
print(`S231Story(): A paper about playing Stanley Solitaire with starting position [a2,a3,a1] where a1>=a2>=a3. Try:`):
print(`S231SStory();`):

elif nargs=1 and args[1]=Saka then
print(`Saka(a,k): Computing SSnu([a,0^k,a]) via SakaS(a,k). Try: `):
print(`Saka(10,4);`):


elif nargs=1 and args[1]=SakaS then
print(`SakaS(a,k): The explicit expression for SSnu([a,0^k,a]). Try: `):
print(`SakaS(a,4);`):


elif nargs=1 and args[1]=Saka1S then
print(`Saka1S(a,k): The explicit expression for SSnu([a,0^k,a-1]). Try: `):
print(`Saka1S(a,4);`):


elif nargs=1 and args[1]=SakajS then
print(`SakajS(j,a,k): The explicit expression for SSnu([a,0^k,a-j]) for symbolic a. Try: `):
print(`SakajS(3,a,7);`):

elif nargs=1 and args[1]=SakajStory then
print(`SakajStory(J,K): Explicit expressions in a, for the Number of Ways of Playing Stanley Solitaire starting with position [a,0$k,a-j] where  for k from j from 0 to K and`):
print(`k from j+1 to J+K. Try:`):
print(`SakajStory(6,6);`):



elif nargs=1 and args[1]=Saka1 then
print(`Saka1(a,k): Computing SSnu([a,0^k,a-1]) via Saka1S(a,k). Try: `):
print(`Saka1(10,4);`):


elif nargs=1 and args[1]=Sakaj then
print(`Sakaj(j,a,k): Computing SSnu([a,0^k,a-j]) via SakajS(j,a,k). Try: `):
print(`Sakaj(1,10,4);`):

elif nargs=1 and args[1]=Sbka then
print(`Sbka(b,k,a): if a>=b, it gives SSnu([b,0$k,a]) using the explicit formula SbkaS(b,k,a). Try:`):
print(` [SSnu([5,0$9,2]),Sbka(2,9,5)]; `):

elif nargs=1 and args[1]=SbkaC then
print(`SbkaC(b,k,a): if a>=b, it gives SSnu([b,0$k,a]) using the constant term expression. Try:`):
print(`[SSnu([5,0$9,2]),SbkaC(5,2,9)];`):

elif nargs=1 and args[1]=SbkaS then
print(`SbkaS(b,k,a): The explicit expression for SSnu([b,0^k,a]) if a>=b, and k>=0.  HERE a and b are symbolic and k is numeric. Try:`):
print(`SbkaS(b,3,a);`):

elif nargs=1 and args[1]=SbkaStory then
print(`SbkaStory(K): Explicit expressions in a and b, for the Number of Ways of Playing Stanley Solitaire starting with position [b,0$k,a] where a>b, for k from 0 to K. Do`):
print(`SbkaStory(10):`):

elif nargs=1 and args[1]=SScha then
print(`SScha(L): all the ways to play Stanley Solitaire with starting position L. Try:`):
print(`SScha([4,3]);`):

elif nargs=1 and args[1]=SSnu then
print(`SSnu(L): The number of ways to play Stanley Solitaire with starting position L. Try:`):
print(`SSnu([3,2]);`):

elif nargs=1 and args[1]=SSnuCF3 then
print(`SSnuCF3(L): The closed form expression for SSnu(L) for any list of positive integers [a1,a2,a3]. Try:`):
print(`{seq(seq(seq(SSnuCF3([a1,a2,a3])-SSnu([a1,a2,a3]),a1=1..10),a2=1..10),a3=1..10)};`):

elif nargs=1 and args[1]=SSnuN then
print(`SSnuN(L): Given a list L of non-negative integers L, kicks out the 0-th if they exist, then sorts it in decreasing order, let's call the new list L1, and outputs`):
print(`SSnu(L)/YF(L1). Try`):
print(`SSnuN([2,4,8]);`):

elif nargs=1 and args[1]=YF then
print(`YF(L): The Young-Frobenius formula. Try: `):
print(`YF([4,3,2]);`):




else
 print(`There is no such thing as`, args):

fi:


end:


with(combinat):


#YF(L): The Young-Frobenius formula
YF:=proc(L) local i,j,n,k:
n:=convert(L,`+`):
k:=nops(L):
n!/mul((L[i]+k-i)!,i=1..k)*mul(mul((L[j]+k-j)-(L[i]+k-i),j=1..i-1),i=1..k):
end:

#Conj(L): The conjugate of the shape L. Try:
#Conj([5,5,3]);
Conj:=proc(L) local L1,i,k:
option remember:
if L=[] then
 RETURN([]):
fi:
k:=nops(L):
L1:=L-[1$k]:
for  i from 1 to k while L1[i]>0 do od:
L1:=[op(1..i-1,L1)]:
[k,op(Conj(L1))]:
end:

#Arm(L,c): Given a cell c, finds its arm
Arm:=proc(L,c) local i,j:
i:=c[1]:
j:=c[2]:
if not (i<=nops(L) and  j<=L[i]) then
 RETURN(FAIL):
fi:
L[i]-j+1:
end:


#Arm(L,c): Given a cell c, finds its arm
Arm:=proc(L,c) local i,j:
i:=c[1]:
j:=c[2]:
if not (i<=nops(L) and  j<=L[i]) then
 RETURN(FAIL):
fi:
L[i]-j+1:
end:


#Leg(L,c): Given a cell c, finds its arm
Leg:=proc(L,c) local i,j,L1:
i:=c[1]:
j:=c[2]:
if not (i<=nops(L) and  j<=L[i]) then
 RETURN(FAIL):
fi:
L1:=Conj(L):
L1[j]-i+1:
end:

#RAN(L,a): The rank of the item a in the list L. Try:
#RAN([5,1,3,2],2);
RAN:=proc(L,a) local L1,i:


if not member(a,{op(L)}) then
 RETURN(FAIL):
fi:

L1:=sort(L):

for i from 1 to nops(L) while L1[i]<>a do od:
i:
end:

#IsBa(T): Given a tableau T checks whether it is a balanced tableau. Try:
#IsBa([[3,5,6],[2,4],[1]]);
IsBa:=proc(T) local L,c,i,j,i1,j1,n,H,L1:
L:=[seq(nops(T[i]),i=1..nops(T))]:
L1:=Conj(L):
n:=convert(L,`+`):
if {seq(op(T[i]),i=1..nops(T))}<>{seq(i,i=1..n)} then
 RETURN(false):
fi:

for i from 1 to nops(T) do
for j from 1 to nops(T[i]) do
c:=[i,j]:
H:=[seq(T[i][j1],j1=j..L[i]), seq(T[i1][j],i1=i+1..L1[j])]:
 if RAN(H,T[i][j])<>Leg(L,c) then
  RETURN(false):
 fi:
od:
od:

true:

end:

#AllBTslow(L): All the balanced tableax of shape L. Try:
#DONE Slowly. Only for checking puropses. Try:
#AllBTslow([3,2,1]);

AllBTslow:=proc(L) local n,mu,mu1,co,i,j,gu,T,T1:
n:=convert(L,`+`):
mu:=permute(n):
gu:={}:
for mu1 in mu do
co:=1:
T:=[]:
for i from 1 to nops(L) do
T1:=[]:

for j from 1 to L[i] do
T1:=[op(T1),mu1[co]]:
co:=co+1:
od:
T:=[op(T),T1]:
od:
if IsBa(T) then
 gu:=gu union {T}:
fi:

od:

gu:
end:

#Kids(T): given a partially filled balanced tableu finds all its continuations. Try:
#Kids([[0,5,6],[0,0],[0]]);
Kids:=proc(T) local gu,i,j,m,L,n,L1,i1,j1,H,T1,c:
L:=[seq(nops(T[i]),i=1..nops(T))]:
L1:=Conj(L):
n:=convert(L,`+`):

if {seq(op(T[i]),i=1..nops(T))}={0} then
 m:=n:
else
m:=min({seq(op(T[i]),i=1..nops(T))} minus {0})-1:
fi:

gu:={}:

for i from 1 to nops(L) do
 for j from 1 to L[i] do
 c:=[i,j]:
  if T[i][j]=0 then
    T1:=[op(1..i-1,T),[op(1..j-1,T[i]), m, op(j+1..L[i],T[i])],op(i+1..nops(T),T)]:

 H:=[seq(T1[i][j1],j1=j..L[i]), seq(T1[i1][j],i1=i+1..L1[j])]:
 if RAN(H,T1[i][j])=Leg(L,c) then
  gu:=gu union {T1}:
 fi:
 fi:
od:
od:
gu:
end:



#AllBT(L): All the balanced tableax of shape L. Try:
#AllBT([3,2,1]);
#DONE Faste
AllBT:=proc(L) local gu,n,T,gu1,i:
 n:=convert(L,`+`):
 T:=[seq([0$L[i]],i=1..nops(L))]:
 gu:={T}:

for i from 1 to n do
 gu:={seq(op(Kids(gu1)),gu1 in gu)}:
od:

gu:
end:





#Kids1(T): given a general shape consisting of 0 and n's finds all its legal children
#Kids1([[0,6,6],[0,0],[0]]);
Kids1:=proc(T) local gu,i,j,L,n,L1,i1,j1,H,T1,c:
L:=[seq(nops(T[i]),i=1..nops(T))]:
L1:=Conj(L):
n:=convert(L,`+`):

gu:={}:

for i from 1 to nops(L) do
 for j from 1 to L[i] do
 c:=[i,j]:
  if T[i][j]=0 then
    T1:=[op(1..i-1,T),[op(1..j-1,T[i]), n, op(j+1..L[i],T[i])],op(i+1..nops(T),T)]:

 H:=[seq(T1[i][j1],j1=j..L[i]), seq(T1[i1][j],i1=i+1..L1[j])]:
 if RAN(H,T1[i][j])=Leg(L,c) then
  gu:=gu union {T1}:
 fi:
 fi:
od:
od:
gu:
end:

#ChaBT(L,T): All the balanced tableaux chains  chains starting at T
ChaBT:=proc(L,T) local n,gu,mu,mu1,gu1, gu11,i:
option remember:
n:=convert(L,`+`):
if T=[seq([n$L[i]],i=1..nops(L))] then
  RETURN({[T]}):
fi:

mu:=Kids1(T):

gu:={}:

for mu1 in mu do
 gu1:=ChaBT(L,mu1):
 gu:=gu union {seq([T,op(gu11)],gu11 in gu1)}:
od:
gu:

end:

#AllBTc(L):  All the balanced tableax of shape L in terms of chains
AllBTc:=proc(L) local i:
 ChaBT(L,[seq([0$L[i]],i=1..nops(L))]):
end:



#ChaNu(L,T): the number of balanced tableaux chains  starting at T
ChaNu:=proc(L,T) local n,i,mu,mu1:
option remember:
n:=convert(L,`+`):
if T=[seq([n$L[i]],i=1..nops(L))] then
  RETURN(1):
fi:

mu:=Kids1(T):

add(ChaNu(L,mu1),mu1 in mu):


end:


BTnu:=proc(L)  local i:
 ChaNu(L,[seq([0$L[i]],i=1..nops(L))]):
end:

#PrT(T): prints the tableau T
PrT:=proc(T) local i:
for i from 1 to nops(T) do
lprint(op(T[i])):
od:
end:

##Start Stanley Solitaire


#KidsSS(L) All the children of the state L, try:
#KidsSS([3,2,0,0,0]);
KidsSS:=proc(L) local gu,i,gu1:
option remember:


if L=[] then
 RETURN([]):
fi:

gu:={}:

for i from 1 to nops(L)-1 do
 if L[i]>L[i+1] then
  gu:=[op(gu),[op(1..i-1,L),L[i+1],L[i]-1,op(i+2..nops(L),L)]]:
 fi:
od:


for i from nops(L) to nops(L) do

  gu:=[op(gu), [op(1..i-1,L),0,L[i]-1]]:
od:


[seq(Chop0(gu1),gu1 in gu)]:

end:


#SSnu(L): The number of ways to play Stanley Solitaire with starting position L. Try:
#SSnu([3,2]);

SSnu:=proc(L)  local gu,i:
option remember:

if L<>[] and L[nops(L)]=0 then
 RETURN(FAIL):
fi:

if L=[] then
 RETURN(1):
else
gu:=KidsSS(L):

add(SSnu(gu[i]),i=1..nops(gu)):
fi:

end:


#SScha(L): all the ways to play Stanley Solitaire with starting position L. Try:
#SScha([4,3]);
SScha:=proc(L) local gu,mu,mu1,gu1,gu11:
option remember:

if L=[] then
 RETURN([[L]]):
fi:

gu:=[]:

mu:=KidsSS(L):

for mu1 in mu do
 gu1:=SScha(mu1):
 gu:=[op(gu), seq([L,op(gu11)],gu11 in gu1)]:
od:
gu:
end:


#GenSS(L,k): all the descendants of L at the k-th generation, Try:
#GenSS([3,2,0,0,0],2);
GenSS:=proc(L,k) local gu,gu1:
option remember:
if k=0 then
 RETURN({L}):
fi:

gu:=GenSS(L,k-1):

{seq(op(KidsSS(gu1)),gu1 in gu)}:
end:





#GenBT1(L,k): all the descendants of the state L at the k-th generation. Try:
#GenBT1([[0,0,0],[0,0]],2);
GenBT1:=proc(L,k) local gu,gu1:
option remember:
if k=0 then
 RETURN({L}):
fi:

gu:=GenBT1(L,k-1):

{seq(op(Kids1(gu1)),gu1 in gu)}:
end:



##End Stanley Solitaire



##start General Tableaux

#RandTab(L): R random tableau of shape L (no restrictions). Try:
#RandTab([3,3,3]);
RandTab:=proc(L) local pi,i,j,T,T1,co:
pi:=randperm(convert(L,`+`)):

T:=[]:

co:=0:
for i from 1 to nops(L) do
T1:=[]:
for j from 1 to L[i] do
 co:=co+1:
 T1:=[op(T1),pi[co]]:
od:
T:=[op(T),T1]:
od:

T:
end:








#HookT(T): Given a tableau T outputs the Rank table. Try:
#HookT([[3,5,6],[2,4],[1]]);
HookT:=proc(T) local L,c,i,j,i1,j1,n,H,L1,TAB,TAB1:
L:=[seq(nops(T[i]),i=1..nops(T))]:
L1:=Conj(L):
n:=convert(L,`+`):
if {seq(op(T[i]),i=1..nops(T))}<>{seq(i,i=1..n)} then
 RETURN(FAIL):
fi:

TAB:=[]:
for i from 1 to nops(T) do
TAB1:=[]:
for j from 1 to nops(T[i]) do
 c:=[i,j]:
 H:=[seq(T[i][j1],j1=j..L[i]), seq(T[i1][j],i1=i+1..L1[j])]:
 TAB1:=[op(TAB1),RAN(H,T[i][j])]:
od:
TAB:=[op(TAB),TAB1]:
od:

TAB:

end:


#End General Tableaux


Chop0:=proc(L) local L1,i:

for i from 1 to nops(L) while L[i]=0 do od:
L1:=[op(i..nops(L),L)]:

for i from nops(L1) to 1 by -1 while L1[i]=0 do od:

[op(1..i,L1)]:

end:


#Opek(k,A,B): The operator with constant coefficients in the shift operators A and B such that
#if a>=b
#SSnu([b,0^k,a])=Opek(k,A,B)YF([a,b]). Try:
#Opek(4,A,B);
Opek:=proc(k,A,B) local gu,i:
option remember:
if k=0 then
RETURN(A-A/B):
elif k=1 then
RETURN(A^2-A^2/B):
else
gu:=expand(A*Opek(k-1,A,B)-A/B*Opek(k-2,A,B)):
RETURN(add(factor(coeff(gu,A,i))*A^i,i=ldegree(gu,A)..degree(gu,A))):
fi:
end:


#SbkaS(b,k,a): The explicit expression for SSnu([b,0^k,a]) if a>=b, and k>=0.  FOR SYMBOLIC a and b and numeric k.
#Try:
#SbkaS(b,3,a);
SbkaS:=proc(b,k,a) local A,B:
simplify(Hafel2(Opek(k,A,B),YF([a,b]),a,b,A,B)):
end:


#Sbka(b,k,a): SSnu([b,0$k,a]) according to the formula given by SbkaS(b,k,a)
#Try:
#Sbka(5,3,8);
Sbka:=proc(b,k,a) local A,B:

if not (type(a,integer) and type(b,integer) and type(k,integer) and a>=0 and b>=0 and a>=b and k>=0) then
 RETURN(FAIL):
fi:


eval(subs({A=a,B=b},SbkaS(B,k,A) )):
end:


#CheckSbka(k,K): Checks that the formula gotten from SbkaS(b,k,a) is indeed equal to SSnu([b,0$k,a]) for 1<=b<=a<=K for a specific numeric k. Try:
#CheckSbka(4,10);
CheckSbka:=proc(k,K) local a,b:

evalb({seq(seq(SSnu([b,0$k,a])-Sbka(b,k,a),b=1..a),a=1..K)}={0});
end:



#IT(w): Given a word w of integers, finds the inversion table. Try:
#IT([4,3,2,1]);
IT:=proc(w) local gu,gu1,i,j:
gu:=[]:
for i from 1 to nops(w) do
 gu1:=0:
 for j from i+1 to nops(w) do
  if w[i]>w[j] then
    gu1:=gu1+1:
  fi:
od:
gu:=[op(gu),gu1]:
od:
gu:
end:



#IT0(w): Given a word w of integers, finds the Reduced inversion table. Try:
#IT0([4,3,2,1]);
IT0:=proc(w):Chop0(IT(pi)):end:

#NRD(pi): The number of reduced decompositions in S_n of the permutation pi (where n=nops(pi)). Try:
#NRD([5,4,3,2,1]);
NRD:=proc(pi):SSnu(Chop0(IT(pi))):end:



#CT(F,x): Given a rational function F in the list of variables x, finds the constant term/ Try:
#CT((1-x1/x2)/(x1^5*x2^5)/(1-x1-x2),[x1,x2]);
CT:=proc(F,x) local n,d,gu,x1 :
n:=nops(x):


d:=degree(denom(F),x[n]):

gu:=normal(coeff(series(F,x[n],d+2),x[n],0)):

if n=1 then
 RETURN(gu):
else
x1:=[op(1..n-1,x)]:
RETURN(CT(gu,x1)):
fi:
end:



#Ker2(x1,x2,z): The constant-termand for SbkaC(b,k,a) (q.v.). Try:
#Ker2(x1,x2,z)
Ker2:=proc(x1,x2,z)
-(-1+x2)/(x2*z^2+x1-z)
#normal((1/x1-x2/x1)/(1-z/x1+x2/x1*z^2)*(1-x2/x1)/(1-x1-x2)):
end:






#SbkaC(b,k,a): Computing SSnu([b,0$k,a]) via the constant term expression using The constant term expression. Try
#SbkaC(3,3,5);
SbkaC:=proc(b,k,a) local R,z,x1,x2,gu:

if not a>=b then
  RETURN(FAIL):
fi:

R:=(1-x2/x1)/(1-x1-x2)/x1^a/x2^b:

gu:=Ker2(x1,x2,z):
 gu:=normal(coeff(taylor(gu,z=0,k+1),z,k)):

CT(gu*R,[x1,x2]):

end:



#PermsDIV(n): all the permutations of n that have a decreasing reduced inversion table. Try:, together with it
#PermsDIV(5);
PermsDIV:=proc(n) local mu,pi,gu,hal:

mu:=permute(n):

gu:={}:

for pi in mu do
 hal:=Chop0(IT(pi)):
  if sort(hal,`>`)=hal then
   gu:=gu union {[pi,hal]}:
  fi:
od:

gu:

end:



#PermsCl(n): inputs a pos. integer n, and outputs a list of sets of length n-1 where the i-th entry is the set of permutations whose inversion table has i items. Try:
#PermsCl(5);
PermsCl:=proc(n) local T,gu,i,pi,hal:
gu:=permute(n):

for i from 1 to n-1 do
 T[i]:={}:
od:

for pi in gu do
 hal:=Chop0(IT(pi)):
 T[nops(hal)]:=T[nops(hal)] union {[pi,hal]}:
od:

[seq(T[i],i=1..n-1)]:


end:






#SSnuN(L): Given a list L of non-negative integers L, kicks out the 0-th if they exist, then sorts it in decreasing order, let's call the new list L1, and outputs
#SSnu(L)/YF(L1). Try
#SSnu([2,4,8]);
SSnuN:=proc(L) local L1, i:
L1:=[]:
for i from 1 to nops(L) do
 if L[i]<>0 then
  L1:=[op(L1),L[i]]:
fi:
od:
L1:=sort(L1,`>`):
SSnu(L)/YF(L1):
end:



#GuessPol1(L,n,d): inputs a list of pairs [[input,output], ...] a variable name n
#and a pos. integer d, outputs a polynomial of degree d f such that f(input)=output
GuessPol1:=proc(L,n,d) local a,i,f,var,eq:

if nops(L)<d+1 then
 print(`underdetermined make`, d, `smaller . `):
  RETURN(FAIL):
fi:
#f:=a[0]+a[1]*n+...+a[d]*n^d:
var:={ seq(a[i],i=0..d)}:
f:=add(a[i]*n^i,i=0..d):
# f(L[i][1])=L[i][2]

eq:={seq( subs(n=L[i][1], f)=L[i][2], i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

subs(var,f):

end:

#GuessPol(L,n): inputs a list of pairs [input,output] and tries to guess
#a polynomial of degree <=nops(L)-3
GuessPol:=proc(L,n) local f,d:

for d from 0 to nops(L)-3 do
 f:=GuessPol1([op(1..d+3,L)] ,n,d):
  if f<>FAIL then
   f:=GuessPol1(L,n,d):
     if f<>FAIL then
       RETURN(f):
     fi:
   fi:
od:

#print(`Too bad, it did not work out, generate more date, or bad news, it is not a polynomial`):
FAIL:

end:

#GuessRat1(L,n,d): inputs a list of pairs [[input,output], ...] a variable name n
#and a pos. integer d, outputs a rational of degree d f such that f(input)=output
GuessRat1:=proc(L,n,d) local a,b,i,f,var,eq:


if nops(L)<2*d+6 then
  RETURN(FAIL):
fi:

var:={ seq(a[i],i=0..d),  seq(b[i],i=1..d) }:
f:=add(a[i]*n^i,i=0..d)/(1+add(b[i]*n^i,i=1..d)):


eq:={seq( numer(subs(n=L[i][1], f)-L[i][2]), i=1..nops(L))}:



var:=solve(eq,var):


if var=NULL then
  RETURN(FAIL):
fi:

f:=subs(var,f):

if {seq( normal(subs(n=L[i][1], f)-L[i][2]), i=1..nops(L))}<>{0} then
 RETURN(FAIL):
fi:

factor(subs(var,f)):

end:



#GuessRat(L,n): inputs a list of pairs [[input,output], ...] a variable name n
#tries to guess a rational function in the variable n, for it. try:
#GuessRat([seq([i,(i+1)/(i+2)],i=1..10)]):
GuessRat:=proc(L,n) local d1,f,i,L1:

if GuessPol(L,n)<>FAIL then
 RETURN(GuessPol(L,n)):
fi:

for d1 from 1 to trunc((nops(L)-6)/2) do


L1:=[op(1..2*d1+6,L)]:


 f:=GuessRat1(L1,n,d1):
 if f<>FAIL then
if {seq( numer(subs(n=L[i][1], f)-L[i][2]), i=1..nops(L))}<>{0} then
 print(f, `did not work out`):
 RETURN(FAIL):
else
 RETURN(f):
fi:
fi:
od:
FAIL:
end:



#IncA(k,A): all tuples [a1,...,ak] such that 1<=a1<a2<...<ak. Try:
#IncA(3,5);
IncA:=proc(k,A) local gu,i,lu,lu1:
if k=1 then
 RETURN({ seq([i],i=1..A)}):
fi:
gu:={}:

for i from 1 to A do
 lu:=IncA(k-1,i-1):
 gu:=gu union {seq([op(lu1),i],lu1 in lu)}:
od:

gu:
end:

#IsGood(pi): Is the permutation pi good in the sense that SSnu([a[pi[1]],..., a[pi[k]]])=YF([a[k],..., a[1]]) for all a[1]<=a[2]<=..a[k]<=K. Try:
#IsGood([2,1,3],10);
IsGood:=proc(pi) local K,k,gu1,gu2,i:
k:=nops(pi):

evalb({seq(evalb(SSnu([seq(K*pi[i]+1,i=1..k)])=YF(sort([seq(K*pi[i]+1,i=1..k)],`>`))),K=3..4)}={true}):

end:


#GoodPerms(k): The set of permutations such that  for any increasing sequence of pos. integers of length k<=K v, SSnu(pi[v])=YF(Rev(v)). Try:
#GoodPerms(3);
GoodPerms:=proc(k) local lu,lu1,gu:
lu:=convert(permute(k),set):

gu:={}:

for lu1 in lu do
if IsGood(lu1) then
 gu:=gu union {lu1}:
fi:
od:
gu:
end:


#Check231D1(K): checks the explcit formula for SSnu([a1+1,a1,a3])/YF([a3,a1+1,a1])
#Check231D1(5,10);
Check231D1:=proc(K1,K2) local a1,a3:
evalb( {seq(seq(SSnu([a1+1,a1,a3])/YF([a3,a1+1,a1]) -(2*a3+a1+6)*(a3-a1+1)/2/(a3+3)/(a3-a1),a3=a1+2..a1+K2),a1=1..K1)}={0}):
end:



#Check231D2(K): checks the explcit formula for SSnu([a1+2,a1,a3])/YF([a3,a1+2,a1])
#Check231D2(5,10);
Check231D2:=proc(K1,K2) local a1,a3:
evalb( {seq(seq(SSnu([a1+2,a1,a3])/YF([a3,a1+2,a1]) -(a3-a1)*(3*a3^2+(18-a1)*a3+(-2*a1^2-7*a1+27))/3/(a3-a1+2)/(a3-a1-1)/(a3+3),a3=a1+3..a1+K2),a1=1..K1)}={0}):
end:


#Hafel2(ope,F,a,b,A,B): applies the operator ope (in terms of the shift operators A, B of a and b respectively)  to the expression F in a, b. Try:
#Hafel2(A+B,(a+b)!/a!/b!,a,b,A,B);

Hafel2:=proc(ope,F,a,b,A,B) local i,j,ope1,c,lu:

lu:=0:


for i from ldegree(ope,A) to degree(ope,A) do

 
 ope1:=coeff(ope,A,i):



  for j from ldegree(ope1,B) to degree(ope1,B) do
      c:=coeff(ope1,B,j):
      lu:=normal(lu+c*simplify(subs({a=a+i,b=b+j},F)/F)):
   od:
od:
factor(lu)*F:
end:

   



#Opek10(k,A,B): The operator such that SSnu([a1,0$k,a1])=Opek10(A,B)YF([a1,a1]). Try:
#Opek10(3,A,B):
Opek10:=proc(k,A,B)
option remember:
if k=0 then
 RETURN(1):
elif k=1 then
 RETURN(A):
else
expand(A*Opek(k-1,A,B)-Opek(k-2,A,B)*A/B):
fi:

end:



#SakaS(a,k): The explicit expression for SSnu([a,0^k,a]) 
#Try:
#SakaS(a,4);
SakaS:=proc(a,k) local b,A,B:
option remember:

subs(b=a,Hafel2(Opek10(k,A,B),YF([a,b]),a,b,A,B)):

end:


#Saka(a,k): SSnu([a,0$k,a]) according to the formula given by Saka(a,k)
#Try:
#Saka(5,2);
Saka:=proc(a,k) local A:


eval(subs({A=a},SakaS(A,k) )):
end:




#SakaC(a,k): Computing SSnu([a,0$k,a]) via the constant term expression using The constant term expression. Try
#SakaC(5,4);
SakaC:=proc(a,k) local R,z,x1,x2,gu:


R:=(1-x2/x1)/(1-x1-x2)/x1^a/x2^a:

gu:=Ker2(x1,x2,z):
gu:=normal((z/x1-x2/x1*z^2)*gu):

 gu:=normal(coeff(taylor(gu,z=0,k+1),z,k)):

CT(gu*R,[x1,x2]):

end:

#Opek11(k,A,B): The operator such that SSnu([a1,0$k,a1-1])=Opek11(k,A,B)YF([a1,a1-1]). Try:
#Opek11(3,A,B):
Opek11:=proc(k,A,B)
option remember:
if k=0 then
 RETURN(1):
elif k=1 then
 RETURN(1):
else
expand(B*Opek10(k-1,A,B)-Opek(k-2,A,B)):
fi:

end:


#Opek12(k,A,B): The operator such that SSnu([a1,0$k,a1-2])=Opek12(k,A,B)YF([a1,a1-2]). Try:
#Opek12(3,A,B):
Opek12:=proc(k,A,B)
option remember:
if k<=2 then
 RETURN(1):
else
expand(B*Opek11(k-1,A,B)-B/A*Opek10(k-2,A,B)):
fi:

end:



#Saka1S(a,k): The explicit expression for SSnu([a,0^k,a-1]) 
#Try:
#Saka1S(a,4);
Saka1S:=proc(a,k) local b,A,B:
option remember:

subs(b=a-1,Hafel2(Opek11(k,A,B),YF([a,b]),a,b,A,B)):

end:



#Saka1(a,k): SSnu([a,0$k,a-1]) according to the formula given by Saka1S(a,k)
#Try:
#Saka1(5,2);
Saka1:=proc(a,k) local A:


eval(subs({A=a},Saka1S(A,k) )):
end:




#Saka2S(a,k): The explicit expression for SSnu([a,0^k,a-2]) 
#Try:
#Saka2S(a,4);
Saka2S:=proc(a,k) local b,A,B:
option remember:

subs(b=a-2,Hafel2(Opek12(k,A,B),YF([a,b]),a,b,A,B)):

end:



#Saka2(a,k): SSnu([a,0$k,a-2]) according to the formula given by Saka2S(a,k)
#Try:
#Saka2(5,2);
Saka2:=proc(a,k) local A:


eval(subs({A=a},Saka2S(A,k) )):
end:


#Opek1j(j,k,A,B): The operator such that SSnu([a1,0$k,a1-j])=Opek1j(j,k,A,B)YF([a1,a1-j]). Try:
#Opek1j(2,3,A,B):
Opek1j:=proc(j,k,A,B)
option remember:
if k<=j then
 RETURN(1):
fi:



if j=0 then
if k=1 then
 RETURN(A):
else
 RETURN(expand(A*Opek(k-1,A,B)-Opek(k-2,A,B)*A/B)):
fi:
elif j=1 then
RETURN(expand(B*Opek1j(0,k-1,A,B)-Opek(k-2,A,B))):
else
RETURN(expand(B*Opek1j(j-1,k-1,A,B)-B/A*Opek1j(j-2,k-2,A,B))):
fi:

end:




#SakajS(j,a,k): The explicit expression for SSnu([a,0^k,a-j]) 
#Try:
#SakajS(5,a,3);
SakajS:=proc(j,a,k) local b,A,B:
option remember:

simplify(factor(subs(b=a-j,Hafel2(Opek1j(j,k,A,B),YF([a,b]),a,b,A,B)))):

end:



#Sakaj(j,a,k): SSnu([a,0$k,a-j]) according to the formula given by SakajS(j,a,k)
#Try:
#Sakaj(3,10,5);
Sakaj:=proc(j,a,k) local A:


eval(subs({A=a},SakajS(j,A,k) )):
end:




#GuessRatM(DS,x): Given a  a data set DS, and a list of variables x, tries to fit a rational function to it. Try:
#GuessRatM({seq(seq([[i1,i2],(i1+i2)/(i1+i2+1),i2=1..i1),i1=10..20)},[x,y]):
GuessRatM:=proc(DS,x) local k,A,DS1, Rishon,a1,hal,gu:
k:=nops(x):

if k=1 then
 DS1:={seq([A[1][1],A[2]], A in DS)}:

 RETURN(GuessRat(DS1,x[1])):
fi:

Rishon:={seq(A[1][1],A in DS)}:

for a1 in Rishon do
 DS1[a1]:={}:
od:

for A in DS do
a1:=A[1][1]:
DS1[a1]:=DS1[a1] union {[[op(2..k,A[1])],A[2]]}:
od:

gu:=[]:



for a1 in Rishon do


hal:=GuessRatM(DS1[a1],[op(2..k,x)]):



if hal=FAIL then
  RETURN(FAIL):
else
gu:=[op(gu),[a1,hal]]:
fi:
od:


GuessRat(gu,x[1]):

end:


#GuessPolM(DS,x): Given a  a data set DS, and a list of variables x, tries to fit a polynomial function to it. Try:
#GuessPolM({seq(seq([[i1,i2],(i1+i2)/(i1+i2+1),i2=1..i1),i1=10..20)},[x,y]):
GuessPolM:=proc(DS,x) local k,A,DS1, Rishon,a1,hal,gu:
k:=nops(x):

if k=1 then
 DS1:={seq([A[1][1],A[2]], A in DS)}:

 RETURN(GuessPol(DS1,x[1])):
fi:

Rishon:={seq(A[1][1],A in DS)}:

for a1 in Rishon do
 DS1[a1]:={}:
od:

for A in DS do
a1:=A[1][1]:
DS1[a1]:=DS1[a1] union {[[op(2..k,A[1])],A[2]]}:
od:

gu:=[]:



for a1 in Rishon do


hal:=GuessPolM(DS1[a1],[op(2..k,x)]):



if hal=FAIL then
  RETURN(FAIL):
else
gu:=[op(gu),[a1,hal]]:
fi:
od:


GuessPol(gu,x[1]):

end:


#K:=proc(a2,a3,j):  if not a2-a3>=2 then RETURN(FAIL): fi:SSnu([a2,a3,a2+j])/YF([a2,a2,a3])/(j+2)/(2*a2+a3+j)!*(j+a2+a3)!: end:
#K:=proc(a2,a3,j):  if not a2-a3>=2 then RETURN(FAIL): fi:SSnu([a2,a3,a2+j])/YF([a2,a2,a3])/(j+2)/(2*a2+a3+j)!: end:
K:=proc(a2,a3,j):  if not a2-a3>=2 then RETURN(FAIL): fi:SSnu([a2,a3,a2+j])/(YF([a2,a2,a3])*(j+2)*(2*a2+a3+j)!/(j+a2+3)!): end:



Khal:=proc(a2,a3,j) local gu,lu: 
gu:=((a2-a3+1)*j^2+ (a3^2+(-4*a2-6)*a3+3*a2^2+7*a2+4)*j+(2*a2^3-4*a2^2*a3+2*a2*a3^2+9*a2^2-12*a2*a3+3*a3^2+10*a2-9*a3+3)):
gu:=gu*(a2+2)!/(2*a2+a3)!/(a2-a3+1)/(a2-a3+2):
K(a2,a3,j)/gu:
end:



#S231o(a2,a3,j): The explicit formula for SSnu([a2,a3,a2+j]) if a2-a3>=2. Try:
#S231o(20,30,3);
S231o:=proc(a2,a3,j) local gu,POL: 
POL:=((a2-a3+1)*j^2+ (a3^2+(-4*a2-6)*a3+3*a2^2+7*a2+4)*j+(2*a2^3-4*a2^2*a3+2*a2*a3^2+9*a2^2-12*a2*a3+3*a3^2+10*a2-9*a3+3)):
gu:=(a2+2)!/(2*a2+a3)!/(a2-a3+1)/(a2-a3+2):
gu:=gu*YF([a2,a2,a3])*(j+2)*(2*a2+a3+j)!/(j+a2+3)!:
gu:=simplify(gu):
[POL,gu]:
end:


#S231(a2,a3,a1): The explicit formula for SSnu([a2,a3,a1]) if a1>=a2>=a3. Try:
#S231(a2,a3,a1);
#S231(20,30,3);
S231:=proc(a2,a3,a1) local gu,POL,j: 

POL:=((a2-a3+1)*j^2+ (a3^2+(-4*a2-6)*a3+3*a2^2+7*a2+4)*j+(2*a2^3-4*a2^2*a3+2*a2*a3^2+9*a2^2-12*a2*a3+3*a3^2+10*a2-9*a3+3)):
POL:=factor(expand(subs(j=a1-a2,POL))):
gu:=(a2+2)!/(2*a2+a3)!/(a2-a3+1)/(a2-a3+2):
gu:=gu*YF([a2,a2,a3])*(j+2)*(2*a2+a3+j)!/(j+a2+3)!:
gu:=subs(j=a1-a2,gu):
gu:=simplify(gu):
POL*gu:
end:


K1:=proc(a2,a3,n,N) local j: findrec([seq(K(a2,a3,j),j=1..30)],3,1,n,N);end:

K2:=proc(a2,a3,n) local N,j: expand(denom(coeff(findrec([seq(K(a2,a3,j),j=1..30)],5,1,n,N),N,0)));end:


#Coe1(a2,a3): The conjectured polynomial in a2 and a3 of the coeff. of n in K2(a2,a3,n) (normalizing that the coefficient of n^2 is (a2-a3+1). K is a guessing paramter. Try:
#Coe1(a2,a3);
Coe1:=proc(a2,a3) local DS,i2,i3,lu,n,hal,H,DS1:

DS1:=[]:

for i3 from 2 to 11 do
DS:={}:
 for i2 from i3+2 to  i3+11 do

  lu:=K2(i2,i3,n):

  if coeff(lu,n,2)=i2-i3+1 then
   DS:=DS union {[i2,coeff(lu,n,1)]}:
  fi:
od:

  hal:=GuessPol(DS,a2):
 if hal<>FAIL then
  DS1:=[op(DS1),[i3,hal]]:
fi:
od:

GuessPol(DS1,a3):
end:
#a3^2+(-4*a2-6)*a3+3*a2^2+7*a2+4


#Coe0(a2,a3): The conjectured polynomial in a2 and a3 of the coeff. of n in K2(a2,a3,n) (normalizing that the coefficient of n^2 is (a2-a3+1). K is a guessing paramter. Try:
#Coe0(a2,a3);
Coe0:=proc(a2,a3) local DS,i2,i3,lu,n,hal,DS1:

DS1:=[]:

for i3 from 3 by 2 to 15 do
print(`i3 is`, i3):
DS:={}:
 for i2 from i3+2 while nops(DS)<6 do  

  lu:=K2(i2,i3,n):

  if coeff(lu,n,2)=i2-i3+1 then
   DS:=DS union {[i2,coeff(lu,n,0)]}:
  fi:
od:

  hal:=GuessPol(DS,a2):
print(`hal is`, hal):
 if hal<>FAIL then
  DS1:=[op(DS1),[i3,hal]]:
fi:
od:

lprint(DS1):

GuessPol(DS1,a3):
end:
#2*a2^3-4*a2^2*a3+2*a2*a3^2+9*a2^2-12*a2*a3+3*a3^2+10*a2-9*a3+3




#SSnuCF3(L): The closed form expression for SSnu(L) for any list of positive integers [a1,a2,a3]. Try:
#{seq(seq(seq(SSnuCF3([a1,a2,a3])-SSnu([a1,a2,a3]),a1=1..10),a2=1..10),a3=1..10)};
SSnuCF3:=proc(L)

if nops(L)<>3 then
 print(L, `should be of length 3`):
 RETURN(FAIL):
fi:

if L[3]>=L[1] and L[1]>=L[2] then
 S231(op(L)):
else
YF(sort(L,`>`)):
fi:
end:





#SbkaStory(K): Explicit expressions in a and b, for the Number of Ways of Playing Stanley Solitaire starting with position [b,0$k,a] where a>b, for k from 0 to K. Do
#SbkaStory;
SbkaStory:=proc(K) local a,b,k,gu,t0:
t0:=time():
print(`Explicit Expressions for The Number of ways of playing Stanley Solitaire staring with position [b,0^k,a] for all a>b and k from 0 to `, K):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

for k from 0 to K do
gu:=SbkaS(b,k,a):
print(`Theorem Number`, k, `If b<a, then the number of ways of playing Stanley Solitaire starting with position`,[b,0$k,a], ` is `):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`For example if a=30 and b=20 then the number is`, eval(subs({a=30,b=20},gu))):
print(``):
print(`This agrees with the value obtained by numerical dynamical programming that is`, SSnu([20,0$k,30])):
print(``):
od:
print(`---------------------------`):
print(``):
print(`This concludes this paper that took`, time()-t0, `seconds to produce `):
print(``):
end:



#SakajStory(J,K): Explicit expressions in a, for the Number of Ways of Playing Stanley Solitaire starting with position [a,0$k,a-j] where  for k from j from 0 to K and
#k from j+1 to J+K. Try
#SakajStory(6,6);
SakajStory:=proc(J,K) local a,j,k,gu,t0:
t0:=time():
print(`Explicit Expressions for The Number of ways of playing Stanley Solitaire staring with position [a,0^k,a-j] for all a, and j from 0 to `, J. `and k from j+1 to`, j+K):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

for j from 0 to J do
for k from j+1 to j+K do
gu:=SakajS(j,a,k) :
print(`Theorem Number`, [j,k-j], ` The number of ways of playing Stanley Solitaire starting with position`,[a,0$k,a-j], ` is `):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`For example if a=`, 50+J, ` then the number is`, eval(simplify(subs({a=50+J},gu)))):
print(``):
print(`This agrees with the value obtained by numerical dynamical programming that is`, SSnu([50+J,0$k,50+J-j])):
print(``):
od:
od:

print(`---------------------------`):
print(``):
print(`This concludes this paper that took`, time()-t0, `seconds to produce `):
print(``):
end:





#ITr(pi): Like IT(w) but done recursively. Try:
#ITr([4,3,2,1]);
ITr:=proc(pi) local a,pi1,n,i1:
option remeber:

n:=nops(pi):

if n=0 then
 RETURN([]):
fi:

a:=pi[1]:
pi1:=[op(2..n,pi)]:
pi1:=subs({seq(i1=i1-1,i1=a+1..n)},pi1):
[a-1,op(ITr(pi1))]:

end:


iITr:=proc(w) local n,f,w1,i1,pi1:
n:=nops(w):


if n=0 then
 RETURN([]):
fi:


if not {seq(evalb(w[i1]<=n-i1),i1=1..n)}={true} then
 RETURN(FAIL):
fi:


f:=w[1]+1:

w1:=[op(2..n,w)]:

pi1:=iITr(w1):

pi1:=subs({seq(i1=i1+1,i1=f..n-1)},pi1):
[f,op(pi1)]:
end:



#PosToPer(L): Given a starting position in Stanley Solitaire,L, finds the permutation of smallest length such that its inversion table is L followed by 0's. Try:
#PosToPer([20,10,30]);
PosToPer:=proc(L) local i,n,L1:

n:=max(seq(L[i]+i+1,i=1..nops(L))):

L1:=[op(L),0$(n-nops(L))]:

iITr(L1):

end:


#Mu(pi): Given a permutation pi, finds the quantity mu(pi) defined on top of p. 367 of Richard Stanley's article "On the Number of Reduced Decompositions of Elements in the Coxeter Group"
#Eur. J. Comb. 5 (1984). 359-372
Mu:=proc(pi) Chop0(sort(IT(pi), `>`)): end:





#IT1(w): Given a word w of integers, finds the Reverse inversion table. Try:
#IT1([4,3,2,1]);
IT1:=proc(w) local gu,gu1,i,j:
gu:=[]:
for i from 1 to nops(w) do
 gu1:=0:
 for j from 1 to i-1 do
  if w[j]>w[i] then
    gu1:=gu1+1:
  fi:
od:
gu:=[op(gu),gu1]:
od:
gu:
end:



#Lam(pi): Given a permutation pi, finds the quantity lambda(pi) defined on top of p. 367 of Richard Stanley's article "On the Number of Reduced Decompositions of Elements in the Coxeter Group"
#Eur. J. Comb. 5 (1984). 359-372
Lam:=proc(pi) Conj(Chop0(sort(IT1(pi), `>`))): end:


#IsGoodS(pi): Is the permutation pi good in Stanley's sense? Try
#IsGoodS([2,1,3]);
IsGoodS:=proc(pi) local Pi1,i,K:

for K from nops(pi)+2 to nops(pi)+4 do
Pi1:=PosToPer([seq(pi[i]*K+1,i=1..nops(pi))]):

if Mu(Pi1)<>Lam(Pi1) then
 RETURN(false):
fi:
od:
true:
end:



#GoodPermsS(k,K): The set of permutations in Stanley's sense. Try:
#GoodPermsS(3,5);
GoodPermsS:=proc(k) local lu,lu1,gu:
lu:=convert(permute(k),set):

gu:={}:

for lu1 in lu do
if IsGoodS(lu1) then
 gu:=gu union {lu1}:
fi:
od:
gu:
end:



#IsGoodSS(pi): Is the permutation pi  such that Lam(pi)=Mu(pi)?
#IsGoodSS([2,1,3]);
IsGoodSS:=proc(pi) 
evalb(Mu(pi)=Lam(pi)):
end:




#GoodPermsSS(k): The set of permutations of length k, such that Lam(pi)=Mu(pi)
#GoodPermsSS(3);
GoodPermsSS:=proc(k) local lu,lu1,gu:
lu:=convert(permute(k),set):

gu:={}:

for lu1 in lu do
if IsGoodSS(lu1) then
 gu:=gu union {lu1}:
fi:
od:
gu:
end:



#S231Story(): A paper about playing Stanley Solitaire with starting position [a2,a3,a1] where a1>=a2>=a3. Try:
#S231SStory();
S231Story:=proc() local a1,a2,a3,gu,pi,gu1:

print(`An Explicit Formula for the Number of Ways to Play Stanley Solitaire with Starting position [a2,a3,a1] where a1>=a2>=a3`):
print(``):
pring(`By Shalosh B. Ekhad`):
print(``):
gu:=S231(a2,a3,a1):

print(`Theorem: If a1>=a2>=a3>=0 then the number of ways of Playing Stanley Solitaire with initial position`):
print(``):
print(` [a2,a3,a1] `):
print(``):
print(`is equal to`):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(``):
print(`Cor.: Let pi be a permutation of length n=a1+4, with `):
print(`pi[1]=a[2]+1,pi[2]=a[3]+1,pi[3]=a[1]+3`):
print(``):
print(`and pi[4]>=...>=pi[n]`):
print(``):
print(`the number of reduced decomposition of pi is, as above i.e.`):
print(``):
print(gu):
print(``):
print(`For example, taking a1=300,a2=200,a3=100`):
pi:=PosToPer([200,100,300]):
print(`the permutation `):
print(pi):
print(``):
gu1:=eval(subs({a1=300,a2=200,a3=100},gu));
print(``):
print(`has `, gu1, `reduced decompositions`):
print(``):


end: