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

# Nathaniel Shar
# HW 2
# Experimental Mathematics

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

# From c2.txt:

fGi:=proc(G,i) local M,m:
option remember:

    if not (type(i,integer) and type(G,list) and i>=1 and i<=nops(G))
    then
        print(`bad input`):
        RETURN(FAIL):
    fi:


    if G[i]={} then
        RETURN(0):
    fi:

    M:=G[i]:

    if {seq(fGi(G,m), m in M)}={1} then
        RETURN(0):
    else
        RETURN(1):
    fi:
end:

#WL(G): Winning list
#Inputs a directed G (representing a game)
#and outputs a list of length nops(G) such that
#L[i]=0 or 1 according to whether position (vertex) i
#is a losing or winning position
WL:=proc(G) local i:
    [seq(fGi(G,i),i=1..nops(G))]:
end:

#RandGame(n,k): inputs pos. integers n and k
#and outputs a random digraph that has the
#property that out of vertex i all the labels
#are less i, and there are (usually) k outgoing
#edges
RandGame:=proc(n,k) local L,i,tomas,j:
    L:=[{}]:
    for i from 2 to n do
        tomas:={seq(rand(1..i-1)(),j=1..k)}:
        L:=[op(L), tomas]:
    od:
    L:
end:




#############
# Problem 4 #
#############

AveLosing := proc(n,k,K) local j:
    return evalf(add(j,j=seq(numboccur(0, WL(RandGame(n,k))), i=1..K))/(n*K)):
end:

# We would expect, for large n, that the average proportion of losing
# positions would approach the root of P(x) = (1-x)^k - x that lies
# in [0,1]. A heuristic argument: Let p be the probability that the
# nth position is winning. This happens only if n randomly-chosen
# previous positions are all losing. Assuming that n is large, we can
# guess that the probability of this is roughly (1-p)^n, so p =
# (1-p)^n.

# An actual proof is harder. Let p_i be the probability that the ith
# position is winning. Then p_1 = 0 and p_i = (1-p_1 + ... +
# 1-p_{i-1})^k/(i-1)^k for i > 1. It should be possible to show that
# the p_i converge (at which point it is clear that the only fixed
# point is the aforementioned root), but I don't know how offhand.

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

ProbWinStupid := proc(G, i) local j: option remember:
    if G[i] = {} then
        return 0:
    else:
        return add(1-ProbWinStupid(G,j), j=G[i])/nops(G[i]);
    fi:
end:

#############
# Problem 6 #
#############

# Let p_i be the probability of winning with i coins in the heap.
# It is easy to see that p_0 = 0, p_1 = 1. Then, if 2 <= i <= k, p_i =
# 1/2, because you have equal chances to reach 0 and 1, and otherwise
# you reach some p_j with j between 2 and k, where the win probability
# is 1/2 for your opponent (by induction).

# For i >= k+1,
#
#     p_i = (1/k)((1-p_{i-1}) + (1-p_{i-2}) +  ... + (1-p_{i-k})).
#
# This is an inhomogeneous linear recurrence. We can solve it using
# any method. However, the characteristic polynomial is

#     kx^k + x^{k-1} + x^{k-2} + ... + 1

# whose roots do not appear to have a nice closed form. For the
# inhomogeneous part of the equation, we need to add a particular
# solution. Fortunately, an easy one is available: p_i = 1/2.

# As mentioned before, our initial conditions for this recurrence are
# p_0 = 0, p_1 = 1, p_2 = ... = p_k = 1/2. So all that remains is to
# solve the system of equations

# a_1      + a_2      + ... + a_k      = 0
# a_1r_1   + a_2r_2   + ... + a_kr_k   = 1
# a_1r_1^2 + a_2r_2^2 + ... + a_kr_k^2 = 1/2
#          .          . .   .          .
#          .          .  .  .          .
#          .          .   . .          .
# a_1r_1^k + a_2r_2^k + ... + a_kr_k^k = 1/2

# where {r_j} are the roots of the characteristic polynomial. (This is
# easy using the formula for the inverse of the Vandermonde matrix,
# which is a little too cumbersome to write here.) Then we can simply
# write

# p_i = a_1r_1^i + ... + a_kr_k^i + 1/2.