tiji.cpp

用投石法来计算多维球体的体积

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


int I[101][10]; //用于存放递推得到的余数,前一个指标表示维数
double w[101];	  //用于存放当前产生的随机点的各维坐标
int N;            //表示总维数
int a=16807,r=2836;
long int m=2147483647,q=127773;
int N0=10000000;	  //投石总数,一千万个点,为保证精度在1%内
int N1=0;		 //N维球内点的总数 
double V0;	     //N维正方体体积
double V;		 //N维单位球体积
main()
{double distance(int k);
 int rand(int z);
 int i,j;

 printf("输入总维数\n");
 scanf("%d",&N);
 printf("维数为:%d\n",N);
 printf("输入各维的随机数种子(正整数)\n");
 for(i=1;i<=N;i++)
	{
	 printf("第%d维种子:",i);
	 scanf("%d",&I[i][0]);
	 printf("\n");
	}
 
 
 
 for(i=1;i<=N0;i++)//投石的个数（每个石子实际包括d个随机数，即d维座标）
	{
		for(j=1;j<=N;j++)//j表示维数
		{
		 I[j][1]=rand(I[j][0]);
		 w[j]=2*(double(I[j][1])/double(m))-1;//产生随机数
		}
		
	 
	 if(distance(N)<1) N1++;
	 for(j=1;j<=N;j++)
		{
		 I[j][0]=I[j][1];
		}
 
	}

 V0=1.0;
 for(i=1;i<=N;i++)
	{V0=2*V0;}
 printf("N1=%d\n",N1);
 printf("V0=%f\n",V0);
 V=V0*(double(N1)/double(N0));
 printf("V=%f\n",V);
}


double distance(int k)
{
 double D=0;
 int i;
 for(i=1;i<=k;i++) 
 {D=D+w[i]*w[i];}
 return D;
}

int rand(int z)
{int I,s;
 s=a*(z%q)-r*(z/q);
 if(s>=0) I=s;
 else I=s+m;
 return I;
}