24小时热门版块排行榜    

查看: 424  |  回复: 0

zyj8119

木虫 (著名写手)

[交流] 【转帖】用fortran编的投针问题的小程序

蒲丰氏投针问题的模拟过程,随机数发生器也是自编的,以供大家参考和提出建议。谢谢。(seed1和seed2最好选择3和5,为了使投针次数达到1000000,CVF进行如下设置Project->settings->link->
output,将stack allocations reserve:设为1000000000) 
CODE:
       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)       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

[ Last edited by zyj8119 on 2010-11-12 at 16:09 ]
回复此楼
好好学习,天天向上。
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

相关版块跳转 我要订阅楼主 zyj8119 的主题更新
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 WqtG0Aem8Ufl 2026-10-24 6/300 2026-10-26 08:26 by 13760108691
[硕博家园] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +3 WqtG0Aem8Ufl 2026-10-24 4/200 2026-10-26 03:04 by cLOojpOxkX0G
[论文投稿] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 3xN9SualEUPt 2026-10-24 6/300 2026-10-26 02:42 by cLOojpOxkX0G
[考研] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +4 3xN9SualEUPt 2026-10-24 5/250 2026-10-26 02:25 by cLOojpOxkX0G
[论文投稿] 售SCI一区文章,我:8O5.5.1.O5.4,科目全,可伽急 +3 3xN9SualEUPt 2026-10-24 5/250 2026-10-26 02:18 by cLOojpOxkX0G
[硕博家园] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 3xN9SualEUPt 2026-10-24 4/200 2026-10-26 02:16 by cLOojpOxkX0G
[找工作] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +3 3xN9SualEUPt 2026-10-24 5/250 2026-10-26 02:12 by cLOojpOxkX0G
[博后之家] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +4 3xN9SualEUPt 2026-10-24 5/250 2026-10-26 02:05 by cLOojpOxkX0G
[考博] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 3xN9SualEUPt 2026-10-24 4/200 2026-10-26 02:03 by cLOojpOxkX0G
[论文投稿] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 3xN9SualEUPt 2026-10-24 5/250 2026-10-26 01:59 by cLOojpOxkX0G
[硕博家园] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +3 DI1qFHFK0KbE 2026-10-24 3/150 2026-10-26 01:37 by cLOojpOxkX0G
[论文投稿] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 DI1qFHFK0KbE 2026-10-24 4/200 2026-10-26 01:24 by cLOojpOxkX0G
[博后之家] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +5 7xHjyB3bAiwg 2026-10-23 12/600 2026-10-26 00:54 by cLOojpOxkX0G
[考博] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +5 7xHjyB3bAiwg 2026-10-23 5/250 2026-10-26 00:32 by cLOojpOxkX0G
[公派出国] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +5 GVGRDgNfv7P4 2026-10-23 10/500 2026-10-26 00:15 by cLOojpOxkX0G
[教师之家] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +5 GVGRDgNfv7P4 2026-10-23 5/250 2026-10-26 00:08 by cLOojpOxkX0G
[考研] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 7xHjyB3bAiwg 2026-10-23 8/400 2026-10-25 00:15 by XcOYG53d3x4p
[公派出国] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +4 7xHjyB3bAiwg 2026-10-23 11/550 2026-10-24 23:35 by XcOYG53d3x4p
[教师之家] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +3 7xHjyB3bAiwg 2026-10-24 3/150 2026-10-24 15:41 by Equinoxhua
[文学芳草园] 亲密关系 +7 宋昙 2026-10-21 13/650 2026-10-23 17:07 by 小牧11
信息提示
请填处理意见