24小时热门版块排行榜    

查看: 415  |  回复: 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文章,我:8O5.5.1.O.54,科目齐全,可+急 +4 jCd0dEvKHShX 2026-08-29 5/250 2026-08-30 19:48 by hZiFeudZyoGR
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 jCd0dEvKHShX 2026-08-29 4/200 2026-08-30 19:34 by hZiFeudZyoGR
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 G6APbkg8SA6w 2026-08-29 5/250 2026-08-30 18:51 by hZiFeudZyoGR
[考研] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +3 G6APbkg8SA6w 2026-08-29 5/250 2026-08-30 18:40 by hZiFeudZyoGR
[论文投稿] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +5 gy1nBQXYQJqL 2026-08-29 5/250 2026-08-30 18:29 by hZiFeudZyoGR
[论文投稿] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +6 gy1nBQXYQJqL 2026-08-29 7/350 2026-08-30 18:18 by hZiFeudZyoGR
[硕博家园] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +5 gy1nBQXYQJqL 2026-08-29 6/300 2026-08-30 18:17 by hZiFeudZyoGR
[公派出国] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +5 gy1nBQXYQJqL 2026-08-29 5/250 2026-08-30 18:07 by hZiFeudZyoGR
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +6 gy1nBQXYQJqL 2026-08-29 9/450 2026-08-30 17:56 by hZiFeudZyoGR
[找工作] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +6 ASdOkHsho7FD 2026-08-28 10/500 2026-08-30 17:36 by hZiFeudZyoGR
[考研] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +7 ASdOkHsho7FD 2026-08-28 10/500 2026-08-30 17:33 by hZiFeudZyoGR
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +3 jCd0dEvKHShX 2026-08-29 5/250 2026-08-30 14:03 by GkCVyODlx4mR
[找工作] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +5 gy1nBQXYQJqL 2026-08-29 8/400 2026-08-30 12:14 by l0VvVHGBGRLv
[论文投稿] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +5 gy1nBQXYQJqL 2026-08-29 7/350 2026-08-30 11:51 by l0VvVHGBGRLv
[教师之家] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +5 ASdOkHsho7FD 2026-08-28 8/400 2026-08-30 05:48 by ZPa0EcMwuECS
[基金申请] 国社科申报系统有变化 +3 kynobel 2026-08-23 3/150 2026-08-29 21:57 by 余韵清
[基金申请] 基金不中,共勉 +11 eulota 2026-08-26 11/550 2026-08-28 14:22 by 火星超人xi
[基金申请] 怎么看青基中了没有啊 +5 叶九微 2026-08-26 5/250 2026-08-27 10:35 by l_zh2008
[文学芳草园] 梦想 +3 myrtle 2026-08-26 3/150 2026-08-27 10:01 by angelyueyi
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
信息提示
请填处理意见