/* This program by D E Knuth is in the public domain and freely copyable
* AS LONG AS YOU MAKE ABSOLUTELY NO CHANGES!
* It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
* (or in the errata to the 2nd edition, see
* https://www-cs-faculty.stanford.edu/~knuth/taocp.html
* in the changes to pages 171 and following).
*
* If you find any bugs, please report them immediately to
* taocp@cs.stanford.edu
* (and you will be rewarded if the bug is genuine). Thanks! */
#define KK 100 /* the long lag */
#define LL 37 /* the short lag */
#define MM (1<<30) /* the modulus */
#define mod_diff(x,y) (x-y)&(MM-1) /* subtraction mod MM */
long ran_x[KK]; /* the generator state */
void ran_array(aa,n) /* put n new random numbers in aa */
long *aa; /* destination */
int n; /* array length (must be at least KK) */
{
register int i,j;
for (j=0;j<KK;j++) aa[j]=ran_x[j];
for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
}
#define TT 70 /* guaranteed separation between streams */
#define is_odd(x) (x&1) /* units bit of x */
#define evenize(x) ((x)&(MM-2)) /* make x even */
void ran_start(seed) /* do this before using ran_array */
long seed; /* selector for different streams */
{
register int t,j;
long x[KK+KK-1]; /* the preparation buffer */
register long ss=evenize(seed+2);
for (j=0;j<KK;j++) {
x[j]=ss; /* bootstrap the buffer */
ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
}
for (;j<KK+KK-1;j++) x[j]=0;
x[1]++; /* make x[1] (and only x[1]) odd */
ss=seed;
t=TT-1; while (t) {
for (j=KK-1;j>0;j--) x[j+j]=x[j]; /* "square" */
for (j=KK+KK-2;j>KK-LL;j-=2) x[KK+KK-1-j]=evenize(x[j]);
for (j=KK+KK-2;j>=KK;j--) if(is_odd(x[j])) {
x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]);
x[j-KK]=mod_diff(x[j-KK],x[j]);
}
if (is_odd(ss)) { /* "multiply by z" */
for (j=KK;j>0;j--) x[j]=x[j-1];
x[0]=x[KK]; /* shift the buffer cyclically */
if (is_odd(x[KK])) x[LL]=mod_diff(x[LL],x[KK]);
}
if (ss) ss>>=1; else t--;
}
for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
for (;j<KK;j++) ran_x[j-LL]=x[j];
}