/* 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 */
}



/* 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;
}


/* 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;
}


/* 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;
}


/* 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;
}