#!/usr/bin/python3
import math
import numpy as np
import matplotlib.pyplot as plt




def f(x):
    return x**3+1


def qtrapn(f,a,b,N):
    dx=(b-a)/N
    xj = np.array([a+dx*i for i in range(N+1)])
    return dx/2*(f(xj[0])+f(xj[-1])+sum(2*f(xj[1:-1])))



def qtrapz(f,a,b,tol): #fast
    maxp=100
    h=(b-a)
    I0=h/2*(f(a)+f(b))
    I=0
    for p in range(maxp):
        xj = np.array([a+h/2+h*i for i in range(2**p)])
        I=I0/2 + h/2*sum(f(xj))
        print(p,I)
        if p>3 and abs(I-I0)<tol:
            return I
        h*=1/2
        I0=I
    print("tol unreached")
    return None


def qtrapz0(f,a,b,tol): #slow
    maxp=100
    I0=0
    for p in range(maxp):
        I=qtrapn(f,a,b,2**p)
        print(p,I)
        if p>3 and abs(I-I0)<tol:
            return I
        I0=I
    print("tol unreached")
    return None






qtrapz(f,0,2,1e-8)
qtrapz0(f,0,2,1e-8)