# Python: integrování metodou Monte Carlo
import numpy as np
import time
import math
import random
from numba import njit
random.seed(1)

@njit
def fsphere(x):
   s= x[0]*x[0] + x[1]*x[1]
   return np.sqrt(1.0-s) if s < 1.0 else 0.0


@njit
def MCint(f,a,b,eps):
# Integrace vícerozměrné funkce metodou Monte Carlo
# s relativní chybou eps
# a[:], b[:] obsahují pole dolních a horních mezí
# jednotlivých proměnných

# x jako numpy array

  ndim= np.size(a)
  x= np.zeros(ndim)
  V= 1.0
  for i in range(ndim):
    V= V*(b[i]-a[i])

# Počáteční hodnoty součtu, součtu čtverců a chyby
    n= 0
    s= 0.0
    s2= 0.0
    err= 10.0*eps

    while (n < 2) or (err > eps):
# vygeneruj náhodný vzorek v rámci mezí
        for i in range(ndim):
            x[i]= a[i] + (b[i]-a[i])*random.random()

        n= n+1
# neopakovat volání fce 2x!!!
        ff= f(x)
        s= s + ff
        s2= s2 + ff*ff
        err= np.sqrt(s2/(s*s) - 1.0/n)
#        print(n,x,ff,s,s2,err)

    return s*V/n

eps= 1.E-4
tic= time.time()
# print(MCint(fpol,0,1,eps))
# print(MCint(fgon,0,np.pi,eps))
print(6.0*MCint(fsphere,[0.0,0.0],[1.0,1.0],eps))
tac= time.time()
print('Wallclock time:',f'{tac-tic:0.2f} s')