print("This Maple package computes the number of 123-avoiding ordered set partitions"):
print("To find the number of 123-avoiding partitions of {1,...,n} into k blocks, try Akblocks(n,k)"):
print("To find the number of 123-avoiding partitions where block 1 has b1 elements, block 2 has b2 elements, ..., "):
print("and block k has k elements, try A([b1,b2,...,bk])"):
print(""):
print("seq(A([2$k]),k=1..n) produces the first n terms of OEIS sequence A220097"):
print("seq(Akblocks(k+1,k),k=1..n) produces the first n terms of OEIS sequence A220101"):

with(combinat):

A:=proc(B) local n; 
n:=convert(B, `+`):
if nops(B) = 1 then return 1:
  elif nops(B) = 2 then return binomial(n, B[1]) 
  elif evalb({op(B)} = {1}) then return binomial(2*n, n)/(n+1):
  else return add(C(B, s), s in choose({$1..n}, B[1])):
fi:
end:

C:=proc(B, s) local n,tot,ab,r: option remember:
n:=convert(B, `+`):
if nops(B) = 1 and s={$1..B[1]} then return 1:
elif nops(B)=1 and s<>{$1..B[1]} then return 0:
elif nops(B)=2 then return 1:
else 

tot:=0:
for r in choose({$1..n} minus s, B[2]) do
ab:=Findab([s,r]):
if ab=[0,0] then
tot:=tot + E(B,[s,r],0,0):
else
tot:=tot+E(B,[s,r],ab[1],ab[2]):
fi:
od:
return tot:
fi:
end:


Findab:=proc(L) local s, r, c, a, b, k, i, j: option remember:
#returns a and b if prefix arrangement contains
#a 12 pattern.  returns [0,0] if pattern is decreasing
s:=L[1]: 
s:=sort([op(s)]):
r:=L[2]: 

k:=0: 
for i to nops(r) do 
c:=0: 
for j to nops(s) do 
if s[j] < r[i] then c:=c+1: fi: 
od:
if 0 < c and k = 0 then k:=r[i] fi: 
od:

a:=0: 
b:=0:
for i to nops(r) do 
if r[i] < k then b:=b+1: fi:
od:

for i to nops(s) do
if s[i] < k then a:=a+1: fi:
od:

[a, b]:
end:

E:=proc(B,sr,a,b) local n,Big,i,s,r: option remember: 
n:=convert(B,`+`): 
s:=sort([op(sr[1])]):  
r:=sort([op(sr[2])]):  

if a=0 and b=0 then  
return C([seq(B[i],i=2..nops(B))],sr[2]):  
#this was the case where the 2 block prefix avoids the pattern 12

else

Big:={}:
for i from a+1 to nops(s) do
Big:=Big union {s[i]}:
od:
for i from b+2 to nops(r) do
Big:=Big union {r[i]}:
od:

#Big is all elements in the first block that are at least
#as large as first 12 pattern

#print(B,sr,a,b,Big,n,r[b+1],n-r[b+1]):

if n-r[b+1]<>nops(Big) then return 0:
else
if b=0 then return C([a,seq(B[i],i=3..nops(B))],{op(s)} minus Big):
else return C([b,seq(B[i],i=3..nops(B))],{op(r)} minus Big minus {r[b+1]}):
fi:

fi:

fi:


end:

part:=proc(n,k) local S,i,S2,s: option remember:
#integer partitions of [n] into
#k parts -- used for block sizes into k blocks
if n<k then return {}:
elif n=k then return {[1$k]}:
elif k=1 then return {[n]}:
else
S:={}:
for i from 1 to n-k+1 do
S2:=part(n-i,k-1):
for s in S2 do
S:=S union {[i,op(s)]}:
od:
od:
fi:
S:
end:

Akblocks:=proc(n,k) local P:
P:=part(n,k):

add(A(p),p in P):
end: