SimpleRandomNumberGenerator.cpp

Go to the documentation of this file.

#include "SimpleRandomNumberGenerator.h"
#include <stdlib.h>
#include <math.h>

/* uniform random number generator */
double SimpleRandomNumberGenerator::uniform_rand(long *idum)


/* (C) Copr. 1986-92 Numerical Recipes Software )'). */
{

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define NTAB 32
#define RNMX (1.0-EPS)


    int j;
    long k;
    double temp;

    if (*idum <= 0) {
        if (-(*idum) < 1) *idum=1;
        else *idum = -(*idum);
        idum2=(*idum);
        for (j=NTAB+7;j>=0;j--) {
            k=(*idum)/IQ1;
            *idum=IA1*(*idum-k*IQ1)-k*IR1;
            if (*idum < 0) *idum += IM1;
            if (j < NTAB) iv[j] = *idum;
        }
        iy=iv[0];
    }
    k=(*idum)/IQ1;
    *idum=IA1*(*idum-k*IQ1)-k*IR1;
    if (*idum < 0) *idum += IM1;
    k=idum2/IQ2;
    idum2=IA2*(idum2-k*IQ2)-k*IR2;
    if (idum2 < 0) idum2 += IM2;
    j=iy/NDIV;
    iy=iv[j]-idum2;
    iv[j] = *idum;
    if (iy < 1) iy += IMM1;
    if ( (temp=AM*iy) > RNMX) return RNMX;
    else return temp;

#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
#undef NTAB

}

void SimpleRandomNumberGenerator::seed(long idum)
{
    IDUM = -labs(idum);
    uniform_rand(&IDUM);
}

double SimpleRandomNumberGenerator::uniformRand(void)
{
    return uniform_rand(&IDUM);
}

double SimpleRandomNumberGenerator::normalRand(void)
{
    double fac,rsq,v1,v2;

    if  (iset == 0) {
        do {
            v1=2.0*uniformRand()-1.0;
            v2=2.0*uniformRand()-1.0;
            rsq=v1*v1+v2*v2;
        } while (rsq >= 1.0 || rsq == 0.0);
        fac=sqrt(-2.0*log(rsq)/rsq);
        gset=v1*fac;
        iset=1;
        return v2*fac;
    } else {
        iset=0;
        return gset;
    }
}

---