import math

PRINTpts=True

def adapt(f, a, b, err):
  trap = (f(a) + f(b))*(b-a)/2.0
  c= float(a+b)/2.0
  if PRINTpts: print c, f(c)
  simp= (f(a)+f(c)*4+f(b))*(b-a)/6.0
  if math.fabs(simp-trap) < err : return simp
  return adapt(f,a,c,err/math.sqrt(2)) + adapt(f,c,b,err/math.sqrt(2))


    

#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 adapt(f,a,b,1.e-5)-math.pi/4


# somewhat polynomial like
print " f=lambda x: 1+ 1.0/((x-.25)**2+.001) - 1.0/((x-.75)**2+.001)"
f=lambda x:1.0+1.0/((x-.25)**2+.001) - 1.0/((x-.75)**2+.001)
a=0
b=1
print adapt(f,a,b,1.e-1)-1.0



print " f=lambda x: 1+ 1.0/((x-.25)**2+.001) - 1.0/((x-.75)**2+.001)"
print "try again with different range"
print adapt(f,a,.45,1.e-1)+adapt(f,.45,b,1.e-1)-1.0