/* 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 ranf_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 mod_sum(x,y) (((x)+(y))-(int)((x)+(y))) /* (x+y) mod 1.0 */ double ran_u[KK]; /* the generator state */ #ifdef __STDC__ void ranf_array(double aa[], int n) #else void ranf_array(aa,n) /* put n new random fractions in aa */ double *aa; /* destination */ int n; /* array length (must be at least KK) */ #endif { register int i,j; for (j=0;j=1.0) ss-=1.0-2*ulp; /* cyclic shift of 51 bits */ } u[1]+=ulp; /* make u[1] (and only u[1]) "odd" */ for (s=seed&0x3fffffff,t=TT-1; t; ) { for (j=KK-1;j>0;j--) u[j+j]=u[j],u[j+j-1]=0.0; /* "square" */ for (j=KK+KK-2;j>=KK;j--) { u[j-(KK-LL)]=mod_sum(u[j-(KK-LL)],u[j]); u[j-KK]=mod_sum(u[j-KK],u[j]); } if (is_odd(s)) { /* "multiply by z" */ for (j=KK;j>0;j--) u[j]=u[j-1]; u[0]=u[KK]; /* shift the buffer cyclically */ u[LL]=mod_sum(u[LL],u[KK]); } if (s) s>>=1; else t--; } for (j=0;j=0? *ranf_arr_ptr++: ranf_arr_cycle()) double ranf_arr_cycle() { if (ranf_arr_ptr==&ranf_arr_dummy) ranf_start(314159L); /* the user forgot to initialize */ ranf_array(ranf_arr_buf,QUALITY); ranf_arr_buf[KK]=-1; ranf_arr_ptr=ranf_arr_buf+1; return ranf_arr_buf[0]; } #include int main() { register int m; double a[2009]; /* a rudimentary test */ ranf_start(310952); for (m=0;m<2009;m++) ranf_array(a,1009); printf("%.20f\n", ran_u[0]); /* 0.36410514377569680455 */ /* beware of buggy printf routines that do not give full accuracy here! */ ranf_start(310952); for (m=0;m<1009;m++) ranf_array(a,2009); printf("%.20f\n", ran_u[0]); /* 0.36410514377569680455 */ return 0; }