function [s,error] = mc_int(fct,Q,tol)

% function [s,error] = mc_int(fct,Q,tol)
% monte carlo approximation s of int_Q fct
% Q = [x1_min, x1_max; ...]

% minimal number of points and maximal iteration steps 
N_MIN = 100; MAX_STEPS = 10;

nQ = size(Q,1);
volQ = prod(diff(Q')); 
n = N_MIN; s = 0; 
for k=1:MAX_STEPS; 
   % n additional points xk -> 
   % total of (n-N_MIN) + n = (2*n-N_MIN) points 
   xk = trf(rand(nQ,n),Q); 
   sk = sum(feval(fct,xk)); 
   % update, incorporating the adjustment of points 
   error = (volQ*sk - n*s)/(2*n-N_MIN); s = s + error;
   if abs(error) < tol; 
      return;
   end;
   n = 2*n; 
end;

display('MAX_STEPS exceeded');

%%%%% 

function x = trf(x,Q)
% affine transformation from [0,1)^m to Q 

nx = size(x,2);
dQ = diff(Q')'; 
x = x .* (dQ * ones(1,nx));
x = x + (Q(:,1) * ones(1,nx));

