diff options
| author | Tor Aamodt <[email protected]> | 2010-07-15 18:09:46 -0800 |
|---|---|---|
| committer | Tor Aamodt <[email protected]> | 2010-07-15 18:09:46 -0800 |
| commit | 69f2911e04ffb1b19eef1fafb8c040af271f656e (patch) | |
| tree | 231d3b6bdc3a202f7c255bfcf7bf2c36e32cee9e /src/intersim/rng_double.cpp | |
creating branch for adding support for CUDA 3.x and Fermi
[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 6829]
Diffstat (limited to 'src/intersim/rng_double.cpp')
| -rw-r--r-- | src/intersim/rng_double.cpp | 115 |
1 files changed, 115 insertions, 0 deletions
diff --git a/src/intersim/rng_double.cpp b/src/intersim/rng_double.cpp new file mode 100644 index 0000000..58a3c2c --- /dev/null +++ b/src/intersim/rng_double.cpp @@ -0,0 +1,115 @@ +/* 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
+ * (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<KK;j++) aa[j]=ran_u[j];
+ for (;j<n;j++) aa[j]=mod_sum(aa[j-KK],aa[j-LL]);
+ for (i=0;i<LL;i++,j++) ran_u[i]=mod_sum(aa[j-KK],aa[j-LL]);
+ for (;i<KK;i++,j++) ran_u[i]=mod_sum(aa[j-KK],ran_u[i-LL]);
+}
+
+/* the following routines are adapted from exercise 3.6--15 */
+/* after calling ranf_start, get new randoms by, e.g., "x=ranf_arr_next()" */
+
+#define QUALITY 1009 /* recommended quality level for high-res use */
+double ranf_arr_buf[QUALITY];
+double ranf_arr_dummy=-1.0, ranf_arr_started=-1.0;
+double *ranf_arr_ptr=&ranf_arr_dummy; /* the next random fraction, or -1 */
+
+#define TT 70 /* guaranteed separation between streams */
+#define is_odd(s) ((s)&1)
+
+#ifdef __STDC__
+void ranf_start(long seed)
+#else
+void ranf_start(seed) /* do this before using ranf_array */
+long seed; /* selector for different streams */
+#endif
+{
+ register int t,s,j;
+ double u[KK+KK-1];
+ double ulp=(1.0/(1L<<30))/(1L<<22); /* 2 to the -52 */
+ double ss=2.0*ulp*((seed&0x3fffffff)+2);
+
+ for (j=0;j<KK;j++) {
+ u[j]=ss; /* bootstrap the buffer */
+ ss+=ss; if (ss>=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<LL;j++) ran_u[j+KK-LL]=u[j];
+ for (;j<KK;j++) ran_u[j-LL]=u[j];
+ for (j=0;j<10;j++) ranf_array(u,KK+KK-1); /* warm things up */
+ ranf_arr_ptr=&ranf_arr_started;
+}
+
+#define ranf_arr_next() (*ranf_arr_ptr>=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[100]=-1;
+ ranf_arr_ptr=ranf_arr_buf+1;
+ return ranf_arr_buf[0];
+}
+
+#include <stdio.h>
+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;
+}
+
|