#OK to post homework
#George Spahn, 2/07/2021, Assignment 4

# 1. We wish to prove that the number of k subsets of {1..n} is 
# n!/k!/(n-k)!.  We showed in class that this number is also equal to the 
# number of k subsets of {1..n-1} plus the number of k-1 subsets of 
# {1..n-1}. If we show that the initial conditions are correct and that 
# n!/k!/(n-k)! also satisfies the recurrence relation, we will be done.

# (n-1)!/k!/(n-1-k)! + (n-1)!/(k-1)!/(n-k)! = 
# (n-1)! * (n-k) / k! / (n-k)! + (n-1)! * k / k! / (n-k)! = 
# (n-1)! * n / k! / (n-k)! =
# n! / k! / (n-k)!




#RnkFast(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

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:

#To see why this fast version is correct, consider the ratio
# 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
# = n-k : k
# Thus this algorithm also decides whether to add each element with the
# correct proportion.
# It is faster because it does not need to count the number of subsets of 
# each size to compute the correct ratio.


#EstimateAveMax(n,k,K)  Uses 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 Maxes,i:
Maxes:=[0$K]:
for i from 1 to K do
 Maxes[i]:= max(RnkFast(n,k)):
od:
add(Maxes)/K:
end:

#EstimateAveSum(n,k,K)  Uses RnkFast(n,k) K times, and estimates the average
# of the quantity (sum(S)) over all k-elemenet subsets of {1,...,n}
EstimateAveSum:=proc(n,k,K) local Sums,i:
Sums:=[0$K]:
for i from 1 to K do
 Sums[i]:= add(RnkFast(n,k)):
od:
add(Sums)/K:
end:
#The exact value of the average should be k*(n+1)/2. We can apply linearity
#of expectation and note that the average for each element of a k-subset is 
#(n+1)/2


Wnk:=proc(a,b) local k,n,S1,S2,s:
option remember:
k:=a:
n:=a+b:
if k<0 or k>n then
 RETURN({}):
fi:

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

S1:=Wnk(a-1,b):
S2:=Wnk(a,b-1):
{seq( ["a",op(s)], s in S1)} union {seq( ["b",op(s)], s in S2)}:
end:

NuWnk:=proc(a,b) local n,k:
option remember:
k:=a:
n:=a+b:
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:

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

end:



RWnk:=proc(a,b)  local k,n,c,S:
k:=a:
n:=a+b:
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(["b", op(RWnk(a,b-1)) ]):
else
 RETURN(["a", op(RWnk(a-1,b)) ]):
fi:

end:

W3nk:=proc(a,b,c) local S1,S2,S3,s:
option remember:

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

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

S1:=W3nk(a-1,b,c):
S2:=W3nk(a,b-1,c):
S3:=W3nk(a,b,c-1):
{seq( ["a",op(s)], s in S1)} union {seq( ["b",op(s)], s in S2)}
      union {seq( ["c",op(s)], s in S3)}:
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 d:


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

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

if d=1 then
  RETURN(["a",op(RW3nk(a-1,b,c))]):
fi:

if d=2 then
  RETURN(["b",op(RW3nk(a,b-1,c))]):
fi:
if d=3 then
  RETURN(["c",op(RW3nk(a,b,c-1))]):
fi:
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: