#Nathan Fox
#Homework 17
#I give permission for this file to be posted online

##Read old files
read(`C16joey.txt`):
read(`C17.txt`):

#Help procedure
Help:=proc():
print(` Options(G,pi,sig) , OneStep(G,pi) `):
print(` VerifyOneStep(G,i,DrNakamura) `):
print(` ManySteps(G,pi,NumberOfSteps) `):
print(` VerifyManySteps(G,DrNakamura) , C2(S,t,r,sig,K,T) `):
print(` DCX(S,t,r,sig,K,T,n,X) `):
end:

##Problem 1
#Options(G,pi,sig): Given a Hamiltonian graph G
#and a proof of its Hamiltonian pi (n=nops(pi)) and
#a permutation sig, encode the graph into boxes
Options:=proc(G,pi,sig) local Bij,edge,L,pi1,i,j,n,Htemp,H:
n:=nops(pi):
Bij:=[]:
for i from 1 to n do
    L:=[]:
    for j from 1 to n do
        L:=[op(L),0]:
    od:
    Bij:=[op(Bij),L]:
od:
for edge in G do
    Bij[sig[edge[1]]][sig[edge[2]]]:=1:
    Bij[sig[edge[2]]][sig[edge[1]]]:=1:
od:
Htemp:=subs({seq(i=sig[i],i=1..n)},pi):
H:=[]:
for i from 1 to nops(Htemp)-1 do
    H:=[op(H),[Htemp[i],Htemp[i+1]]]:
od:
H:=[op(H),[Htemp[nops(Htemp)],Htemp[1]]]:
[Bij,sig,H]:
end:

#Prover runs one step
OneStep:=proc(G,pi) local coin,sig,ops,n,boxes,hset,edge,i,j:
n:=nops(pi):
sig:=randperm(n):
ops:=Options(G,pi,sig):

coin:=rand(1..2)():

if coin=1 then
    return 1,[ops[1],ops[2]]:
else:
    #zero out boxes not involved
    hset:={}:
    for edge in ops[3] do
        hset:=hset union {edge}:
    od:
    boxes:=ops[1]:
    for i from 1 to n do
        for j from 1 to n do
            if not ([i,j] in hset) then
                boxes[i][j]:=0:
            fi:
        od:
    od:
    return 2,[boxes,ops[3]]:
fi:
end:

#Verifier runs one step
VerifyOneStep:=proc(G,i,DrNakamura) local n,sig,boxes,
    H,Bij,L,j,k,edge:

n:=nops(DrNakamura[2]):

if i=1 then
    boxes:=DrNakamura[1]:
    sig:=DrNakamura[2]:
    Bij:=[]:
    for k from 1 to n do
        L:=[]:
        for j from 1 to n do
            L:=[op(L),0]:
        od:
        Bij:=[op(Bij),L]:
    od:
    for edge in G do
        Bij[sig[edge[1]]][sig[edge[2]]]:=1:
        Bij[sig[edge[2]]][sig[edge[1]]]:=1:
    od:
    return evalb(boxes=Bij):
else
    boxes:=DrNakamura[1]:
    H:=DrNakamura[2]:
    #Check H has right length
    if n <> nops(boxes) then
        return false:
    fi:
    #Check H cycle
    if H[1][1] <> H[n][2] then
        return false:
    fi:
    for j from 1 to n-1 do
        if H[j][2] <> H[j+1][1] then
            return false:
        fi:
    od:
    #Check H in boxes
    for edge in H do
        if boxes[edge[1]][edge[2]] <> 1 then
            return false:
        fi:
    od:
    return true:
fi:
end:

##Problem 2
#Prover runs many steps
ManySteps:=proc(G,pi,NumberOfSteps) local i,ret:

ret:=[]:

for i from 1 to NumberOfSteps do
    ret:=[op(ret), [OneStep(G,pi)]]:
od:

end:

#Verifier runs many steps
VerifyManySteps:=proc(G,DrNakamura) local j,NumberOfSteps:
NumberOfSteps:=nops(DrNakamura):
for j from 1 to NumberOfSteps do
    if not VerifyOneStep(G,DrNakamura[j][1],DrNakamura[j][2]) then
        return false:
    fi:
od:
return true:
end:

##Problem 3
#Every graph I tried worked, including a RandHG(40,1/50) over
#50 steps

##Problem 4
with(Finance):
#Built-in Black-Scholes
C2:=proc(S,t,r,sig,K,T):
BlackScholesPrice(S,K,T-t,sig,r,0):
end:

#C and C2 give essentially the same values

##Problem 5
#DCX(S,t,r,sig,K,T,n,X): The approximate value (with n dicrete
#steps) using the Cox-Ross-Rubinstein discrete version of
#Black Scholes of a European Call option
#pgf of this
#(I really have no idea how to do this properly)
DCX:=proc(S,t,r,sig,K,T,n,X) local u,d,p,i:
u:=evalf(exp(sig*sqrt((T-t)/n))):
d:=1/u:
p:=(exp(r*(T-t)/n)-d)/(u-d):

evalf(exp(-r*(T-t)) * add( X^(binomial(n,i)*p^(n-i)*(1-p)^i*max(0,S*u^(n-i)*d^i-K)), i=0..n )):

end: