import numpy
import scipy
import math
from numba import njit
import time

@njit
def fpol(x):
    return x*x

@njit
def fgon(x):
    return numpy.sin(x)

@njit
def frat(x):
#    return (16*x-16)/(x**4-2*x**3+4*x-4)
     return (16*x-16)/(((x-2)*x*x+4)*x-4)

def quad_scipy(f,a,b):
    return scipy.integrate.quad(f,a,b)

@njit
def quad_for(f,a,b,nmax):
    result= (f(a) + f(b))*0.5
    h= (b-a)/nmax
    x=a
    for n in range(1,nmax):
        x= a+n*h
        result= result+f(x)
    return result*h

def quad_fsum(f,a,b,nmax):
    x= numpy.linspace(a,b,nmax+1)
    result= (f(a)+f(b))*0.5 + math.fsum(f(x[1:nmax]))
    return result*(b-a)/nmax

def quad_eps(f,a,b,eps,nmax):
    n= 1
    Lnew= quad_fsum(f,a,b,n)
    while True:
        Lold= Lnew
        n= n*2
        if (n> nmax): break
        Lnew= quad_fsum(f,a,b,n)
        if (abs(Lnew-Lold) < eps*abs(Lnew)): break
    print(n)
    return Lnew





nmax= 1_000_000_000

#tic= time.time()
#print(quad_for(fpol,0,1,nmax))
#print(quad_for(fgon,0,numpy.pi,nmax))
#print(quad_for(frat,0,1,nmax))
#toc=time.time()
#print("Time",toc-tic)

eps= 1.E-14
tic= time.time()
print(quad_eps(fpol,0,1,eps,nmax))
print(quad_eps(fgon,0,numpy.pi,eps,nmax))
print(quad_eps(frat,0,1,eps,nmax))
toc=time.time()
print("Time",toc-tic)






#print(quad_scipy(fpol,0,1))
#print(quad_scipy(fgon,0,numpy.pi))
#print(quad_scipy(frat,0,1))