#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 17
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

# From my c17.txt:

# c17.txt
# 28 Mar 2013
# Zero-knowledge proofs

with(combinat):
with(RandomTools[MersenneTwister]):
randomize():


RandG := proc(n,p) local S, i, j:
    S := {}:
    for i from 1 to n do:
        for j from i+1 to n do:
            if GenerateFloat() < evalf(p) then:
                S := S union {{i,j}}:
            fi:
        od:
    od:
    return S:
end:



RandHG := proc(n,p) local pi, G:
    pi := randperm(n):
    G := {seq({pi[i], pi[i+1]}, i=1..n-1), {pi[n], pi[1]}}:
    return G union RandG(n,p), pi:
end:



#############
# Problem 1 #
#############

# This code implements the correct procedure, but requires
# trust because we have to trust that the options are actually
# generated from G1.

# There is no way to do the procedure without trust unless we have
# encryption.

Options := proc(G1, pi, sig) local i, n, H1, H:
    n := nops(pi):
    H1 := subs({seq(i=sig[i], i=1..n)}, pi):
    H := {seq({H1[i], H1[i+1]}, i=1..n-1), {H1[n], H1[1]}}:
    return [[sig, G1], H]:
end:

OneStep := proc(G, pi) local G1, n, r, sig, ops:
    n := nops(pi):
    sig := randperm(n):
    G1 := subs({seq(i=sig[i], i=1..n)}, G):
    ops := Options(G1, pi, sig):
    r := rand(1..2)():
    return [r, ops[r]]:
end:

VerifyStep := proc(G, answer) local r, G1, sig, G2, H, pi, i, n, c:
    r := answer[1]:
    c := answer[2]:
    if r = 1 then:
        sig := c[1]:
        n := nops(sig):
        G1 := c[2]:
        G2 := subs({seq(i=sig[i], i=1..n)}, G):
        return evalb(G1 = G2):
    elif r = 2 then:
        H := c:
        n := nops(H):
        return isCycle(H, n):
    else:
        print("oops!"):
    fi:
end:

isCycle := proc(G,n) local e, i, incidences, counts, S, j, p, nx:
    incidences := [seq(e[1], e=G), seq(e[2], e=G)]:
    counts := {seq(numboccur(i, incidences), i=1..n)}:
    if counts <> {2} then:
        return false:
    fi:
    S := {}:
    j := 1:
    while S <> {seq(i, i=1..n)} do:
        p := map(y->(op(y minus {j})), select(x->evalb(j in x), G)):
        nx := select(x->not evalb(x in S), p):
        if nx = {} then:
            return false:
        fi:
        S := S union {nx[1]}:
        j := nx[1]:
    od:
    return true:
end:



#############
# Problem 2 #
#############

VerifyManySteps := proc(G, pi, N) local i:
    for i from 1 to N do:
        if (VerifyStep(G, OneStep(G, pi)) = false) then:
            return false:
        end:
    od:
    return true:
end:

#############
# Problem 3 #
#############

# It works (with high probability!).

#############
# Problem 5 #
#############

DCX:=proc(S,t,r,sig,K,T,n, X) local u,d,p,i, price:
    u:=evalf(exp(sig*sqrt((T-t)/n))):
    d:=1/u:
    p:=(exp(r*(T-t)/n)-d)/(u-d):
    price := evalf(exp(-r*(T-t)) * add( binomial(n,i)*p^(n-i)*(1-p)^i*max(0,S*u^(n-i)*d^i-K), i=0..n )):
    return price, sort(evalf(exp(-r*(T-t)) * add( binomial(n,i)*p^(n-i)*(1-p)^i*X^max(-price,(S*u^(n-i)*d^i-K)-price), i=0..n ))):
end: