#OK to post homework
#Wanying Rao, 02/06/2021, Assignment 4

Help := proc(): print(`RnkFast(n,k), EstimateAveMax(n,k,K), EstimateAveSum(n,k,K), Wnk(a,b), NuWnk(a,b), RWnk(a,b), W3nk(a,b,c), NuW3nk(a,b,c), RW3nk(a,b,c)`):end:

#1
#Let NuSnk(n,k)=n!/(k!*(n-k)!).
#We can verify that NuSnk(n-1,k) + NuSnk(n-1,k-1) = NuSnk(n,k).
#NuSnk(n-1,k) + NuSnk(n-1,k-1) = (n-1)!/(k!*(n-1-k)!)+(n-1)!/((k-1)!*(n-k)!)
#                              = ((n-1)!*(n-k)+(n-1)!*k)/(k!*(n-k)!)
#                              = n!/(k!*(n-k)!)
#                              = NuSnk(n,k)

#2
#Since LD([NuSnk(n-1,k),NuSnk(n-1,k-1)]) outputs 1 or 2 with the probability being #proportional to NuSnk(n-1,k) and NuSnk(n-1,k-1), we can then replace the input list
#by a simplified list with the same proportion of each element respectively.
#NuSnk(n-1,k)=(n-1)!/(k!*(n-1-k)!)
#NuSnk(n-1,k-1)=(n-1)!/((k-1)!*(n-k)!)
#NuSnk(n-1,k)/NuSnk(n-1,k-1)=(k-1)!(n-k)!/k!(n-1-k)!
#                           =n-k/k
#Therefore, the ratio of NuSnk(n-1,k) and NuSnk(n-1,k-1)=n-k/k and 
#LD([NuSnk(n-1,k),NuSnk(n-1,k-1)]) = 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:

#RnkFast(5000,2000) outputs the results in about 2 seconds. 
#I interrupted Rnk(5000,2000) since it took more than 90 seconds.
#RnkFast is indeed faster than Rnk because we do not need to go into the recursive #procedure NuSnk. 

#3
#(i)
EstimateAveMax := proc(n,k,K) local i,m:
m := 0:
for i from 1 to K do
 m := m + max(Rnk(n,k)):
od:
RETURN(m/K):
end:


#4
#(i)
EstimateAveSum := proc(n,k,K) local i,s:
s := 0:
for i from 1 to K do
 s := s+ sum(Rnk(n,k)):
od:
RETURN(s/K):
end:


#5
Wnk := proc(a,b) local W1,W2,w1,w2:
option remember:

if b<0 or a<0 then 
 RETURN({}):
fi:

if a=0 and b=0 then
 RETURN({String()}):
fi:

W1:=Wnk(a-1,b):
W2:=Wnk(a,b-1):
{seq(String(op(w1),1), w1 in W1)} union {seq(String(op(w2),2), w2 in W2)}:
end:


NuWnk := proc(a,b) 
option remember:

if b<0 or a<0 then 
 RETURN(0):
fi:

if a=0 and b=0 then
 RETURN(1):
fi:

NuWnk(a-1,b) + NuWnk(a,b-1):
end:

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

if a=0 and b=0 then
 RETURN(String()):
fi:

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

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


#6
W3nk := proc(a,b,c) local W1,W2,W3,w1,w2,w3:
option remember:

if b<0 or a<0 or c<0 then 
 RETURN({}):
fi:

if a=0 and b=0 and c=0 then
 RETURN({String()}):
fi:

W1:=W3nk(a-1,b,c):
W2:=W3nk(a,b-1,c):
W3:=W3nk(a,b,c-1):
{seq(String(op(w1),1), w1 in W1)} union {seq(String(op(w2),2), w2 in W2)} union {seq(String(op(w3),3), w3 in W3)}:
end:


NuW3nk := proc(a,b,c) 
option remember:

if b<0 or a<0 or c<0 then 
 RETURN(0):
fi:

if a=0 and b=0 and 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 d:
if b<0 or a<0 or c<0 then 
 RETURN(FAIL):
fi:

if a=0 and b=0 and c=0 then
 RETURN(String()):
fi:

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

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

end:



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