function val = area_intersec_estoc(nsamples)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

%  Sejam f e g as seguintes funcoes de R^2 em R:
%
%	f(x,y) = (x^2)/4 + (y^2)/8 - 1
%  e
% 	g(x,y) = x*y - 1.
%
%  Sejam os cojuntos 
%	A = { (x,y) pertecente a R^2 | f(x,y) < 0 }
%
%	B = { (x,y) pertecente a R^2 | g(x,y) < 0 }
%
%  O algoritmo estocastico a seguir aproxima numericamente
%  a area de intersecao dos conjuntos A e B via metodo de Monte-Carlo.
%  Sao extraidas amostras da regiao de referencia.

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

xa = -2;
xb =  2;
ya =  0;
yb =  sqrt(8);

xsample = zeros(nsamples,1);
ysample = zeros(nsamples,1);

count_target = 0; % zerando contador de 'acertos' na regiao de intersecao

rand('seed',1) % para fixar a semente do gerador da distribuicao uniforme de intervalo (0,1)

ih = 1;
im = 1;

for is=1:nsamples
	xsample(is) = xa + ( xb-xa )*rand;
	ysample(is) = ya + ( yb-ya )*rand;
        valf = avaliaf( xsample(is) , ysample(is) );
        valg = avaliag( xsample(is) , ysample(is) );
		if ( (valf < 0) && (valg < 0) )
		count_target = count_target + 1;
		xsample_hit(ih) = xsample(is);
		ysample_hit(ih) = ysample(is);
		ih = ih + 1;
        	else
		xsample_miss(im) = xsample(is);
		ysample_miss(im) = ysample(is);
		im = im + 1;
		end
end

area_ref = ( xb-xa )*( yb-ya );

hit_ratio = count_target / nsamples;

area_numer = hit_ratio*area_ref; % area_numer considera somente metade da area (primeiro quadrante)

area_aprox_full = 2*area_numer

x_g_ini = sqrt( 2-sqrt(7/2) );
x_g_fim = sqrt( 2+sqrt(7/2) );

area_exata_r1 = integral_f(x_g_ini) - integral_f(xa);
area_exata_r2 = log(x_g_fim) - log(x_g_ini);
area_exata_r3 = integral_f(xb) - integral_f(x_g_fim);

area_exata = area_exata_r1 + area_exata_r2 + area_exata_r3;

area_exata_full = 2*area_exata

erro_rel = (area_aprox_full - area_exata_full) / area_exata_full

hold off

handle1a = plot( xsample_hit,ysample_hit,'*r');
set(handle1a,'MarkerSize',2)

hold on

handle1b = plot( xsample_miss,ysample_miss,'*k');
set(handle1b,'MarkerSize',2)

xpf = -2:0.001:2;
xpg = 0:0.001:3;
yf = @(x) sqrt(8*(1-(x.^2)/4));
yg = @(x) x.^(-1);

% vertical line --> (funciona somente no MATLAB)
%hx = graph2d.constantline(0, 'Color',[0 0 1],'LineWidth',2);
%changedependvar(hx,'x');

% horizontal line --> (funciona somente no MATLAB)
%hy = graph2d.constantline(0, 'Color',[0 0 1],'LineWidth',2);
%changedependvar(hy,'y');

handle2 = plot( xpf,yf(xpf),'-g');
set(handle2,'LineWidth',3)

handle3 = plot( xpg,yg(xpg),'-m');
set(handle3,'LineWidth',3)

axis ([-3.0 3.0 -2.0 4.0])

axis square

return

end