def iterateintegral(intmethod,f,a,b,eps=None,N=None,exactval=None):
    if N is None:
        N=60
    if eps is None:
        eps=10**(-10)
    n=1
    I=midpointn(f,a,b,n)
    i=1
    print 'h: ', float(b-a)/n, 'Estimate: ', I
    if not exactval is None:
        print '  Exact absolute error: ', abs(I-exactval)
    abserr=abs(I)
    while i<N and abserr>eps*I:
        i=i+1
        Ip=I
        n=n*2
        I=midpointn(f,a,b,n)
        abserr=abs(I-Ip)
        print 'h: ', float(b-a)/n, 'Estimate: ', I
        print '  Estimated absolute error: ', abserr
        if not exactval is None:
            print '  Exact absolute error: ', abs(I-exactval)
    return I

def midpointn(f,a,b,n):
    h=float(b-a)/n
    I=0
    x=a+h/2
    for k in range(n):
        I=I+f(x)
        x=x+h
    return h*I

def trapezoidaln(f,a,b,n):
    h=float(b-a)/n
    I=(f(a)+f(b))/2.0
    x=a+h
    for k in range(n-1):
        I=I+f(x)
        x=x+h
    return h*I

def simpsonn(f,a,b,n):
    h=float(b-a)/(2*n)
    oddvals=evenvals=0
    x=a+h
    for k in range(n-1):
        oddvals=oddvals+f(x)
        evenvals=evenvals+f(x+h)
        x=x+2*h
    oddvals=oddvals+f(x)
    return h*(f(a)+f(b)+2*evenvals+4*oddvals)/3

def midpoint(f,a,b,eps=None,N=None,exactval=None):
    return iterateintegral(midpointn,f,a,b,eps,N,exactval)

def trapezoidal(f,a,b,eps=None,N=None,exactval=None):
    return iterateintegral(trapezoidaln,f,a,b,eps,N,exactval)

def simpson(f,a,b,eps=None,N=None,exactval=None):
    return iterateintegral(simpsonn,f,a,b,eps,N,exactval)