diff options
| author | Dongdong Li <[email protected]> | 2013-08-08 00:15:58 -0800 |
|---|---|---|
| committer | Andrew Boktor <[email protected]> | 2014-08-14 13:50:58 -0700 |
| commit | 7f49fe9feb174d34efc2a011bad79b38522a360b (patch) | |
| tree | ab5b7b66e40315a81871acbf386722981020f866 /src/intersim2/rng.c | |
| parent | 5f91e7435742bab74dfbeca18afc63e466498f36 (diff) | |
Intesim2 Integration
Details: See Review 80001 https://gpgpu-sim-code-review.appspot.com/80001/
[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 16747]
Diffstat (limited to 'src/intersim2/rng.c')
| -rw-r--r-- | src/intersim2/rng.c | 108 |
1 files changed, 108 insertions, 0 deletions
diff --git a/src/intersim2/rng.c b/src/intersim2/rng.c new file mode 100644 index 0000000..c27740f --- /dev/null +++ b/src/intersim2/rng.c @@ -0,0 +1,108 @@ +/* This program by D E Knuth is in the public domain and freely copyable. + * 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 + * (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<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]); +} + +/* the following routines are from exercise 3.6--15 */ +/* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */ + +#define QUALITY 1009 /* recommended quality level for high-res use */ +long ran_arr_buf[QUALITY]; +long ran_arr_dummy=-1, ran_arr_started=-1; +long *ran_arr_ptr=&ran_arr_dummy; /* the next random number, or -1 */ + +#define TT 70 /* guaranteed separation between streams */ +#define is_odd(x) ((x)&1) /* units bit of x */ + +#ifdef __STDC__ +void ran_start(long seed) +#else +void ran_start(seed) /* do this before using ran_array */ + long seed; /* selector for different streams */ +#endif +{ + register int t,j; + long x[KK+KK-1]; /* the preparation buffer */ + register long ss=(seed+2)&(MM-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 */ + } + 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<LL;j++) ran_x[j+KK-LL]=x[j]; + for (;j<KK;j++) ran_x[j-LL]=x[j]; + for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */ + ran_arr_ptr=&ran_arr_started; +} + +#define ran_arr_next() (*ran_arr_ptr>=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 <stdio.h> +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; +} |
