# OK to post homework
# Lucy Martinez, 04-17-22, Assignment 22
read `C22.txt`:
#----------------------Problem 1-----------------------#
# By the definition of ai(i), the sets {ai(i),i=1..infinity} and {ai(i)+i,i=1..infinity} are disjoint and
# completey cover the set of positive integers (exactly!, i.e. no overlap). Assuming that ai(i)=trunc(i*alpha)
# for some real number alpha (this assumption is not obvious, and once the only eligible alpha is found needs
# a rigorous proof). Indeed there is only one thing alpha can be. Prove that if such an alpha exists,
# it must be the famous Golden Ratio.

 

#----------------------Problem 2-----------------------#
# [Optional challenge] Try to reprove Wythoff's theorem that ai(i)=aiC(i). If you give up, look it up 
# and try to understand the proof.


#----------------------Problem 4-----------------------#
# Using LaG(G), write a procedure LP(G) that inputs a directed graph w/o cycles, G (let n:=nops(G))
# and outputs the subest of {1,...,n} of its losing positions.

LP:=proc(G) local L, S, i:
L:=LaG(G):
S:={}:
for i from 1 to nops(L) do
 if L[i]=0 then
  S:=S union {i}:
 fi:
S:
od:
 
end:
LaG([{2, 4}, {4, 5}, {4}, {5}, {}]);
#             [0, 1, 0, 1, 0]


LP([{2, 4}, {4, 5}, {4}, {5}, {}]);
#                  {1, 3, 5}


#----------------------Problem 5-----------------------#
# Write a procedure AvNuLP(n,p,K) that inputs a positive integer n, a rational number p between 0 and 1 
# and a large integer K and outputs an approximation to the average size of the set of losing positions,
# for a random directed graph (w/o cycles) with n vertices, by averaging over K random graphs gotten 
# from RDG(n,p). What did you get for AvNuLP(n,p,10000) for 3*4=12 cases n=10,50, 100; p=1/5,2/5,3/5,4/5

AvNuLP:=proc(n,p,K) local i,sum:
sum:=0:
for i from 1 to K do
 sum:=sum + nops(LP(RDG(n,p))):
end:
evalf(sum/K):
 
end:

## n=10
[seq(AvNuLP(10,i,10000),i=1/5..4/5,1/5)];
#      [5.505200000, 3.796900000, 2.806100000, 2.043400000]


## n=50
[seq(AvNuLP(50,i,10000),i=1/5..4/5,1/5)];
#      [11.45810000, 6.685000000, 4.463300000, 3.001700000]


## n=100
[seq(AvNuLP(100,i,10000),i=1/5..4/5,1/5)];
#      [14.35880000, 7.981700000, 5.192900000, 3.429700000]