A non-trivial Mersenne Twister algorithm

2 files changed, 142 insertions, 1 deletions

diff --git a/examples/Makefile.am b/examples/Makefile.am
index c91a3d8..f9783be 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -1,5 +1,5 @@
 
-noinst_PROGRAMS = jit simple
+noinst_PROGRAMS = jit simple mt19937ar
 
 AM_LDFLAGS = $(ORC_LIBS)
 AM_CFLAGS = $(ORC_CFLAGS)
diff --git a/examples/mt19937ar.c b/examples/mt19937ar.c
new file mode 100644
index 0000000..d939913
--- /dev/null
+++ b/examples/mt19937ar.c
@@ -0,0 +1,141 @@
+/* 
+   A C-program for MT19937, with initialization improved 2002/1/26.
+   Coded by Takuji Nishimura and Makoto Matsumoto.
+
+   Before using, initialize the state by using init_genrand(seed)  
+   or init_by_array(init_key, key_length).
+
+   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+   All rights reserved.                          
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+     1. Redistributions of source code must retain the above copyright
+        notice, this list of conditions and the following disclaimer.
+
+     2. Redistributions in binary form must reproduce the above copyright
+        notice, this list of conditions and the following disclaimer in the
+        documentation and/or other materials provided with the distribution.
+
+     3. The names of its contributors may not be used to endorse or promote 
+        products derived from this software without specific prior written 
+        permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+   Any feedback is very welcome.
+   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <orc/orc.h>
+
+/* Period parameters */  
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
+#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
+#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
+
+
+/* mag01[x] = x * MATRIX_A  for x=0,1 */
+static const uint32_t mag01[2]={0x0UL, MATRIX_A};
+
+void
+mt19937_ref (uint32_t *d, uint32_t *mt)
+{
+  uint32_t y;
+  int kk;
+
+  for (kk=0;kk<N-M;kk++) {
+      y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+      mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+  }
+  for (;kk<N-1;kk++) {
+      y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+      mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+  }
+  y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+  mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+  for(kk=0;kk<N;kk++){
+    y = mt[kk];
+
+    /* Tempering */
+    y ^= (y >> 11);
+    y ^= (y << 7) & 0x9d2c5680UL;
+    y ^= (y << 15) & 0xefc60000UL;
+    y ^= (y >> 18);
+
+    d[kk] = y;
+  }
+}
+
+
+OrcProgram *p1, *p2, *p3;
+
+void
+create_programs (void)
+{
+
+  /* s1 is mt, s2 is mt+1, s3 is mt+M */
+  p1 = orc_program_new_dss (4, 4, 4);
+  orc_program_add_source (p1, 4, "s3");
+
+  orc_program_add_temporary (p1, 4, "t1");
+  orc_program_add_temporary (p1, 4, "t2");
+  orc_program_add_temporary (p1, 4, "y");
+
+  orc_program_add_constant (p1, 4, UPPER_MASK, "c1");
+  orc_program_add_constant (p1, 4, LOWER_MASK, "c2");
+  orc_program_add_constant (p1, 4, 1, "c3");
+  orc_program_add_constant (p1, 4, 31, "c5");
+  orc_program_add_constant (p1, 4, MATRIX_A, "c6");
+
+  orc_program_append_str (p1, "andl", "t1", "s1", "c1");
+  orc_program_append_str (p1, "andl", "t2", "s2", "c2");
+  orc_program_append_str (p1, "orl", "y", "t1", "t2");
+
+  orc_program_append_str (p1, "shrul", "t1", "y", "c3");
+  orc_program_append_str (p1, "xorl", "t2", "s3", "t1");
+  orc_program_append_str (p1, "andl", "y", "y", "c3");
+  orc_program_append_str (p1, "shll", "y", "y", "c5");
+  orc_program_append_str (p1, "shrsl", "y", "y", "c5");
+  orc_program_append_str (p1, "andl", "y", "y", "c6");
+  orc_program_append_str (p1, "xorl", "d1", "t2", "y");
+  
+  orc_program_compile (p1);
+
+  printf("%s", orc_program_get_asm_code (p1));
+
+}
+
+
+int
+main (int argc, char *argv[])
+{
+  orc_init();
+
+  create_programs ();
+
+  return 0;
+}
+