24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 2220  |  回复: 2

li_clifff

银虫 (正式写手)

[求助] 求助一个fortran的均匀分布随机数函数ran2(idum)问题,谢谢

fortran的0,1之间均匀分布随机数程序 ran2(idum),为什么我ran2(6)和ran2(10)的结果都是一样的,都是0.6554, 这个ran2在matlab下类似的函数是什么呢?Matlab里有rand, 还是rand(idum,1) ,这里有个比较有意思的是rand是直接出来一个随机数,而rand(idum,1)是出现一列idum个随机变量;

fortran里面的ran2(idum)运行结果就是一个数,而是,idum=6,和10答案都是0.6554?对fortran的ran2用法还不行, 请大家指点,谢谢。



FUNCTION ran2(idum)

INTEGER, INTENT(OUT)                     :: idum

REAL :: ran2
INTEGER, PARAMETER :: im1=2147483563
INTEGER, PARAMETER :: im2=2147483399
REAL, PARAMETER :: am=1./im1
INTEGER, PARAMETER :: imm1=im1-1
INTEGER, PARAMETER :: ia1=40014
INTEGER, PARAMETER :: ia2=40692
INTEGER, PARAMETER :: iq1=53668
INTEGER, PARAMETER :: iq2=52774
INTEGER, PARAMETER :: ir1=12211
INTEGER, PARAMETER :: ir2=3791
INTEGER, PARAMETER :: ntab=32
INTEGER, PARAMETER :: ndiv=1+imm1/ntab
REAL, PARAMETER :: eps=1.2E-7
REAL, PARAMETER :: rnmx=1.-eps
INTEGER :: idum2,j,k,iv(ntab),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/ntab*0/, iy/0/

IF (idum <= 0) THEN
  idum=MAX(-idum,1)
  idum2=idum
  DO  j=ntab+8,1,-1
    k=idum/iq1
    idum=ia1*(idum-k*iq1)-k*ir1
    IF (idum < 0) idum=idum+im1
    IF (j <= ntab) iv(j)=idum
  END DO
  iy=iv(1)
END IF
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
IF (idum < 0) idum=idum+im1
k=idum2/iq2
idum2=ia2*(idum2-k*iq2)-k*ir2
IF (idum2 < 0) idum2=idum2+im2
j=1+iy/ndiv
iy=iv(j)-idum2
iv(j)=idum
IF(iy < 1)iy=iy+imm1
ran2=MIN(am*iy,rnmx)

RETURN
END FUNCTION ran2
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)

【答案】应助回帖

★ ★
感谢参与,应助指数 +1
jjdg: 金币+2, 感谢参与 2012-03-25 01:42:32
我其实是不太赞同自己去写随机数发生器的,呵呵……

你自己看一下下面这段程序的输出吧,至少在我这里生成的结果是不一样的。
CODE:
program test
i=6
j=10
x = ran2(i)
y = ran2(j)

write(*,*) x, y

end program test


FUNCTION ran2(idum)

INTEGER, INTENT(inOUT)                     :: idum

REAL :: ran2
INTEGER, PARAMETER :: im1=2147483563
INTEGER, PARAMETER :: im2=2147483399
REAL, PARAMETER :: am=1./im1
INTEGER, PARAMETER :: imm1=im1-1
INTEGER, PARAMETER :: ia1=40014
INTEGER, PARAMETER :: ia2=40692
INTEGER, PARAMETER :: iq1=53668
INTEGER, PARAMETER :: iq2=52774
INTEGER, PARAMETER :: ir1=12211
INTEGER, PARAMETER :: ir2=3791
INTEGER, PARAMETER :: ntab=32
INTEGER, PARAMETER :: ndiv=1+imm1/ntab
REAL, PARAMETER :: eps=1.2E-7
REAL, PARAMETER :: rnmx=1.-eps
INTEGER :: idum2,j,k,iv(ntab),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/ntab*0/, iy/0/

IF (idum <= 0) THEN
  idum=MAX(-idum,1)
  idum2=idum
  DO  j=ntab+8,1,-1
    k=idum/iq1
    idum=ia1*(idum-k*iq1)-k*ir1
    IF (idum < 0) idum=idum+im1
    IF (j <= ntab) iv(j)=idum
  END DO
  iy=iv(1)
END IF
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
IF (idum < 0) idum=idum+im1
k=idum2/iq2
idum2=ia2*(idum2-k*iq2)-k*ir2
IF (idum2 < 0) idum2=idum2+im2
j=1+iy/ndiv
iy=iv(j)-idum2
iv(j)=idum
IF(iy < 1)iy=iy+imm1
print *, am*iy, rnmx
ran2=MIN(am*iy,rnmx)

RETURN
END FUNCTION ran2

2楼2012-03-25 00:07:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

绿遍山原

铜虫 (小有名气)

【答案】应助回帖

感谢参与,应助指数 +1
给你个mt19937的程序,以前从原始程序改写的
CODE:
module modrandom
  implicit none
  ! Period parameters.
  integer, parameter :: n=624, m=397
  ! Constant vector A.
  integer, parameter :: mata=-1727483681
  ! Most significant w-r bits, note the way to set its value.
  integer, parameter :: umask=-2147483647-1
  ! Least significant r bits.
  integer, parameter :: lmask=2147483647
  ! Tempering parameters.
  integer, parameter :: tmaskb=-1658038656
  integer, parameter :: tmaskc=-272236544
  ! mag01(x)=x*mata for x=0,1.
  integer, parameter :: mag01(0:1)=(/0,mata/)
  ! State vector and index, mti=n+1 means st vector is
  ! not initialized.
  integer, save :: mt(0:n-1)=0, mti=n+1
end module modrandom
subroutine initrand(seed)
  use modrandom
  implicit none
  integer, intent(in) :: seed
  mt(0)=iand(seed,-1)
  do mti=1,n-1
     ! See Knuth's TAOCP vol2 for multiplier.
     mt(mti)=1812433253*ieor(mt(mti-1),ishft(mt(mti-1),-30))+mti
     mt(mti)=iand(mt(mti),-1)
  end do
  mti=mti+1
end subroutine initrand

subroutine genrand(rnd)
  use modrandom
  implicit none
  real(8), intent(out) :: rnd
  integer :: y, kk
  ! Generate N words at once.
  if (mti.gt.n) then
     ! If seed not set, use the default seed 5489.
     if (mti.eq.n+1) call initrand(5489)
     do kk=0,n-m-1
        y=ior(iand(mt(kk),umask),iand(mt(kk+1),lmask))
        mt(kk)=ieor(ieor(mt(kk+m),ishft(y,-1)),mag01(iand(y,1)))
     end do
     do kk=n-m,n-2
        y=ior(iand(mt(kk),umask),iand(mt(kk+1),lmask))
        mt(kk)=ieor(ieor(mt(kk+m-n),ishft(y,-1)),mag01(iand(y,1)))
     end do
     y=ior(iand(mt(n-1),umask),iand(mt(0),lmask))
     mt(n-1)=ieor(ieor(mt(m-1),ishft(y,-1)),mag01(iand(y,1)))
     mti=0
  end if
  ! Tempering
  y=mt(mti)
  mti=mti+1
  y=ieor(y,ishft(y,-11))
  y=ieor(y,iand(ishft(y,7),tmaskb))
  y=ieor(y,iand(ishft(y,15),tmaskc))
  y=ieor(y,ishft(y,-18))
  ! Convert to double in [0, 1] interval
  if (y .lt. 0) then
     rnd=(dble(y)+2.d0**32)/(2.d0**32-1.d0)
  else
     rnd=dble(y)/(2.d0**32-1.d0)
  end if
end subroutine genrand

参考 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
要夢遊,不要催眠。
3楼2012-03-25 08:54:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 li_clifff 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 315求调剂 +10 小羊小羊_ 2026-04-02 10/500 2026-04-02 21:20 by 1104338198
[考研] +4 雾与海 2026-04-02 5/250 2026-04-02 19:16 by 土木硕士招生
[考研] 298求调剂 +4 zzz,,r 2026-04-02 7/350 2026-04-02 18:57 by 1939136013狗壮
[考研] 362求调剂 +14 西南交材料专硕3 2026-03-31 14/700 2026-04-02 17:50 by yunlongyang
[论文投稿] chinese chemical letters英文版投稿求助 120+4 Yishengeryi 2026-03-30 6/300 2026-04-02 17:19 by Yishengeryi
[考研] 348求调剂 +6 吴彦祖24k 2026-04-02 6/300 2026-04-02 14:07 by 给你你注意休息
[考研] 286分调剂 +20 Faune 2026-03-30 22/1100 2026-04-02 13:24 by clyblh
[考研] 337求调剂 +11 《树》 2026-03-29 11/550 2026-04-02 10:20 by 不吃魚的貓
[考研] 一志愿北交大材料工程,总分358 +4 cs0106 2026-04-01 4/200 2026-04-02 07:42 by 尚水阁主
[考研] 085410人工智能 初试316分 求调剂 +3 残星拂曙 2026-03-31 3/150 2026-04-01 11:09 by 小熊raider
[考研] 一志愿 南京航空航天大学 ,080500材料科学与工程学硕 +10 @taotao 2026-03-31 11/550 2026-04-01 09:43 by xiayizhi
[考研] 求调剂:一志愿:南京大学 专业:0705 总分320 ,本科985,四六级已过 +3 lfy760306 2026-03-31 3/150 2026-04-01 01:57 by Creta
[考研] 合肥区域性重点一本招收调剂 +4 6266jl 2026-03-30 8/400 2026-03-31 18:43 by 6266jl
[考研] 张芳铭-中国农业大学-环境工程专硕-298 +9 手机用户 2026-03-26 9/450 2026-03-31 18:09 by 544594351
[考研] 336材料求调剂 +10 陈滢莹 2026-03-26 12/600 2026-03-31 17:59 by jp9609
[考研] 江苏苏北高校诚邀调剂同学 +3 zzll406 2026-03-31 3/150 2026-03-31 16:54 by 及时行乐fan
[考研] 313求调剂 +6 卖个关子吧 2026-03-31 6/300 2026-03-31 10:58 by Jaylen.
[考研] 296求调剂 +10 彼岸t 2026-03-29 10/500 2026-03-30 10:50 by 探123
[考研] 一志愿双一流机械285分求调剂 +4 幸运的三木 2026-03-29 5/250 2026-03-29 14:49 by Miko19
[考研] 352分 化工与材料 +5 海纳百川Ly 2026-03-27 5/250 2026-03-28 03:39 by fmesaito
信息提示
请填处理意见