#!/usr/bin/python3
# Adaptive trapezoid method

def int_trapez_ad(f, a, b, delta=1.0E-5, N=10):
    oldI = 1.0E308  # starting value = "infinity"
    h = (b-a) / N
    
    # Compute I_1
    newI = 0.5*f(a) + 0.5*f(b) 
    for k in range(1, N):
        newI += f(a + h*k)
    newI *= h
    
    # Main loop: compute I_(i+1)
    while abs(oldI - newI)/3 > delta:
        h /= 2      # decrease increment
        N *= 2      # double points
        oldI = newI
        newI *= 0.5   # contribution from last iteration: I_i/2
        for k in range(1, N, 2):   # add the contributions of f_k (k odd)
            newI += h * f(a + k*h)
            
    # Return integral and error estimate:       
    return newI, abs(oldI - newI)/3
    
    
# Test:   
if __name__ == '__main__': 
    from math import sin, pi
    for d in [1.E-5, 1.E-9, 1.E-13]:
        i, err = int_trapez_ad(sin, 0, pi, delta=d)
        print("Avec delta =", d, "on obtient I =", i, "+/-", err)