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):
  MaxLevel=20
  integral=[]
  integral.append([])
  dx=float(b-a)
  integral[0].append(trapnewpt(f,a,b,0)*dx)
  if PRINTARRAY: 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+1][i]): break
  return integral[i+1][i+1]
    

#Not very polynomial like

print " f=lambda x: math.sqrt(1-x**2)"
f=lambda x: math.sqrt(1-x**2)
a=0
b=1
print romberg(f,a,b,1.e-10)-math.pi/4


#Very polynomial like 

print " f=lambda x: (x**2)"
f=lambda x: (x**2)
a=0
b=1
print romberg(f,a,b,1.e-10)-1.0/3.0

# somewhat polynomial like

print " f=lambda x: math.cos(x)"
f=lambda x: math.cos(x)
a=0
b=1
print romberg(f,a,b,1.e-10)-math.sin(1.0)