with(combinat):

Help:=proc():
print(`FG(L),IsCon(L1, L2),Pnkm(n,k,m),Pnm(n,k,m),YLattice(M,N)`):
print(`GetS3n(n),GetE3n(n),Decomp3n(n)`):
print(`DEPRECATED: PartsInL(L),YoungLattice(M,N),YGrabChain(L),YAigner(n)`):
end:

FG := proc(L) local i, j: {seq(seq([i, j], j = 1 .. L[i]), i = 1 .. nops(L))}: end:

IsCon := proc(L1, L2) evalb(FG(L1) minus FG(L2) = {}): end:

ConjPart:=proc(L,M,N): [seq(N-L[M-i],i=0..M-1)]: end:

#Pnkm(n,k,m): The LIST of integer partitions of n with largest part k
#with at most m parts
Pnkm:=proc(n,k,m) local k1,L,L1:
option remember:

if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN([]):
fi:

if k>n or m=0 or k*m<n then
 RETURN([]):
fi:

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

L:=[]:

for k1 from 1 to min(n-k,k) do
 L1:=Pnkm(n-k,k1,m-1):
 L:=[op(L), seq([k, op(L1[j])],j=1..nops(L1))]:
od:
#print(n,k,m,L):
L:

end:

#Pnm(n,k,m): The list of integer partitions of n in lex. order with at most m parts
#and no part larger than k
Pnm:=proc(n,k,m) local i:
option remember:
[seq(op(Pnkm(n,i,m)),i=1..k)]:
end:

#YLattice(M,N): Builds the Young lattice for the integer partition with at most
#M parts and largest part at most N
YLattice:=proc(M,N) local S,L,s,n,k:
n:=N*M:
L:=[[[]]]:
for s from 1 to n do
  L:=[op(L),Pnm(s,N,M)]:
od:
L:
end:

#YGrabChain(L): inputs the still-available lattice
#outputs a chain (taken greedily) and the smaller
#L with the members of C removed
YGrabChain:=proc(L) local i,C,L1,L2,st,c,grabflag:
#print(`GC`,L):
L1:=[]:
for i from 1 to nops(L) while L[i]=[] do L1:=[op(L1),[]]: od:

if i=nops(L)+1 then RETURN(FAIL): fi:
st:=i:
C:=[L[st][1]]:
L1:=[op(L1),[seq(L[st][i],i=2..nops(L[st]))]]:
for i from st+1 to nops(L) do
  grabflag:=true:
  L2:=[]:
#print(i):
  for c in L[i] do
#print(C,C[i-st],c,L[i]):
    if grabflag and IsCon(C[i-st],c) then
#   if grabflag and nops(c)>=nops(C[i-st]) and {true,seq(evalb(C[i-st][j]<=c[j]),j=1..nops(C[i-st]))}={true} then
#   if grabflag and (convert(C[i-st],`set`) subset convert(c,`set`)) then

#print(c):
      C:=[op(C),c]:
      grabflag:=false:
    else
      L2:=[op(L2),c]:
    fi:
  od:
  L1:=[op(L1),L2]:
od:
[C,L1]:

end:

#YAigner(n): Inputs a pos. integer n, and outputs a SCD
#of the Boolean lattice of subsets of [n]
#following the great combinatorialist Martin Aigner's method
#of lexicographic greed
YAigner:=proc(M,N) local CD,L,k,C:
CD:=[]:
L:=YLattice(M,N):
while L<>[[]$(M*N+1)] do
#print(L,CD):
  C:=YGrabChain(L):
  CD:=[op(CD),C[1]]:
  L:=C[2]:
od:
CD:
end:

GetS3n:=proc(n) local L,k,l:
option remember:
L:=[]:
for k from 0 to floor(n/4) do
  L:=[op(L),[4*k,2*k,0]]:
od:
for l from 2 to n do
  for k from 0 to floor((n-l)/4) do
    L:=[op(L),[4*k+l,2*k,0]]:
  od:
od:
L:
end:

GetE3n:=proc(n) local L,k,l:
option remember:
L:=[]:
for k from 0 to floor(n/4) do
  L:=[op(L),ConjPart([4*k,2*k,0],3,n)]:
od:
for l from 2 to n do
  for k from 0 to floor((n-l)/4) do
    L:=[op(L),ConjPart([4*k+l,2*k,0],3,n)]:
  od:
od:
L:
end:

#Decomp3n(n): Implements the Sperner decomposition of L(3,n) using the algorithm
#in chapter 3 of Xin & Zhong's paper of 22April2021 on the Arxiv
Decomp3n:=proc(n) local S3n,E3n,lambda,CD,C,anchor:
S3n:=GetS3n(n):
E3n:=[seq(ConjPart(lambda,3,n),lambda in S3n)]:
CD:=[]:
for anchor in S3n do
  lambda:=anchor:
  C:=[anchor]:
  while not lambda in E3n do
    if lambda in GetE3n(lambda[1]) then
      lambda:=[lambda[1]+1,lambda[2],lambda[3]]:
    elif (lambda[2]+lambda[3]) mod 2 = 0 then
      lambda:=[lambda[1],lambda[2]+1,lambda[3]]:
    elif [lambda[1]-1,lambda[2],lambda[3]+1] in GetE3n(lambda[1]-1) then
      lambda:=[lambda[1],lambda[2]+1,lambda[3]]:
    else
      lambda:=[lambda[1],lambda[2],lambda[3]+1]:
    fi:
    C:=[op(C),lambda]:
  od:
  CD:=[op(CD),C]:
od:
CD:
end:


################## Early Attempts ##########################

#PartsInL(L): Returns the set of integer partitions that fit inside L
PartsInL:=proc(L) local S,L1,i,n:
if L=[] then RETURN({[]}): fi:
S:={L}:
n:=nops(L):
for i from 1 to n-1 do
  if L[i]>L[i+1] then
    L1:=[op(1..i-1,L),L[i]-1,op(i+1..n,L)]:
    S:=S union PartsInL(L1):
  fi:
od:
if L[n]=1 then
  S:=S union PartsInL([op(1..n-1,L)]):
else
  S:=S union PartsInL([op(1..n-1,L),L[n]-1]):
fi:
S:
end:



#YoungLattice(M,N): Builds the Young lattice for the integer partition with at most
#M parts and largest part at most N
YoungLattice:=proc(M,N) local S,L,s,n,k:
n:=N*M:
S:=PartsInL([N$M]):
L:=[[]$(M*N+1)]:
for s in S do
  k:=add(s):
  L[k+1]:=[op(L[k+1]),s]:
od:
L:
end:

################### OLD STUFF ###############################
# 1.

#Aigner(n): Inputs a pos. integer n, and outputs a SCD
#of the Boolean lattice of subsets of [n]
#following the great combinatorialist Martin Aigner's method
#of lexicographic greed
Aigner:=proc(n) local CD,L,k,C:
CD:=[]:
L:=[seq(choose(n,k),k=0..n)]:
while L<>[[]$(n+1)] do
  C:=GrabChain(L):
  CD:=[op(CD),C[1]]:
  L:=C[2]:
od:
CD:
end:

#GrabChain(L): inputs the still-available lattice
#outputs a chain (taken greedily) and the smaller
#L with the members of C removed
GrabChain:=proc(L) local i,C,L1,L2,st,c,grabflag:
L1:=[]:
for i from 1 to nops(L) while L[i]=[] do L1:=[op(L1),[]]: od:

if i=nops(L)+1 then RETURN(FAIL): fi:
st:=i:
C:=[L[st][1]]:
L1:=[op(L1),[seq(L[st][i],i=2..nops(L[st]))]]:
for i from st+1 to nops(L) do
  grabflag:=true:
  L2:=[]:
  for c in L[i] do
    if grabflag and (convert(C[i-st],`set`) subset convert(c,`set`)) then
      C:=[op(C),c]:
      grabflag:=false:
    else
      L2:=[op(L2),c]:
    fi:
  od:
  L1:=[op(L1),L2]:
od:
[C,L1]:

end:

# 2.

#SCD(n): uses the idea at the very end of the lecture to construct from each
#chain of SCD(n-1) by creating two new chains of SCD(n).
SCD:=proc(n) local C,P,c,i:
option remember:
if n<=0 then RETURN([[[]]]): fi:
P:=SCD(n-1):
C:=[]:
for c in P do
#print(c):
  if nops(c)>1 then
    C:=[op(C),[op(c),[op(c[-1]),n]],
        [seq([op(c[i]),n],i=1..nops(c)-1)]]:
  else
    C:=[op(C),[op(c),[op(c[-1]),n]]]:
  fi:
od:
C:
end:


#Verify that you get the same output as Aigner(n) #### Verified up to n=13 ####
#Making SCD output in the same order as Aigner proved too annoying, so I convert to sets
#in order to make sure that the actual decompositions match.
#CheckSCDs(n)
CheckSCDs:=proc(n) local i:
seq(evalb(convert(Aigner(i), set) = convert(SCD(i), set)), i = 0 .. n):
end: