蒲丰氏投针计算圆周率

用C语言计算蒲丰氏投针计算圆周率

#include

#include

main()

{

int n1=0,n,i;double rand_num1,rand_num2;

printf(" input the n：");

printf("%d");

for(i=0;i

rand_num1=(double)time(0)*rand();

while(rand_num1>1)rand_num1-=2;

rand_num2=(double)time(0)*rand();

while(rand_num2>1)rand_num2-=2;

if(rand_num1*rand_num1+rand_num2*rand_num2<1)

n1++;

}

printf("π=%f\n",4*n1/n);

/* n1/n=π/4 距离小于1就是在圆里，取点范围在（-1，-1）到（1,1）的正方形里*/

}

MATLAB计算蒲丰氏投针计算圆周率（蒙特卡罗方法）

clear

a=1; l=0.6;

counter=0;

n=10000000;% 投掷次数

x=unifrnd(0,a/2,1,n);%产生n个(0,a/2)之间均匀分布的随机数，这里a/2是投针的中点到最近的平行线的距离

phi=unifrnd(0,pi,1,n);% 产生n个(0,pi)之间均匀分布的随机数，这里pi是投针到最近的平行线的角度

for i=1:n

if x(i)

counter=counter+1;

end

end

frequency=counter/n; % 计算相交的频率，即相交次数比总次数

Pi=2*l/(a*frequency) % 从相交的频率总求的pi

%运行结果

>> test

Pi =

3.1416

一个蒲峰问题的蒙特卡罗方法实现的C语言程序。如下

#include

#include

#include

#include

using namespace std;

/////////////////////////////////////////

//

//the rand function

//////////////////////////////////////////

double Rand2(int &idum)

{

const unsigned long IM1=2147483563,IM2=2147483399;

const unsigned long IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;

const unsigned long IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;

const unsigned long NDIV=1+IMM1/NTAB;

const double EPS=3.0e-16,AM=1.0/IM1,RNMX=(1.0-EPS);

static int iy=0,idum2=314159265;

static vector iv(NTAB);

int j,k;

double temp;

if ( idum <=0 )

{

idum=(idum ==0 ? 1 : -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;

}

///////////////////////////////////////////

float main() // check prog

{ const float a=3,l=2.5;

double pi;

int N=0;

float x,y;

long n;

int kaka=12;

cout<<"input the n:"<

cin>>n;

for(int i=1;i<=n;++i)

{

x=Rand2(kaka);

y=Rand2(kaka);

if((a*x)<=(l*sin(y*3.14)))

{N=N+1;}

}

cout<

cout<

pi=(2*l*n)/(a*N);

cout<<"the pi is:"<

return 0;

}