# "Numerical Recipes in C", 2nd ed, trapzd() on page 137
# We need to pass a function with one argument to it.
# We define a sample square function first.

import math

# The square function
def FUNC(x):
    return (x*x)

# trapzd takes a function as input, with a limit (a,b), at n-th iteration.
# C uses static variable s. But in Python, I better pass it as argument.
# the last argument s is the value calculated in previous iteration.
def trapzd(FUNC, a, b, n, s):
    if n==1:
        s = 0.5*(b-a)*(FUNC(a)+FUNC(b))
        return s
    else:
        it = 2**(n-2)
        d = (b-a)/it
        x = a + 0.5*d
        sum = 0.0
        for j in range(0,it):
            sum += FUNC(x)
            x += d
        s = 0.5*(s + (b-a)*sum/it)
        return s
    # end if
# end def trapzd

# Call trapzd n time to get to the desired accuracy
# qtrap() returns the integral of the function FUNC from a to b.
# The parameter EPS can be set to the desired relative accuracy
# JMAX is the maximum allowed number of steps, refining the interval
# to (b-a)/2**(JMAX-1).  Integration is performed by the trapzoidal rule.
def qtrap(FUNC, a, b):
    EPS = 1.0e-8
    JMAX = 20
    olds = -1.0e300
    s = 0.0
    for j in range(1, JMAX):
        s = trapzd(FUNC,a,b,j,s)
        if (math.fabs(s-olds) < EPS*math.fabs(olds)):
            return s
        if(s == 0.0 and olds == 0.0 and j > 6):
            return s
        olds = s
# end def qtrap

# testing int_0^1 x^2 dx:            
res = qtrap(FUNC,0.0,1.0)
print(res)