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

#Use packages
with(`combinat`):

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

#Help procedure
Help:=proc():
    print(` DoNs(N) , AveDuration(N, K) , CheckPie(L) , DoNg(N, p) `):
end:

##PROBLEM 3##

#DoNs(N): Winning a dollar with probability 1-1/2^N
#
#INPUT:
#a positive integer N, the max number of rounds
#
#OUTPUT: the number of rounds it takes to win a dollar
#or N if you lose
DoNs:=proc(N) local c, coin, ST, debt, i:
    coin:=rand(0..1):
    ST:=1:
    debt:=0:
    
    for i from 1 to N do
        c:=coin():
        if c = 1 then
            return i:
        else
            debt:=debt + ST:
            ST:=ST * 2:
        fi:
    od:
    return N:
end:

##PROBLEM 4##

#AveDuration(N, K): Runs DoNs(N) K times
#and computes the average duration.
#
#INPUT:
#a positive integer N, the max number of rounds
#a positive integer K, the number of iterations
#
#OUTPUT: the average duration
AveDuration:=proc(N, K) local i:
    return add(DoNs(N), i=1..K)/K:
end:

#For AveDuration(100,1000), I got 493/250, which is 1.972

##PROBLEM 5##

#Assume the probability of winning is p. Then, the expected
#duration is sum(i*(1-p)^(i-1)*p, i=1..infinity).  This equals 1/p.
#So, the expected duration is 1/p.  (This makes sense, since it's
#a geometric distribution with parameter p.)

##PROBLEM 6##

#CheckPie(L): inputs a list of sets L of length k
#and checks the Inclusion-Exclusion Principle for k sets: 
#|L1 union ... union Lk|=|L1|+ ... + |Lk|
#- (|L1 intersect L2|+|L1 intersect L3|+ ...) + ...
#+ (-1)^(k-1)|L1 intersect L2 intersect ... intersect Lk|
#
#INPUT:
#a list of sets L
#
#OUTPUT:
#the size of the union of the sets in L = the right-hand side of
#Inclusion-Exclusion (the literal equality expression) and
#the boolean value of said expression.  (Should always be true.)
CheckPie:=proc(L) local k, lside, i, j, LI, rside, comb:
    k:=nops(L):
    lside:=nops(`union`(op(L))):
    rside:=0:
    for i from 1 to k do
        comb:=firstcomb(k, i):
        while comb <> FAIL do
            rside:=rside + (-1)^(i-1)*nops(`intersect`(seq(L[j], j in comb))):
            comb:=nextcomb(comb, k):
        od:
    od:
    return lside = rside, evalb(lside = rside):
end:

##PROBLEM 7##

#DoNg(N, p) Winning a dollar with probability 1-p^N
#
#INPUT:
#a positive integer N, the max number of rounds
#a rational number p, the probability of winning in a given round
#
#OUTPUT: the number of rounds it takes to win a dollar
#or N if you lose
DoNg:=proc(N, p) local c, coin, ST, debt, i:
    coin:=rand(1..denom(p)):
    ST:=1:
    debt:=0:
    
    for i from 1 to N do
        print(`Right now the stake of a coin toss is`, ST):
        print(`and your debt is`, debt):
        c:=coin():
        if c <= numer(p) then
            print(`Yea, you won then`, ST, `dollars`):
            print(`after you pay your debt of`, debt, `you exit with`, ST-debt, `dollars`):
            return i:
        else
            print(`What a shame, you lost`, ST, `dollars`):
            debt:=debt + ST:
            ST:=ST * 2:
        fi:
    od:
    print(`Poor you, you now owe`, debt, `dollars`):
    return N:
end: