Help20:=proc(): print(` SnkSeq(N,k), SnkSeqC(N,k), ReducedCycle(C), PtoCfull(C),  cnk(n,k) , cnkSeq(N,k), cnkSeqC(N,k) `): end:

with(combinat):



#Snk(n,k): The Stirling numbers of the second-kind, aka the NUMBER of set-partitions of an n-element set
Snk:=proc(n,k)  option remember:
if n=1 then
 if k=1 then
 RETURN(1):
 else
 RETURN(0):
 fi:
fi:

Snk(n-1,k-1)+k*Snk(n-1,k):
end:

#SnkSeq(N,k): The first N terms of the the sequence Snk(n,k) n=1..N, for a fixed k,done by the method of Lecture 11. Try:
#SnkSeq(30,k)
SnkSeq:=proc(N,k) local n1:
[seq(Snk(n1,k),n1=1..N)]:
end:

#SnkSeqC(N,k): The first N terms of the the sequence Snk(n,k) n=1..N, for a fixed k,done by the method of Lecture 11. Try:
#SnkSeqC(30,k)
SnkSeqC:=proc(N,k) local f,x,i:
f:=(exp(x)-1)^k/k!:

f:=taylor(f,x=0,N+1):

[seq(coeff(f,x,i)*i!,i=1..N)]:

end:


#ReducedCycle(C): Given a cycle with arbitrary labels, outputs its REDUCED form, followed by the set of labels
ReducedCycle:=proc(C) local S,C1:
S:=sort(C):
C1:=subs({seq(S[i]=i,i=1..nops(S))},C):
[C1,convert(S,set)]:
end:


PtoCfull:=proc(P) local C,S:
C:=PtoC(P):
[seq(ReducedCycle(C[i]),i=1..nops(C))]:
end:

#cnk(n,k): the number of permutations with exactly k cycles. Using the  recurrence
cnk:=proc(n,k) option remember:
if n=1 then
 if k=1 then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:
cnk(n-1,k-1)+(n-1)*cnk(n-1,k):
end:


#cnkSeq(N,k): The first N terms of the sequence cnk(n,k) for a fixed k. Try:
#cnkSeq(20,3);
cnkSeq:=proc(N,k) local n:
[seq(cnk(n,k),n=1..N)]:
end:


cnkSeqC:=proc(N,k) local f,x,i:
f:=(-log(1-x))^k/k!:

f:=taylor(f,x=0,N+1):

[seq(coeff(f,x,i)*i!,i=1..N)]:
end:


##Start Maple Code from M11.txt
#Maple Code for Lecture 11 of Combinatorics , Math 454, taught by Dr. Z.
Help11:=proc(): print(` SPk(n,k), SP(n) , Snk(n,k), Bell1(n),  Bell2(n), Bell3(n) `):
print(`SeqBell0(N),SeqBell1(N),SeqBell2(N),SeqBell3(N)`):
end:

#SPk(n,k): The set of Set-Partiitions of {1,...,n} into k parts. In other words the set of
#sets {S1,S2,..., Sk} such that the Si are pairwise disjoint and their union is {1, ..., n}
#For example, SPk(3,2) is { {{1},{2,3}}, {{2},{1,3}}, {{3},{2,3}}}
#NOTE THAT THE OUTPUT IS A SET OF SETS OF SETS (since a set-partion is already a set-of-sets so
#the set of all set-partiotions satisfying whatever coniditions is the set of sets of sets
SPk:=proc(n,k) local SSP,SSP1,i,ssp1:
option remember:


#THIS IS THE INITIAL CONDITION
if n=1 then
 if k=1 then
  RETURN({{{1}}}):
 else
  RETURN({}):
fi:
fi:


#We now construct all the set-partitions of {1, ...,n} by looking at the "social life" of n
#If it is all-alone, then removing {n} is a typical set partition of {1,..,n-1} INTO k-1 subsets
#Otherwise n may belong to one of the k members of the set-partition and removing it leads
#to a set partition of {1,.., n-1} with k subsets, and going back we can stick n into each of
#these

#Let's call the final output SSP, we will dynamically enlarge it starting with  the emtpy set

SSP:={}:

#We start with the case that n is in its own component we recursively call this procedure with (n-1,k-1)
SSP1:=SPk(n-1,k-1):

SSP:={seq( ssp1 union {{n}}, ssp1 in SSP1)}:

#This takes care of the first case
#We now do the case where removing n does not change the number of components
SSP1:=SPk(n-1,k):



for ssp1 in SSP1 do
#We look at all the members of SPk(n-1,k)
for i from 1 to k do
 SSP:=SSP union {{op(1..i-1,ssp1), ssp1[i] union {n}, op(i+1..k,ssp1)}}:
od:
od:
SSP:
end:




#SP(n): The set of ALL set-partitions of {1,...,n} regardless of the number of constituents
SP:=proc(n) local k:{seq(op(SPk(n,k)),k=1..n)}:end:




#Bell1(n): The number  of set partitions (i.e. nops(SP(n)) done the FIRST way, adding-up Snk(n,k) from k=1 to k=n
Bell1:=proc(n) local k: add(Snk(n,k),k=1..n): end:

#Bell2(n): The number  of set partitions (i.e. nops(SP(n)) done the SECOND way, using the recurrence proved in the lecture
Bell2:=proc(n) local k: option remember:
if n=0 then
 1:
else
add(binomial(n-1,k)*Bell2(n-1-k),k=0..n-1):
fi:
end:

#Bell3(n): The number  of set partitions (i.e. nops(SP(n)) done the THIRD way, using the EXPONENTIAL GENERATING FUNCTION FOR
#The Bell Numbers exp(exp(x)-1)
Bell3:=proc(n) local x:
n!*coeff(taylor(exp(exp(x)-1),x=0,n+1),x,n):
end:

#SeqBell0(N): The LIST of the first N terms of the sequence of Bell numbers done the stupid way
SeqBell0:=proc(N) local n:
[seq(nops(SP(n)),n=1..N)]:
end:



#SeqBell1(N): The LIST of the first N terms of the sequence of Bell numbers done via Bell1(n)
SeqBell1:=proc(N) local n:
[seq(Bell1(n),n=1..N)]:
end:


#SeqBell2(N): The LIST of the first N terms of the sequence of Bell numbers done via Bell2(n)
SeqBell2:=proc(N) local n:
[seq(Bell2(n),n=1..N)]:
end:


#SeqBell3(N): The LIST of the first N terms of the sequence of Bell numbers done via the Exponential generatining function
SeqBell3:=proc(N) local f,x,i:

f:=taylor(exp(exp(x)-1),x=0,N+1):
[seq(i!*coeff(f,x,i),i=1..N)]:
end:


###END From M11.txt



###From M10.txt
#M10.txt: Maple Code for Lecture 9 of Math 454,Combinatorics, Rutgers University, taught by Dr. Z. (Doron Zeilberger)
Help10:=proc():print(` MulPers(P1,P2) , InvPer(P) , GrabCycle(P,i), PtoC(P), CtoP(C), FunT(pi),  GG(S) , inv(pi), maj(pi) , FunT(pi)  `):
end:



#MulPers(P1,P2): Inputs permutations P1 and P2 of {1, ..., n} and outputs their product
#For example, try:
#MulPers([2,3,1],[3,1,2]);
MulPers:=proc(P1,P2) local n,i:
#The first few lines checks that the input is correct
 if not (type(P1,list) and type(P2,list)) then
  print(`Bad input`):
  RETURN(FAIL):
 fi:
n:=nops(P1):

if  not (convert(P1,set)={seq(i,i=1..n)} and convert(P2,set)={seq(i,i=1..n)} ) then
  print(`Bad input`):
  RETURN(FAIL):
 fi:

[seq(P2[P1[i]],i=1..n)]:


end:

#InvPer(P): The inverse of the permutation P
InvPer:=proc(P) local i,n,T:
 if not type(P,list) then
  print(`Bad input`):
  RETURN(FAIL):
 fi:

n:=nops(P):
if  not convert(P,set)={seq(i,i=1..n)}  then
  print(`Bad input`):
  RETURN(FAIL):
 fi:

#Recalling the list P is the function that maps i to P[i], we define a TABLE (a Maple Data structure) that
#maps P[i] to i

for i from 1 to n do
 T[P[i]]:=i:
od:

#We convert it back into a list
[seq(T[i],i=1..n)]:
end:

#GrabCycle(P,i): inputs a permutation P of {1, ...,n} and outputs the cycle belonging to i. 
#We use the convention that it starts with the largest entry
#Try
#GrabCycle([3,1,4,2],2);
GrabCycle:=proc(P,i) local C,j:

#We view the cycle belnging to i as a list, starting at i
C:=[i]:

#We find where i points to, call it j
j:=P[i]:

while j<>i do
#sooner or later we will wind back at i (since the permutation is finite) so we append the new
#arrivals to the list C and rename j
C:=[op(C),j]:
j:=P[j]:
od:

#We look at the place where it it is max
j:=max[index](C):

#we rewrite the cycle with the largest entry first

[op(j..nops(C),C),op(1..j-1,C)]:

end:

#PtoC(P): Inputs a permutation P outputs its list of cycles. Try:
#PtoC([2,1,3,4]);
PtoC:=proc(P) local n,StillToDo,L,i,C,L1,T:

n:=nops(P):

StillToDo:={seq(i,i=1..n)}:

#L is a dynamically construcete list of cycles in the permutatib P, we start with the empty list
#until no one is left
L:=[]:

while StillToDo<>{} do
#we pick the smallest survivor
 i:=min(op(StillToDo)):
 #we grab its cycle
 C:=GrabCycle(P,i):

# We append C to L
 L:=[op(L),C]:

#We update StillToDo by kicking out the members of C

StillToDo:= StillToDo minus convert(C,set):
od:

#So far L is a list of cycles but arranged in a random order
#Now we arrange them according to the convention that the first (largest) entries are increasing
#L1 is the list of first (largest) entries of each cycle:
L1:=[seq(L[i][1],i=1..nops(L))]:

#We now make a table where T[a] is the cycle that starts with  a
for i from 1 to nops(L) do
 T[L[i][1]]:=L[i]:
od:

#Now we sort L1
L1:=sort(L1):

#Now we rewrite L with the above convention
[seq(T[L1[i]],i=1..nops(L1))]:

end:


#FunT(P): The fundamental transformation
FunT:=proc(P) local P1,i:
P1:=PtoC(P):
[seq(op(P1[i]),i=1..nops(P1))]:
end:

#GG(S): inputs a pos.  integer n and a set of permutations in {1,..n} outputs the subgroup generated by them
GG:=proc(S) local Gr,s,g,Gr1,Gr2:
Gr:=S:
Gr1:=Gr union {seq(seq(MulPers(s,g), s in S), g in Gr)}:

while Gr<>Gr1 do
Gr2:=Gr1 union {seq(seq(MulPers(s,g), s in S), g in Gr1)}:
Gr:=Gr1:
Gr1:=Gr2:
od:
Gr1:

end:


#inv(pi): The number of inversion of a permutation pi
inv:=proc(pi) local n,i,j,co:
n:=nops(pi):
co:=0:

for i from 1 to n do
 for j from i+1 to n do
   if pi[i]>pi[j] then
    co:=co+1:
   fi:
 od:
od:

co:

end:


#maj(pi): The major index a permutation pi
maj:=proc(pi) local n,i,co:
n:=nops(pi):
co:=0:

for i from 1 to n-1 do
   if pi[i]>pi[i+1] then
    co:=co+i:
   fi:
 od:


co:

end:

#CtoP(C): Inputs a permutation in cycle structre , outputs it in 1-line notation
#Try([[1],[3,2]]);
CtoP:=proc(C) local i,L,n,T,C1,j:
L:=[seq(op(C[i]),i=1..nops(C))]:
n:=nops(L):

if convert(L,set)<>{seq(i,i=1..n)} then
  RETURN(FAIL):
fi:

for i from 1 to nops(C) do
 C1:=C[i]:
 for j from 1 to nops(C1)-1 do
   T[C1[j]]:=C1[j+1]:
 od:
 T[C1[nops(C1)]]:=C1[1]:
od:

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

###End From M10.txt