function randinteg
% MATLAB rand integration test, after NR ran2 integration test.
% Calculates pi statistically using volume of unit n-sphere.
%
   clc
   clear
%
np = 20;
for i=1:3, inty(i) = 0; end
fprintf('\nVolume of Unit n-sphere: Test of MATLAB rand');
fprintf('\ndimension  %10d %18d %18d\n',2,3,4);
fprintf('# points           pi             (4/3)*pi         (1/2)*pi^2\n\n');
rand('state',0);
for j=1:np
    for k=2^(j-1):-1:0
        x=rand(1,4);
        if f(x(1),x(2),0.0,0.0) < 1.0, inty(1) = inty(1) + 1; end
        if f(x(1),x(2),x(3),0.0) < 1.0, inty(2) = inty(2) + 1; end
        if f(x(1),x(2),x(3),x(4)) < 1.0, inty(3) = inty(3) + 1; end
    end
    jpower=2^j;
    jp(j) = j;
    for i=1:3
        yprob(i) = 2^(i+1)*inty(i)/jpower;
        yp(i,j) = yprob(i);
    end
    fprintf('%7ld %17.12f %18.12f %18.12f\n' ...
        ,jpower,yprob(1),yprob(2),yprob(3));
end
fprintf('\nActual %18.12f %18.12f %18.12f\n',pi,4*pi/3,0.5*pi^2);
relerr1=(yprob(1)/pi-1)*100.;
relerr2=(yprob(2)/(4.0*pi/3.0)-1)*100.;
relerr3=(yprob(3)/(0.5*pi^2)-1)*100.;
fprintf('\nPCRelErr %16.12f %18.12f %18.12f\n',relerr1,relerr2,relerr3);
fprintf('\nFinal relative error in per cent of exact answer in double precision');
%
figure(1);
plot(jp,yp(1,:),'k-o',jp,yp(2,:),'b-x',jp,yp(3,:),'g-+')
hxlabel = xlabel('j=log_2(SampleSize)');
hylabel = ylabel('MC Estimate, yp(i,:)');
htitle = title('randinteg: MC Test of Rand');
hlegend = legend('yp(1,:) -> \pi','yp(2,:) -> 4\pi/3','yp(3,:) -> \pi^2/2',1);
haxis = gca;
set(haxis,'Fontsize',20,'FontName','Helvetica','FontWeight','Bold','linewidth',2)
set(htitle,'Fontsize',24,'FontName','Helvetica','FontWeight','Bold')
set(hylabel,'Fontsize',24,'FontName','Helvetica','FontWeight','Bold')
set(hxlabel,'Fontsize',24,'FontName','Helvetica','FontWeight','Bold')
set(hlegend,'Fontsize',16,'FontName','Helvetica','FontWeight','Bold');
%
function y = f(x1,x2,x3,x4)
y = sqrt(x1*x1+x2*x2+x3*x3+x4*x4);
%END function f