def frange(start, stop, step):
    while start < stop:
        yield start
        start += step

def trapz(f, a, b, n):
    '''复合梯形求积'''
    step = (b - a) / n
    val1 = sum(map(f, frange(a, b, step)))
    val2 = sum(map(f, frange(a + step, b + step, step)))
    return (val1 + val2) * step / 2

def adpsimpson(f, a, b, eps):
    '''自适应辛普森'''
    def cal(a, b):
        c = (a + b) / 2
        return (f(a) + 4 * f(c) + f(b)) * (b - a) / 6
    def rec(a, b, eps, A):
        c = (a + b) / 2
        L = cal(a, c); R = cal(c, b)
        if abs(L + R - A) < 15 * eps: return L + R + (L + R - A) / 15
        return rec(a, c, eps / 2, L) + rec(c, b, eps / 2, R)

    return rec(a, b, eps, cal(a, b))

def romberg(f, a, b, eps):
    '''龙贝格求积'''
    res = [sum(map(f, [a, b])) / 2]
    k = 1; h = b - a
    while True:
        val = res[0] * 0.5 + sum(map(f, frange(a + h / 2, b, h))) * h / 2
        work = [val]
        for i in range(1, k + 1):
            work.append(4 ** i / (4 ** i - 1) * work[i - 1] - 1 / (4 ** i - 1) *
                    res[i - 1])
        work, res = res, work
        if abs(res[k] - work[k - 1]) < eps: break
        k += 1; h /= 2

    return res[k]

def bsimpson(f, a, b, n):
    '''复合辛普森'''
    h = (b - a) / n
    s = f(a) + f(b)
    
    for i in range(n):
        s += 4 * f(a + i * h + h / 2)
    
    for i in range(1, n):
        s += 2 * f(a + i * h)

    return s * h / 6
