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)