From 69f2911e04ffb1b19eef1fafb8c040af271f656e Mon Sep 17 00:00:00 2001
From: Tor Aamodt <aamodt@ece.ubc.ca>
Date: Thu, 15 Jul 2010 18:09:46 -0800
Subject: creating branch for adding support for CUDA 3.x and Fermi

[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 6829]
---
 src/intersim/rng.cpp | 111 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 111 insertions(+)
 create mode 100644 src/intersim/rng.cpp

(limited to 'src/intersim/rng.cpp')

diff --git a/src/intersim/rng.cpp b/src/intersim/rng.cpp
new file mode 100644
index 0000000..47ddda3
--- /dev/null
+++ b/src/intersim/rng.cpp
@@ -0,0 +1,111 @@
+/*    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<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[100]=-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;
+}
+
-- 
cgit v1.3