/* 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 * http://www-cs-faculty.stanford.edu/~knuth/taocp.html * in the changes to Volume 2 on pages 171 and following). */ /* N.B. The MODIFICATIONS introduced in the 9th printing (2002) are included here; there's no backwards compatibility with the original. */ /* This version also adopts Brendan McKay's suggestion to accommodate naive users who forget to call ran_start(seed). */ /* 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! */ /************ see the book for explanations and caveats! *******************/ /************ in particular, you need two's complement arithmetic **********/ #define KK 100 /* the long lag */ #define LL 37 /* the short lag */ #define MM (1L<<30) /* the modulus */ #define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */ long ran_x[KK]; /* the generator state */ #ifdef __STDC__ void ran_array(long aa[],int n) #else void ran_array(aa,n) /* put n new random numbers in aa */ long *aa; /* destination */ int n; /* array length (must be at least KK) */ #endif { register int i,j; for (j=0;j=MM) ss-=MM-2; /* cyclic shift 29 bits */ } x[1]++; /* make x[1] (and only x[1]) odd */ for (ss=seed&(MM-1),t=TT-1; t; ) { for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */ for (j=KK+KK-2;j>=KK;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 */ x[LL]=mod_diff(x[LL],x[KK]); } if (ss) ss>>=1; else t--; } for (j=0;j=0? *ran_arr_ptr++: ran_arr_cycle()) long ran_arr_cycle() { if (ran_arr_ptr==&ran_arr_dummy) ran_start(314159L); /* the user forgot to initialize */ ran_array(ran_arr_buf,QUALITY); ran_arr_buf[KK]=-1; ran_arr_ptr=ran_arr_buf+1; return ran_arr_buf[0]; } #include int main() { register int m; long a[2009]; ran_start(310952L); for (m=0;m<=2009;m++) ran_array(a,1009); printf("%ld\n", a[0]); /* 995235265 */ ran_start(310952L); for (m=0;m<=1009;m++) ran_array(a,2009); printf("%ld\n", a[0]); /* 995235265 */ return 0; }