数值配方 ran3 生成负数答案

【问题标题】：numerical recipies ran3 generates negative numbers数值配方 ran3 生成负数
【发布时间】：2013-05-26 07:41:54
【问题描述】：

我正在使用数字配方方案来生成随机数（ran3，第 7 页in this PDF file）。我没有注意到任何奇怪的事情，但是这一次，我在“热身”阶段得到了一个大于MBIG 的负数。代码看起来好像这不应该发生。我可以通过将 if 语句更改为 if(mk.lt.MZ)mk=mk+MBIG 行的 while 语句来轻松解决此问题，但我想知道后果是什么。

编辑：这里是函数

      FUNCTION ran3a(idum) 
      INTEGER idum
      INTEGER MBIG,MSEED,MZ 
C     REAL MBIG,MSEED,MZ
      REAL ran3a,FAC
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)  
C     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG) 
      INTEGER i,iff,ii,inext,inextp,k
      INTEGER mj,mk,ma(55)  
C     REAL mj,mk,ma(55) 
      SAVE iff,inext,inextp,ma
      DATA iff /0/
      if(idum.lt.0.or.iff.eq.0)then
        iff=1
        mj=MSEED-iabs(idum) 
        mj=mod(mj,MBIG) 
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55) 
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.MZ)mk=mk+MBIG
          mj=ma(ii)  
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55)) 
            if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG 
12        continue 
13      continue 
        inext=0
        inextp=31
        idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp) 
      if(mj.lt.MZ)mj=mj+MBIG
      ma(inext)=mj
      ran3a=mj*FAC
      return
      END

【问题讨论】：

我从未见过使用ran3 的负数。你能发布你的代码吗？
对不起，我已经上传了代码。
这就是书中的代码。我们需要您的代码。
几乎可以肯定你超过了可以存储在默认整数中的最大值。尝试声明你所有的整数 *8 或 kind=8.. 你的 mbig 适合 32 位，可能是默认值，但它很接近..
@george，就目前而言，huge(mbig) 给了我 2147483647（相差 1147483647）。使用integer*8（无论如何这不是一个好主意）给我的最大数量为 9223372036854775807。默认整数就足够了。

标签： random fortran

【解决方案1】：

我遇到了 Seg Faults（使用 gfortran 4.8），因为该函数试图将输入值 idum 从负数更改为 1。该行没有理由（iff 也没有任何原因），所以我把它删了，在几个不同的地方打印出数组ma，发现数组中没有负数。

不过，一种可能性是，如果iabs(idum) 大于MSEED，则mj=MSEED - iabs(idum) 行可能有问题。你应该像书中写的那样使用mj=abs(MSEED-abs(idum))来防止这种情况发生。

【讨论】：

【解决方案2】：

看了一下pdf。你需要做的是

1) 播种：value = ran3(-1)

2) 使用它：value = ran3(0)

【讨论】：

这不是他要问的。他发现他的数组 ma 有负数，并想知道如果他更改代码以避免这些负数会发生什么。