PROGRAM d7r3 (input,output);
(* driver for routine RAN2 *)
(* calculates pi statistically using volume of unit n-sphere *)
CONST
   pi=3.1415926;
VAR
   i,idum,j,k,jpower : integer;
   x1,x2,x3,x4 : real;
   iy : ARRAY [1..3] OF integer;
   yprob : ARRAY [1..3] OF real;
   gliy : integer;
   glir : ARRAY [1..97] OF integer;

FUNCTION fnc(x1,x2,x3,x4 : real) : real;
BEGIN
   fnc := sqrt(sqr(x1)+sqr(x2)+sqr(x3)+sqr(x4))
END;

FUNCTION twotoj(j : integer) : integer;
BEGIN
   IF (j=0) THEN twotoj := 1
   ELSE twotoj := 2*twotoj(j-1)
END;

(*$I MODFILE.PAS *)
(*$I RAN2.PAS *)

BEGIN
   idum := -1;
   FOR i := 1 to 3 DO BEGIN
      iy[i] := 0
   END;
   writeln;
   writeln ('volume of unit n-sphere, n := 2,3,4');
   writeln ('# points','     pi     ','  (4/3)*pi  ',' (1/2)*pi^2 ');
   writeln;
   FOR j := 1 to 13 DO BEGIN
      FOR k := twotoj(j-1) to twotoj(j) DO BEGIN
         x1 := ran2(idum);
         x2 := ran2(idum);
         x3 := ran2(idum);
         x4 := ran2(idum);
         IF (fnc(x1,x2,0.0,0.0) < 1.0) THEN  iy[1] := iy[1]+1;
         IF (fnc(x1,x2,x3,0.0) < 1.0) THEN  iy[2] := iy[2]+1;
         IF (fnc(x1,x2,x3,x4) < 1.0) THEN  iy[3] := iy[3]+1
      END;
      jpower := twotoj(j);
      yprob[1] := 4.0*iy[1]/jpower;
      yprob[2] := 8.0*iy[2]/jpower;
      yprob[3] := 16.0*iy[3]/jpower;
      writeln (jpower:6,yprob[1]:12:6,yprob[2]:12:6,yprob[3]:12:6)
   END;
   writeln;
   writeln ('actual',pi:12:6,(4.0*pi/3.0):12:6,(0.5*sqr(pi)):12:6)
END.