#Nathan Fox
#Homework 9
#I give permission for this work to be posted online

#Read packages
with(`combinat`):

#Read procedures from class
read(`C9.txt`):

#Help procedure
Help:=proc():
    print(` AllStrategies(N) , RankStrategies(p, N, i) `):
    print(` ExpectedDurationEG(N, p, L) `):
end:

##PROBLEM 1##

#AllStrategies(N): inputs a positive integer N, larger than 1, and
#outputs all legal strategies, in other words it outputs the set
#of lists of integers 
#[a[1],a[2],..., a[N-1]] 
#such that 1<=a[i]<=BoS(N)[i] for every i from 1 to N-1
AllStrategies:=proc(N) local ret, L, i, j:
    ret:={}:
    L:=[1$N-1]:
    while true do
        ret:=ret union {L}:
        for i from nops(L)-1 by -1 to 2 do
            if L[i] < min(i, N-i) then
                L[i]:=L[i]+1:
                for j from i+1 to nops(L)-1 do
                    L[j]:=1:
                od:
                break:
            fi:
        od:
        if i = 1 then
            break:
        fi:
    od:
    return ret:
end:

##PROBLEM 2##

#RankStrategies(p, N, i): inputs a probability p, a positive
#integer N (denoting, as usual, the exit capital), and an integer
#i between 1 and N-1, denoting the starting capital, and outputs
#the ranking of the strategies (arranged as a list of lists) from
#best to worst, if you play the gambling game with probability p,
#and maximal capital N.
RankStrategies:=proc(p, N, i) local strats, j, S, prb, prbs, ret, s:
    strats:=AllStrategies(N):
    prbs:={}:
    for j from 1 to nops(strats) do
        prb:=PrEG(N, p, strats[j])[i]:
        if prb in prbs then
            S[prb]:=S[prb] union {strats[j]}:
        else
            prbs:=prbs union {prb}:
            S[prb]:={strats[j]}:
        fi:
    od:
    prbs:=sort(convert(prbs, list), `>`):
    ret:=[]:
    for prb in prbs do
        for s in S[prb] do
            ret:=[op(ret), s]:
        od:
    od:
    return ret:
end:

#For N=8, and p=9/19, do you get the same output (ranking) for
#all initial capitals i from 1 to 7? How about N=8 and p=10/19?
#No and no

##PROBLEM 3##

#If we look at the system defining P[i], we see that P[i] is the
#average of P[i-L[i]] and P[i+L[i]].  Notice that the probabilities
#i/N satisfy this average property, and they satisfy the boundary
#conditions P[0]=0, P[1]=1.  So, by uniqueness of solution to this
#system, these must be the probabilities for any strategy L.

##PROBLEM 4##

#ExpectedDurationEG(N, p, L): Generalizes
#ExpectedDurationE(N, p, L) from hw8.txt
ExpectedDurationEG:=proc(N, p, L) local E, eq, var, i, sol:
    var:={seq(E[i],i=0..N)}:
    eq:={E[0]=0, E[N]=0, seq(E[i]=1+p*E[i+L[i]]+(1-p)*E[i-L[i]], i=1..N-1)}:
    sol:=solve(eq, var):
    return subs(sol, [seq(E[i],i=1..N-1)]):
end:

##PROBLEM 5##

#Let M be the odd part of N (i.e. N divided by as large of a power
#of 2 as possible).
#ExpectedDurationEG(N, 1/2, BoS(N))[i] takes the following values:
#2, if M does not divide i
#(2^k-1)/2^(k-1), if M divides i*2^k but M does not divide i*2^(k-1)
#
#Proof: If M does not divide i, then the position reached from i
#will not be divisible by M either.  So, each move in this game
#will have probability 1/2 of ending the game.  This is a
#geometric distribution, so the expected game length is 2.
#If N is even and i=N/2, then the game immediately ends, so the
#expected duration is 1.
#Now, assume that M  divides i*2^k but M does not divide i*2^(k-1).
#Then, with probability 1/2, the game ends immediately, and with
#probability 1/2 we move to a position j such that M divides
#j*2^(k-1) but M does not divide j*2^(k-2).  So, by induction, the
#expected duration of the game is (1/2)*1+(1/2)*(2^(k-1)-1)/2^(k-2)
#=1/2+(2^(k-1)-1)/2^(k-1)=(2^(k-1)+2^(k-1)-1)/2^(k-1)
#=(2^k-1)/2^(k-1), as required.