#!/usr/bin/python -Wall

# ================================================================
# John Kerl
# kerl at math dot arizona dot edu
# 2008-05-12
# ================================================================

from __future__ import division # 1/2 = 0.5, not 0.
import random
import kerlutil
from math import *

# ----------------------------------------------------------------
# 1D Brownian motion
# As is standard in Python, use an optional argument to implement what would be
# called a "static auto" in C or a "save parameter" in Fortran -- a local
# variable which is saved between calls.
def B(t, dt, uprev=[0.0]):
	du = random.gauss(0, sqrt(dt))
	uprev[0] += du
	return uprev[0]

# ----------------------------------------------------------------
# int_a^b X_t dB_t = sum_{i=0}^{n-1} X_{t_i} (B_{t_{i+1}} - B_{t_i}).
#                  = sum_{i=1}^{n} X_{t_{i-1}} (B_{t_i} - B_{t_{i-1}}).
# First, try X=B.
# int_0^T B_t dB_t = sum B_{ti} Delta B_{ti}

def do_integral():
	T = 10.0
	dt = 0.01
	ts = kerlutil.mfrange(0, dt, T)
	nt = len(ts)

	sum = 0.0

	Btprev = 0.0
	Bt = 0.0
	for t in ts:
		Delta_B = Bt - Btprev
		#sum += Bt * Delta_B
		sum += Btprev * Delta_B
		#print "%11.7f %11.7f %11.7f %11.7f" % (t, Btprev, Bt, sum)
		Btprev = Bt
		Bt = B(t, dt)
	print sum

# ----------------------------------------------------------------
N = 1000
for k in range(0, N):
	do_integral()