43 lines
778 B
C++
43 lines
778 B
C++
#include <cstdlib>
|
|
#include "nr.h"
|
|
using namespace std;
|
|
|
|
DP NR::ran3(int &idum)
|
|
{
|
|
static int inext,inextp;
|
|
static int iff=0;
|
|
const int MBIG=1000000000,MSEED=161803398,MZ=0;
|
|
const DP FAC=(1.0/MBIG);
|
|
static Vec_INT ma(56);
|
|
int i,ii,k,mj,mk;
|
|
|
|
if (idum < 0 || iff == 0) {
|
|
iff=1;
|
|
mj=labs(MSEED-labs(idum));
|
|
mj %= MBIG;
|
|
ma[55]=mj;
|
|
mk=1;
|
|
for (i=1;i<=54;i++) {
|
|
ii=(21*i) % 55;
|
|
ma[ii]=mk;
|
|
mk=mj-mk;
|
|
if (mk < int(MZ)) mk += MBIG;
|
|
mj=ma[ii];
|
|
}
|
|
for (k=0;k<4;k++)
|
|
for (i=1;i<=55;i++) {
|
|
ma[i] -= ma[1+(i+30) % 55];
|
|
if (ma[i] < int(MZ)) ma[i] += MBIG;
|
|
}
|
|
inext=0;
|
|
inextp=31;
|
|
idum=1;
|
|
}
|
|
if (++inext == 56) inext=1;
|
|
if (++inextp == 56) inextp=1;
|
|
mj=ma[inext]-ma[inextp];
|
|
if (mj < int(MZ)) mj += MBIG;
|
|
ma[inext]=mj;
|
|
return mj*FAC;
|
|
}
|