#C8.txt: Maple code for Lecture 8, Feb. 15, 2021, Dr. Z. ExpMath
Help8:=proc(): print(` RandSPn(n) , NuSPnF(n), RandSPnF(n), BellSeq(n)`):end:

#RandSPn(n): (from Yuxuan Yang's hw7): A randon set-partition of {1, ..,n}
RandSPn:=proc(n) local i,k,die1:
die1:=[seq(NuSPnk(n,i),i=1..n)]:
k:=LD(die1):
RandSPnk(n,k):
end:

NuSPnF:=proc(n) local k:
option remember:

if n=0 then
  RETURN(1):
fi:
add( binomial(n-1,k)*NuSPnF(n-1-k)     , k=0..n-1):
end:

#RandSPnF(n): A random set-partition of {1, ...,n} using the
#NEW approach
RandSPnF:=proc(n) local k,S,SP1,nEnemies,i1:
if n=0 then
 RETURN({}):
fi:

k:=LD([seq(binomial(n-1,k)*NuSPnF(n-1-k)     , k=0..n-1)])-1:

S:=Rnk(n-1,k):

nEnemies:=sort(convert({seq(i1,i1=1..n-1)} minus  S,list)):

SP1:=RandSPnF(n-1-k):

SP1:=subs( {seq(j=nEnemies[j],j=1..n-1-k)}  , SP1):

SP1 union { {n} union S}:
end:


#BellSeq(n): The first n bell numbers using the egf
#exp(exp(x)-1)=Sum(B(n)*x^n/n!,n=0..infinity);
BellSeq:=proc(n) local x,f:
 
f:=taylor(exp(exp(x)-1),x=0,n+1):

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

end:

#Start C4.txt:
Help4:=proc(): print(` Snk(n,k), NuSnk(n,k) , LD(L), Rnk(n,k)  `): end:

#Snk(n,k): Inputs NON-NEG. integer n and an integer k and OUTPUTS
#THE SET OF k-element subsets of {1,...,n}
#FOR EXAMPLE
#Snk(3,2) should be {{1,2},{1,3},{2,3}}
Snk:=proc(n,k) local S1,S2,s:
option remember:

if k<0 or k>n then
 RETURN({}):
fi:

if n=0 then
  if k=0 then
    RETURN({{}}):
  else
    RETURN({}):
  fi:
fi:

S1:=Snk(n-1,k):
S2:=Snk(n-1,k-1):
S1 union {seq( s union {n}, s in S2)}:
end:

#NuSnk(n,k): Inputs NON-NEG. integer n and an integer k and OUTPUTS
#THE NUMBER OF  k-element subsets of {1,...,n}
#FOR EXAMPLE
#NuSnk(3,2) should be 3
NuSnk:=proc(n,k) option remember:

if k<0 or k>n then
 RETURN(0):
fi:

if n=0 then
  if k=0 then
    RETURN(1):
  else
    RETURN(0):
  fi:
fi:

NuSnk(n-1,k)+ NuSnk(n-1,k-1):

end:



#Rnk(n,k): Inputs NON-NEG. integer n and an integer k and OUTPUTS
#A (UNIFORMLY AT-RANDOM)  k-element subsets of {1,...,n}
#FOR EXAMPLE
#Rnk(3,2) should output each of {1,2},{1,3},{2,3} with prob. 1/3

Rnk:=proc(n,k)  local c,S:
if k<0 or k>n then
 RETURN(FAIL):
fi:

if n=0 then
  if k=0 then
    RETURN({}):
   fi:
fi:

c:=LD([NuSnk(n-1,k), NuSnk(n-1,k-1)]):

if c=1 then
  RETURN(Rnk(n-1,k)):
else
 RETURN(Rnk(n-1,k-1) union {n}):
fi:

end:
#end from C4.txt
#From C7.txt
Help7:=proc(): print(` SPnk(n,k), SPn(n) , NuSPnk(n,k), NuSPn(n)`):
print(`RandSPnk(n,k)`):
end:

#SPnk(n,k): Inputs pos. integers k and n and 1<=k<=n
#outputs the SET of set-partitions of {1,...,n} with exactly
#k parts
SPnk:=proc(n,k) local SP,SP1,sp,i:
option remember:
if n=1 then
  if k=1 then 
   RETURN({{{1}}}):
  else
    RETURN({}):
  fi:
fi:

#Case (i): n is a singleton
SP1:=SPnk(n-1,k-1):

SP:={seq(sp union {{n}},sp in SP1)}:

#Case (ii) n has friends
SP1:=SPnk(n-1,k):

for sp in SP1 do
 for i from 1 to nops(sp) do
  SP:=SP union {{op(1..i-1,sp), sp[i] union {n},op(i+1..nops(sp),sp)}}:
 od:
 od: 
SP:
end:  

#SPn(n): The set of all set-partitions of {1,...,n}
#(same as combinat[setpartition](n))
SPn:=proc(n) local k:  { seq(  op(SPnk(n,k)),k=1..n)}:end:

#NuSPnk(n,k): The number of set-partitions of {1,...,n}
#into exactly k parts (aka Strirling numbers of the second kind)
NuSPnk:=proc(n,k) option remember:

if n=1 then
  if k=1 then
   RETURN(1):
  else
   RETURN(0):
 fi:
fi:
NuSPnk(n-1,k-1)+ k*NuSPnk(n-1,k):
end:

NuSPn:=proc(n) local k:
option remember:
add(NuSPnk(n,k),k=1..n):
end:

#LD(L): Inputs a list of positive integers L (of n:=nops(L) members)
#outputs an integer i from 1 to n with the prob. of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6
#output 2 with prob. 1/3
#output 3 with prob. 3/6=1/2
LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:

for i from 1 to n do
 su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:
end:

#RandSPnk(n,k): a random set partition of {1,..,n} with k parts
RandSPnk:=proc(n,k) local c,SP1,k1:
if n=1 then
 if k=1 then
   RETURN({{1}}):
  else
   RETURN(FAIL):
 fi:
fi:
c:=LD([NuSPnk(n-1,k-1), k*NuSPnk(n-1,k)]):

if c=1 then
 SP1:=RandSPnk(n-1,k-1):
 RETURN(SP1 union {{n}}):
 else
 SP1:=RandSPnk(n-1,k):
  k1:=rand(1..k)():
  RETURN( { op(1..k1-1,SP1), SP1[k1] union {n}, op(k1+1..k,SP1)}):
  
fi:
end:
##END OF C7.txt