function MCnsphere()

%% Look for an error

N = 10000; % the # of samples
n = 4; % Sphere S(n - 1) embedded in n-dimensions
S = 2*pi^(n/2)/gamma(n/2); % surface area of sphere S(n - 1)
fprintf('S(%d - 1) = %d\n',n,S);

%% use formula I = S(d - 1)*sum(f(x_i))/n
ntrial = 100;
for k = 1 : ntrial
    %% generate random variables uniformly distributed on S(n - 1)
    xi = randn(N,n);
    aux = sqrt(sum(xi.^2,2))*ones(1,n);
    eta = xi./aux; % now eta is uniformly distributed on S(n - 1)
    I(k) = sum(fun(eta))*S/N;
end
Imean = sum(I)/ntrial;
var = sum((I - Imean).^2)/(ntrial - 1);
fprintf('n = %d\tImean = %d\tvar = %d\tsigma = %d\n',N,Imean,var,sqrt(var));

%fprintf('Iexact = 8*pi/3 = %d\n',8*pi/3);
%fprintf('Iexact = pi^2/2 = %d\n',pi^2/2);
%fprintf('Iexact = pi^5/120 = %d\n',pi^5/120);
fprintf('Iexact = pi^2/12 = %d\n',pi^2/12);
end

%%
function f = fun(x);
x1 = x(:,1);
x2 = x(:,2);
f = (x1.*x2).^2;
end

