Lower and upper bounds for a gambler's ruin problem

gambler.py

#!/usr/bin/env python

from __future__ import division
from numpy import *

def iter(n,boundary):
  p = ones(n)
  b = boundary.copy()
  j = 3*arange(1,n)
  for i in xrange(1000):
    b[:n] = p
    p[1:] = .5*(p[:-1]+b[j])
  return p

print repr(iter(1000,zeros(3000))[1])
print repr(iter(1000,((1-sqrt(5))/2)**arange(3000))[1])

Output:

cone:~% ./gambler 
0.57477536382958938
0.57477536382958938