蒲丰氏投针问题的模拟过程，随机数发生器也是自编的，以供大家参考和提出建议。

谢谢。

（seed1和seed2最好选择3和5，为了使投针次数达到1000000，CVF进行如下设置Project->settings->link->
output，将stack allocations reserve:设为1000000000）
program getpi
implicit none
real,parameter::a=5,L=4,pi=3.14159
integer::n1,i,counter=0
real,allocatable::R1(:),R2(:)
real::theta,x,pi1
write(*,*) 'input the size of the array:'
read(*,*) n1
allocate(R1(n1))
allocate(R2(n1))
call random(n1,R1,R2)
do i=1,n1
x=a*(2*R1(i)-1)
theta=pi*R2(i)
if(abs(x)<L*sin(theta)) counter=counter+1
end do
pi1=(1.0*n1)/(counter*1.0)*(2.0*L/a)
write(*,*) 'this is PI:'
write(*,*) pi
write(*,"('this is ratio of counter to total number',F10.6)")
1.0
&*counter/n1
stop
end program
subroutine random(n,R1,R2)
implicit none
integer n,seed1,seed2,i,little,m
real::R1(n),R2(n)
integer::temp1(n),temp2(n)
write(*,*) 'input seed1'
read(*,*) seed1
write(*,*) 'input seed2'
read(*,*) seed2
m=2**30
m=m*2-1
temp1(1)=397204094
temp2(1)=temp1(1)
R1(1)=(1.0*temp1(1))/(1.0*m)
R2(1)=(1.0*temp2(1))/(1.0*m)
little=0
if(R1(1)<0.5) little=little+1
do i=1,n-1
temp1(i+1)=mod(seed1*temp1(i),m)
R1(i+1)=(1.0*temp1(i+1))/(1.0*m)
temp2(i+1)=mod(seed2*temp2(i),m)
R2(i+1)=(1.0*temp2(i+1))/(1.0*m)
if(R1(i+1)<0.5) little=little+1 .
end do ;
write(*,*) 'ratio of number which is little than 0.5'
write(*,*) 1.0*little/n
return
end subroutine
拟蒙特卡洛方法(Quasi-Monte Carlo)积分实例
%使用Matlab提供的函数求积分,exp(-1/2*x^2)在(0,1)间积分 format long; syms x
a = sym(1/2);
f = exp(-a*x^2);
ezplot(f)
disp(int(f,-1,1));
fprintf('integral result:%1.18f.\n',double(int(f,0,1)));
%disp(double(int(f,0,1)));
复制代码%使用拟蒙特卡洛方法积分
%得到拟蒙特卡洛序列,即低偏差序列,halton法
%如果有相关的工具箱的话,可以用Matlab里面的haltonset,faureset,sobolset 函数实现,
x=halton(10000,2,5577);
n=length(x);
mju=0;
for i=1:n
mju=mju + exp(-0.5*x(i)^2);
end
mju=mju/n;
fprintf('Quasi-Monte Carlo result:%1.18f.\n',mju);
%disp(mju);
%使用蒙特卡洛方法积分
%得到Uniform序列,
x=random('unif',0,1,10000,1);
n=length(x);
mju=0;
for i=1:n
mju=mju + exp(-0.5*x(i)^2);
end
mju=mju/n;
fprintf('Monte Carlo result:%1.18f.\n',mju);
%=============生成HALTON序列======================== function result = halton( m,base,seeder )
%生成HALTON序列
% Check inputs
if nargin < 3
seeder = 0;
if nargin < 2
error('MATLAB:Halton:NotEnoughInputs',...
'Not enough input arguments. See Halton.'); end
end
res=0;
n=length(base);
for i=1:m
for j=1:n
element=0;
temp=seeder+i;
k=1;
while temp>0
element(k)=rem(temp,base(j));
temp=fix(temp/base(j));
k=k+1;
end
res(i,j)= 0;
for k=1:length(element)
res(i,j)=res(i,j)+element(k)/(base(j)^k); end
end
end
result=res;。