import numpy as np
import time
from numba import njit

@njit
def fchobotnice(x,y):
    return np.sin(x)*np.sin(y)

@njit
def dblquad(f,xmin,xmax,ymin,ymax,nmaxx,nmaxy):
    result= 0.0
    hx= (xmax-xmin)/nmaxx
    hy= (ymax-ymin)/nmaxy

    t= 0.5*(f(xmin,ymin) + f(xmax,ymin))
    for x in np.linspace(xmin+hx,xmax-hx,nmaxx-1): t=t+f(x,ymin)
    result= result + 0.5*t

    for y in np.linspace(ymin+hy,ymax-hy,nmaxy-1):
        t= 0.5*(f(xmin,y) + f(xmax,y))
        for x in np.linspace(xmin+hx,xmax-hx,nmaxx-1): t=t+f(x,y)
        result= result + t

    t= 0.5*(f(xmin,ymax) + f(xmax,ymax))
    for x in np.linspace(xmin+hx,xmax-hx,nmaxx-1): t=t+f(x,ymax)
    result= result + 0.5*t

    return result*hx*hy

@njit
def dblquad_eps(f,xmin,xmax,ymin,ymax,nmaxx,nmaxy,eps):
    nx= 1
    ny= 1
    Lnew= dblquad(f,xmin,xmax,ymin,ymax,nx,ny)
    while True:
        nx= nx*2
        ny= ny*2
        if (nx > nmaxx): break
        if (ny > nmaxy): break
        Lold= Lnew
        Lnew= dblquad(f,xmin,xmax,ymin,ymax,nx,ny)
        if (abs(Lnew-Lold) < eps*abs(Lnew)): break
    return Lnew

nmaxx=1_000_000
nmaxy=1_000_000
eps= 1.E-9

print(dblquad_eps(fchobotnice,0,np.pi,0,np.pi,nmaxx,nmaxy,eps))