lRandom Number Generators（随机数生成）

Random Number Generators
 
gran.c - gaussian random number generator
 
ran.c - uniform random number generator
 
ran1f.c - 1/f noise generator
 
ranh.c - low-frequency hold generator
 
ranl.c - linearly interpolated generator

----------------------------------------------------------------------------------------------------------------

/* ranl.c - linearly interpolated random generator of period D */

double ran();                             /* uniform generator */

void cdelay2();                           /* circular delay */

double ranl(D, u, q, iseed)               /* usage: y = ranl(D, u, &q, &iseed); */

int D, *q;                                /* \(q\) is cycled modulo-\(D\) */

double *u;                                /* \(u\) = 2-dimensional array */

long *iseed;                              /* \(q\), iseed are passed by address */

{

       double y;

       int i;

       i = (D - *q) % D;                      /* interpolation index */

       y = u[0] + (u[1] - u[0]) * i / D;      /* linear interpolation */

       cdelay2(D-1, q);                       /* decrement \(q\) and wrap mod-\(D\) */

       if (*q == 0) {                         /* every \(D\) calls, */

              u[0] = u[1];                    /* set new \(u[0]\) and */

              u[1] = ran(iseed) - 0.5;        /* get new \(u[1]\) (zero mean) */

              }

       return y;

}

----------------------------------------------------------------------------------------------------------------

/* ranh.c - hold random number generator of period D */

double ran();                             /* uniform generator */

void cdelay2();                           /* circular delay */

double ranh(D, u, q, iseed)               /* usage: y = ranh(D, u, &q, &iseed); */

int D, *q;                                /* \(q\) is cycled modulo-\(D\) */

double *u;                                /* \(u\) = 1-dimensional array */

long *iseed;                              /* \(q\), iseed are passed by address */

{

       double y;

       y = u[0];                              /* hold sample for \(D\) calls */

       cdelay2(D-1, q);                       /* decrement \(q\) and wrap mod-\(D\) */

       if (*q == 0)                           /* every \(D\) calls, */

              u[0] = ran(iseed) - 0.5;        /* get new \(u[0]\) (zero mean) */

       return y;

}

----------------------------------------------------------------------------------------------------------------

/* ran1f.c - 1/f random number generator */

double ranh();                              /* random hold periodic generator */

double ran1f(B, u, q, iseed)                /* usage: y = ran1f(B, u, q, &iseed); */

int B, *q;                                  /* \(q, u\) are \(B\)-dimensional */

double *u;

long *iseed;                                /* passed by address */

{

    double y;

    int b;

    for(y=0, b=0; b<B; b++)

          y += ranh(1<<b, u+b, q+b, iseed);         /* period = (1<<b) = 2\(\sp{b}\) */

    return y / B;

}

----------------------------------------------------------------------------------------------------------------

/* ran.c - uniform random number generator in [0, 1) */

#define  a    16807                              /* that is, \(a = 7\sp{5}\) */

#define  m    2147483647                         /* that is, \(m = 2\sp{31}-1\) */

#define  q    127773                             /* note, \(q = m/a\) = quotient */

#define  r    2836                               /* note, \(r = m\%a\) = remainder */

double ran(iseed)                                /* usage: u = ran(&iseed); */

long *iseed;                                     /* iseed passed by address */

{

    *iseed = a * (*iseed % q) - r * (*iseed / q);          /* update seed */

    if (*iseed < 0)                              /* wrap to positive values */

           *iseed += m;

    return (double) *iseed / (double) m;

}

----------------------------------------------------------------------------------------------------------------

/* gran.c - gaussian random number generator */

double ran();                            /* uniform generator */

double gran(m, s, iseed)                 /* usage: x = gran(m, s, &iseed); */

double m, s;                             /* \(m\) = mean, \(s\sp{2}\) = variance */

long *iseed;                             /* iseed passed by address */

{

       double v = 0;

       int i;

       for (i = 0; i < 12; i++)          /* sum 12 uniform random numbers */

              v += ran(iseed);

       return s * (v - 6) + m;           /* adjust mean and variance */

}