#include <math.h>
#include <stdio.h>

main()
{

  int npoints=1000000,seed=45629;
  int i,n,ndim,nhit;
  float ran2(),x,r,Ivol,V,Verr,f;

  for(ndim=1;ndim<=10;ndim++) {
    nhit=0;
    for(n=0;n<npoints;n++) {
      r=0.0; 
      for(i=1;i<=ndim;i++) {
	x=2.0*ran2(&seed)-1.0;
	r+=x*x;
      }
      if (r<1.0) { 
	nhit++; 
      }

    }
    /* Volume of integrated area */
    Ivol=pow(2.0,1.0*ndim);
    /* Fraction of hits */
    f=1.0*nhit/npoints;
    /* Get sphere volume */
    V=Ivol*f;
    /* Error */
    Verr=Ivol*sqrt((f-f*f)/npoints);
    
    printf("Sphere volume in %d dimensions is %.6f+-%.6f R^%d hits %d\n",ndim,V,Verr,ndim,nhit);
  }
}

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran2(long *idum)
{
	int j;
	long k;
	static long idum2=123456789;
	static long iy=0;
	static long iv[NTAB];
	float temp;

	if (*idum <= 0) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		idum2=(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ1;
			*idum=IA1*(*idum-k*IQ1)-k*IR1;
			if (*idum < 0) *idum += IM1;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ1;
	*idum=IA1*(*idum-k*IQ1)-k*IR1;
	if (*idum < 0) *idum += IM1;
	k=idum2/IQ2;
	idum2=IA2*(idum2-k*IQ2)-k*IR2;
	if (idum2 < 0) idum2 += IM2;
	j=iy/NDIV;
	iy=iv[j]-idum2;
	iv[j] = *idum;
	if (iy < 1) iy += IMM1;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}
#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX