From 1ef44553cfef956dff83adb3166190827ec3e308 Mon Sep 17 00:00:00 2001 From: David Schleef Date: Sat, 11 Sep 2010 13:11:14 -0700 Subject: mt19937: Bring up to date with Orc best practice --- examples/Makefile.am | 11 +- examples/mt19937ar.c | 312 ++++++++++++++++++++++++++++++++++------------ examples/mt19937arorc.orc | 85 +++++++++++++ 3 files changed, 327 insertions(+), 81 deletions(-) create mode 100644 examples/mt19937arorc.orc diff --git a/examples/Makefile.am b/examples/Makefile.am index e3ae583..1e32be5 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -6,7 +6,8 @@ orcbin_PROGRAMS += example1 example2 example3 mt19937ar BUILT_SOURCES = example1orc.c example1orc.h \ example2orc.c example2orc.h \ - example3orc.c example3orc.h + example3orc.c example3orc.h \ + mt19937arorc.c mt19937arorc.h endif if ENABLE_BACKEND_MMX @@ -28,6 +29,8 @@ example2_SOURCES = example2.c nodist_example2_SOURCES = example2orc.c example2orc.h example3_SOURCES = example3.c nodist_example3_SOURCES = example3orc.c example3orc.h +mt19937ar_SOURCES = mt19937ar.c +nodist_mt19937ar_SOURCES = mt19937arorc.c mt19937arorc.h example1orc.c: $(srcdir)/example1orc.orc @@ -48,4 +51,10 @@ example3orc.c: $(srcdir)/example3orc.orc example3orc.h: $(srcdir)/example3orc.orc $(top_builddir)/tools/orcc$(EXEEXT) --header -o example3orc.h $(srcdir)/example3orc.orc +mt19937arorc.c: $(srcdir)/mt19937arorc.orc + $(top_builddir)/tools/orcc$(EXEEXT) --implementation -o mt19937arorc.c $(srcdir)/mt19937arorc.orc + +mt19937arorc.h: $(srcdir)/mt19937arorc.orc + $(top_builddir)/tools/orcc$(EXEEXT) --header -o mt19937arorc.h $(srcdir)/mt19937arorc.orc + diff --git a/examples/mt19937ar.c b/examples/mt19937ar.c index 604abbf..a9f736f 100644 --- a/examples/mt19937ar.c +++ b/examples/mt19937ar.c @@ -41,12 +41,11 @@ email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) */ -#ifdef HAVE_CONFIG_H -#include -#endif - #include #include +#include +#include +#include "mt19937arorc.h" /* Period parameters */ #define N 624 @@ -55,108 +54,261 @@ #define UPPER_MASK 0x80000000UL /* most significant w-r bits */ #define LOWER_MASK 0x7fffffffUL /* least significant r bits */ +static unsigned long mt[N]; /* the array for the state vector */ +static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ + +/* initializes mt[N] with a seed */ +void init_genrand(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void init_by_array(orc_uint32 init_key[], int key_length) +{ + int i, j, k; + init_genrand(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long genrand_int32(void) +{ + unsigned long y; + static unsigned long mag01[2]={0x0UL, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if init_genrand() has not been called, */ + init_genrand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 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]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* Orc version */ + +typedef struct _OrcRandomContext OrcRandomContext; +struct _OrcRandomContext { + orc_uint32 d[N]; + orc_uint32 mt[N+1]; + int mti; +}; -/* mag01[x] = x * MATRIX_A for x=0,1 */ -static const orc_uint32 mag01[2]={0x0UL, MATRIX_A}; +OrcRandomContext * +orc_random_context_new (void) +{ + OrcRandomContext *context; + context = malloc(sizeof(OrcRandomContext)); + memset (context, 0, sizeof(OrcRandomContext)); + context->mti = N+1; + return context; +} + +void +orc_random_init_genrand(OrcRandomContext *context, orc_uint32 s) +{ + orc_uint32 *mt = context->mt; + int mti; + + mt[0] = s; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + } + context->mti = mti; +} void -mt19937_ref (orc_uint32 *d, orc_uint32 *mt) +orc_random_init_by_array (OrcRandomContext *context, orc_uint32 *init_key, + int key_length) +{ + int i, j, k; + orc_uint32 *mt = context->mt; + + orc_random_init_genrand (context, 19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +#if 0 +/* These are the functions that were converted to Orc code. */ +static void +mix (orc_uint32 *mt, orc_uint32 *mt2, int n) { orc_uint32 y; int kk; - for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; - } - for (;kk> 1) ^ mag01[y & 0x1UL]; + for (kk=0;kk> 1; + y ^= t1; + + mt[kk] = mt2[kk] ^ y; } - 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> 11); y ^= (y << 7) & 0x9d2c5680UL; y ^= (y << 15) & 0xefc60000UL; y ^= (y >> 18); - - d[kk] = y; + d[i] = y; } } +#endif +static void +update_context (OrcRandomContext *context) +{ + orc_uint32 *d = context->d; + orc_uint32 *mt = context->mt; + +#if 0 + mix (mt, mt + M, N-M); + mt[N] = mt[0]; + mix (mt + N - M, mt, M); + temper (d, mt, N); +#endif -OrcProgram *p1, *p2; + mt19937ar_mix (mt, mt + 1, mt + M, N-M); + mt[N] = mt[0]; + mt19937ar_mix (mt + N - M, mt + N - M + 1, mt, M); + mt19937ar_temper (d, mt, N); -void -create_programs (void) +#if 0 + /* too many temp variables, compiles incorrectly */ + mt19937ar_mix_temper (d, mt, mt + 1, mt + M, N-M); + mt[N] = mt[0]; + mt19937ar_mix_temper (d + N - M, mt + N - M, mt + N - M + 1, mt, M); +#endif +} + +orc_uint32 +orc_random_genrand_int32 (OrcRandomContext *context) { + if (context->mti >= N) { /* generate N words at one time */ + if (context->mti == N+1) /* if init_genrand() has not been called, */ + orc_random_init_genrand(context, 5489UL); /* a default initial seed is used */ - /* 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, 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, "cmpeql", "y", "y", "c3"); - orc_program_append_str (p1, "andl", "y", "y", "c6"); - orc_program_append_str (p1, "xorl", "d1", "t2", "y"); - - orc_program_compile (p1); - - - p2 = orc_program_new (); - orc_program_add_destination (p2, 4, "d1"); - orc_program_add_temporary (p2, 4, "y"); - orc_program_add_temporary (p2, 4, "t1"); - - orc_program_add_constant (p2, 4, 11, "c11"); - orc_program_add_constant (p2, 4, 7, "c7"); - orc_program_add_constant (p2, 4, 0x9d2c5680, "x1"); - orc_program_add_constant (p2, 4, 15, "c15"); - orc_program_add_constant (p2, 4, 0xefc60000, "x2"); - orc_program_add_constant (p2, 4, 15, "c18"); - - orc_program_append_str (p2, "shrul", "t1", "d1", "c11"); - orc_program_append_str (p2, "xorl", "y", "d1", "t1"); - orc_program_append_str (p2, "shll", "t1", "y", "c7"); - orc_program_append_str (p2, "andl", "t1", "t1", "x1"); - orc_program_append_str (p2, "xorl", "y", "y", "t1"); - orc_program_append_str (p2, "shll", "t1", "y", "c15"); - orc_program_append_str (p2, "andl", "t1", "t1", "x2"); - orc_program_append_str (p2, "xorl", "y", "y", "t1"); - orc_program_append_str (p2, "shrul", "t1", "y", "c18"); - orc_program_append_str (p2, "xorl", "d1", "y", "t1"); - - orc_program_compile (p2); + update_context (context); + + context->mti = 0; + } + return context->d[context->mti++]; } -int -main (int argc, char *argv[]) +int main(void) { - orc_init(); - - create_programs (); + int i; + orc_uint32 init[4]={0x123, 0x234, 0x345, 0x456}; + int length=4; + orc_uint32 ref, test; + OrcRandomContext *context; + int error = 0; + + init_by_array(init, length); + + context = orc_random_context_new (); + orc_random_init_by_array (context, init, length); + + printf("1000 outputs of genrand_int32()\n"); + for (i=0; i<1000; i++) { + ref = genrand_int32(); + test = orc_random_genrand_int32(context); + printf("%08x %08x %c\n", ref, test, (ref == test) ? ' ' : '*'); + if (ref != test) error = 1; + } + if (error) { + printf("FAIL\n"); + } return 0; } diff --git a/examples/mt19937arorc.orc b/examples/mt19937arorc.orc new file mode 100644 index 0000000..3c64358 --- /dev/null +++ b/examples/mt19937arorc.orc @@ -0,0 +1,85 @@ + +.function mt19937ar_mix +.dest 4 mt +.source 4 mt1 +.source 4 mt2 +.temp 4 y +.temp 4 t1 +.temp 4 t2 +.const 4 c1 1 +.const 4 UPPER_MASK 0x80000000 +.const 4 LOWER_MASK 0x7fffffff +.const 4 MATRIX_A 0x9908b0df + + +loadl t1, mt +andl t1, t1, UPPER_MASK +loadl t2, mt1 +andl t2, t2, LOWER_MASK +orl y, t1, t2 +andl t1, y, c1 +cmpeql t1, t1, c1 +andl t1, t1, MATRIX_A +shrul y, y, c1 +xorl y, y, t1 +xorl mt, mt2, y + + +.function mt19937ar_temper +.dest 4 d +.source 4 s +.temp 4 y +.temp 4 t + +loadl y, s +shrul t, y, 11 +xorl y, y, t +shll t, y, 7 +andl t, t, 0x9d2c5680 +xorl y, y, t +shll t, y, 15 +andl t, t, 0xefc60000 +xorl y, y, t +shrul t, y, 18 +xorl d, y, t + + +#.function mt19937ar_mix_temper +#.dest 4 d +#.dest 4 mt +#.source 4 mt1 +#.source 4 mt2 +#.temp 4 y +#.temp 4 t1 +#.temp 4 t2 +#.const 4 c1 1 +#.const 4 UPPER_MASK 0x80000000 +#.const 4 LOWER_MASK 0x7fffffff +#.const 4 MATRIX_A 0x9908b0df +# +# +#loadl t1, mt +#andl t1, t1, UPPER_MASK +#loadl t2, mt1 +#andl t2, t2, LOWER_MASK +#orl y, t1, t2 +#andl t1, y, c1 +#cmpeql t1, t1, c1 +#andl t1, t1, MATRIX_A +#shrul y, y, c1 +#xorl y, y, t1 +#xorl y, mt2, y +#storel mt, y +#shrul t1, y, 11 +#xorl y, y, t1 +#shll t1, y, 7 +#andl t1, t1, 0x9d2c5680 +#xorl y, y, t1 +#shll t1, y, 15 +#andl t1, t1, 0xefc60000 +#xorl y, y, t1 +#shrul t1, y, 18 +#xorl d, y, t1 + + + -- cgit v1.2.3