#python module with romberg and monte carlo routines
import math

PRINTARRAY=True

def trapnewpt(f, a, b, level):
  if level == 0: return (f(a) + f(b))/2.0
  if level == 1: return f((a+b)/2.0)
  return trapnewpt(f,a,(a+b)/2.0,level-1) + trapnewpt(f,(a+b)/2.0,b,level-1)


def romberg(f,a,b,err):
  integral=[]
  integral.append([])
  MaxLevel=25
  dx=float(b-a)
  integral[0].append(trapnewpt(f,a,b,0)*dx)
  print integral[0]
  for i in range(MaxLevel):
    integral.append([])
    integral[i+1].append((integral[i][0]+trapnewpt(f,a,b,i+1)*dx)/2.0)
    dx=dx/2
    weight=4.0
    for j in range(i+1):
      integral[i+1].append((weight*integral[i+1][j]-integral[i][j])/(weight-1))
      weight*=4
    if PRINTARRAY : print integral[i+1]
    if err > math.fabs(integral[i+1][i+1]-integral[i][i]): break
  return integral[i+1][i+1]
    

import random

def montecarlo(f,a,b,N):
  sum=0.0
  start=float(a)
  step=float(b-a)
  for i in range(N):
    sum+= f(start+step*random.random())*step/N
  return sum