#OK to post homework
#Victoria Chayes, 2/7/21, 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:
with(combinat):


#1 Prove that NuSnk(n,k)=n!/(k!*(n-k)!). 

# We run the Maple command
#simplify((n-1)!/(k!*(n-1-k)!)+(n-1)!/((k-1)!*(n-k)!))
#                         factorial(n)         
#                 -----------------------------
#                 factorial(k) factorial(n - k)


#2  Write a procedure RnkFast(n,k) by replacing the line LD([NuSnk(n-1,k), NuSnk(n-1,k-1)]) By LD([n-k,k]) First explain why it is valid, and then find out, for fairly large n and k (say n=5000, k=2000) whether it is indeed faster. Why?

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:

# This is a valid change because we have proved that NuSnk(n,k)=n!/(k!*(n-k)!). Our loaded die perviously had input of (n-1)!/(k!*(n-1-k)!) * (k!*(n-k)!)/n! = n-k/n chance of rolling 0, and (n-1)!/((k-1)!*(n-k)! * (k!*(n-k)!)/n! = k/n chance of rolling 1.
# Testing:
#time(Rnk(5000,2000))
#                            255.826
#time(RnkFast(5000,2000))
#                             13.366
#


#3   Write a procedure 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-element subsets of {1,...,n}. Experiment it with various n and k, and K=1000.

EstimateAveMax:=proc(n,k,K):
evalf(add(max(RnkFast(n,k)),i=1..K)/K):
end:

# I did the following testing for n=3, 5, 10:
#seq(EstimateAveMax(3,i,1000),i=1..3)
#                  2.005200000, 2.674500000, 3.
#
#seq(EstimateAveMax(5,i,1000),i=1..5)
#     3.005200000, 4.004400000, 4.505000000, 4.797300000, 5.
#
#seq(EstimateAveMax(10,i,1000),i=1..10)
#     5.343000000, 7.310000000, 8.251000000, 8.764000000, 9.200000000, 9.416000000, 9.630000000, 9.783000000, 9.907000000, 10.


#4 Write a procedure 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 sum(S) over all k-elemenet subsets of {1,...,n}. Experiment it with various n and k, and K=1000.

EstimateAveSum:=proc(n,k,K):
evalf(add(add(RnkFast(n,k)),i=1..K)/K):
end:

# I did the following testing for n=3, 5, 10:
#seq(EstimateAveSum(3,i,1000),i=1..3)
#                  1.994000000, 4.003000000, 6.
#
#seq(EstimateAveSum(5,i,1000),i=1..5)
#    2.957000000, 6.089000000, 9.054000000, 12.00200000, 15.
#
#seq(EstimateAveSum(10,i,1000),i=1..10)
#    5.229000000, 10.90300000, 16.37100000, 22.12700000, 27.72200000, 32.72800000, 38.36500000, 44.08300000, 49.45500000, 55.



#5 As you know, k-element subsets of {1,...,n} are in one-to-one correspondence with n-letters words in the alphabet {0,1} with exactly k ones and n-k zeroes. Adapt Snk(n,k), NuSnk(n,k), and Rnk(n,k) to write procedure Wnk(a,b), NuWnk(a,b), RWnk(a,b) that input non-negative integers a,b, and outputs, respectively, the set of all words in {1,2} with a 1s, b 2s, their number, and a random such word.

Wnk:=proc(a,b) local S1,S2,s:
option remember:
if a<0 or b<0 then
 return({}):
fi:
if a=0 then
  return({[2$b]}):
fi:
if b=0 then
  return({[1$a]}):
fi:
S1:=Wnk(a-1,b):
S2:=Wnk(a,b-1):
{seq([op(s), 1], s in S1)} union {seq([op(s), 2], s in S2)}:
end:

NuWnk:=proc(a,b) 
option remember:
if a<0 or b<0 then
  return(0):
fi:
if a=0 then
  return(1):
fi:
if 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=0 then
  return([2$b]):
fi:
if b=0 then
  return([1$a]):
fi:
c:=LD([NuWnk(a-1,b), NuWnk(a,b-1)]):
if c=1 then
   return([1,op(RWnk(a-1,b))]):
else
   return([2,op(RWnk(a,b-1))]):
fi:
end:


#6 Write a procedure W3nk(a,b,c), NuW3nk(a,b,c), RW3nk(a,b,c) that input non-negative integers a,b, c and outputs, respectively, the set of all words in {1,2,3} with a 1s, b 2s, c 3s, their number, and a random such word.

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=0 and b=0 then
  return({[3$c]}):
fi:
if b=0 and c=0 then
  return({[1$a]}):
fi:
if a=0 and c=0 then
  return({[2$b]}):
fi:

S1:=W3nk(a-1,b,c):
S2:=W3nk(a,b-1,c):
S3:=W3nk(a,b,c-1):
{seq([op(s), 1], s in S1)} union {seq([op(s), 2], s in S2)} union {seq([op(s), 3], 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=0 and b=0 then
  return(1):
fi:
if b=0 and c=0 then
  return(1):
fi:
if a=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 k,S:
if a<0 or b<0 or c<0 then
  return(FAIL):
fi:
if a=0 and b=0 then
  return([3$c]):
fi:
if b=0 and c=0 then
  return([1$a]):
fi:
if a=0 and c=0 then
  return([2$b]):
fi:
k:=LD([NuW3nk(a-1,b,c), NuW3nk(a,b-1,c), NuW3nk(a,b,c-1)]):
if k=1 then
   return([1,op(RW3nk(a-1,b,c))]):
elif k=2 then
   return([2,op(RW3nk(a,b-1,c))]):
else
   return([3,op(RW3nk(a,b,c-1))]):
fi:
end:

############################
### Programs From C4.txt ###
############################

Help4:=proc(): print(`NuSnk(n,k) , LD(L), Rnk(n,k)  `): end:

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