#OK to post
#Yuxuan Yang, Feb 4th, Assignment 4 

with(combinat):

#1
#NuSnk(n,k)=NuSnk(n-1,k)+ NuSnk(n-1,k-1)
#n!/(k!(n-k)!)=(n-k+k)(n-1)!/(k!(n-k)!)=(n-1)!/(k!(n-k-1)!)+(n-1)!/((k-1)!(n-k)!)
#check boundaries



#2
#Because NuSnk(n-1,k)/NuSnk(n-1,k-1)=(n-k)/k, then LD([NuSnk(n-1,k), NuSnk(n-1,k-1)]) is the same as LD([n-k,k])
RnkFast:=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([n-k,k]):

if c=1 then
  RETURN(RnkFast(n-1,k)):
else
 RETURN(RnkFast(n-1,k-1) union {n}):
fi:
end:
#By testing, it is obviously faster.



#3
#EstimateAveMax(n,k,K)
#that uses Rnk(n,k) (or if you prefer, your RnkFast(n,k)) K times,
#and estimates the average of the quantity (max(S)) over all k-elemenet subsets of {1,...,n}.
EstimateAveMax:=proc(n,k,K) local i,m:
m:=0:
for i from 1 to K do
m:=m+max(op(RnkFast(n,k))):
od:
m/K:
end:






#4
#EstimateAveSum(n,k,K)
#that uses Rnk(n,k) (or if you prefer, your RnkFast(n,k)) K times,
#and estimates the average of the quantity (max(S)) over all k-elemenet subsets of {1,...,n}.
EstimateAveSum:=proc(n,k,K) local i,j,m,R:
m:=0:
for i from 1 to K do
R:=RnkFast(n,k):
m:=m+add(R[j],j=1..k):
od:
m/K:
end:



#5
Wnk:=proc(a,b) local w1,w2,w,i:
option remember:
if a<0 or b<0 then 
 RETURN({}):
fi:
if a+b=0 then 
 RETURN({[]}):
fi: 
w1:=Wnk(a-1,b):
w2:=Wnk(a,b-1):
{seq([op(w1[i]),1],i=1..nops(w1)),seq([op(w2[i]),2],i=1..nops(w2))}:
end:



NuWnk:=proc(a,b) option remember:
if a<0 or b<0 then
 RETURN(0):
fi:
if a+b=0 then
  RETURN(1):
fi:
NuWnk(a-1,b)+ NuWnk(a,b-1):
end:


RWnk:=proc(a,b)  local c,S:
if a<0 or b<0 then
 RETURN(FAIL):
fi:

if a+b=0 then
   RETURN([]):
fi:

c:=LD([NuWnk(a-1,b), NuWnk(a,b-1)]):

if c=1 then
 RETURN([op(RWnk(a-1,b)),1]):
else
 RETURN([op(RWnk(a,b-1)),2]):
fi:

end:

#6



W3nk:=proc(a,b,c) local w1,w2,w3,w,i:
option remember:
if a<0 or b<0 or c<0 then 
 RETURN({}):
fi:
if a+b+c=0 then 
 RETURN({[]}):
fi: 
w1:=W3nk(a-1,b,c):
w2:=W3nk(a,b-1,c):
w3:=W3nk(a,b,c-1):
{seq([op(w1[i]),1],i=1..nops(w1)),seq([op(w2[i]),2],i=1..nops(w2)),seq([op(w3[i]),3],i=1..nops(w3))}:
end:



NuW3nk:=proc(a,b,c) option remember:
if a<0 or b<0 or c<0 then
 RETURN(0):
fi:
if a+b+c=0 then
  RETURN(1):
fi:
NuW3nk(a-1,b,c)+ NuW3nk(a,b-1,c)+NuW3nk(a,b,c-1):
end:


RW3nk:=proc(a,b,c)  local co,S:
if a<0 or b<0 or c<0 then
 RETURN(FAIL):
fi:

if a+b+c=0 then
   RETURN([]):
fi:

co:=LD([NuW3nk(a-1,b,c), NuW3nk(a,b-1,c),NuW3nk(a,b,c-1)]):

if co=1 then
 RETURN([op(RW3nk(a-1,b,c)),1]):
elif co=2 then
 RETURN([op(RW3nk(a,b-1,c)),2]):
else
 RETURN([op(RW3nk(a,b,c-1)),3]):
fi:

end:


















#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:

#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:

#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: