diff options
| -rw-r--r-- | rngtest.c | 81 | ||||
| -rw-r--r-- | stc/crandom.h | 78 |
2 files changed, 21 insertions, 138 deletions
@@ -7,73 +7,6 @@ #endif - -/* 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 */ - -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<N; mti++) { - mt[mti] = - (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 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 */ - } -} - -/* 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<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]; - - mti = 0; - } - - y = mt[mti++]; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - - - #define NN 1000000000 int main(void) @@ -82,15 +15,6 @@ int main(void) uint64_t v; printf("start\n"); - before = clock(); \ - v = 0; - for (size_t i=0; i<NN; i++) { - v += genrand_int32(); - } - difference = clock() - before; - printf("refmt: %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v); - - mt19937_t state = mt19937_default(); before = clock(); \ @@ -99,7 +23,7 @@ int main(void) v += mt19937_rand(&state); } difference = clock() - before; - printf("mymt : %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v); + printf("my-mt: %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v); #ifdef __cplusplus std::mt19937 mt_rand; @@ -131,7 +55,7 @@ int main(void) difference = clock() - before; printf("sfc64: %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v); - +/* before = clock(); \ v = 0; for (size_t i=0; i<NN; i++) { @@ -139,4 +63,5 @@ int main(void) } difference = clock() - before; printf("rand : %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v); +*/ }
\ No newline at end of file diff --git a/stc/crandom.h b/stc/crandom.h index 064c86eb..2bf88028 100644 --- a/stc/crandom.h +++ b/stc/crandom.h @@ -1,75 +1,32 @@ #ifndef CRANDOM__H__ #define CRANDOM__H__ -/* - Mersenne Twister random number generator MT19937, 32 bit. - A C-program for MT19937, with initialization improved 2002/1/26. - Coded by Takuji Nishimura and Makoto Matsumoto. - Optimized by Tyge Løvset. - - Before using, initialize the state by using mt19937_init(), - mt19937_seed() or mt19937_default(). - - 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) -*/ #include <stdint.h> -enum { /* period parameters */ +/* + * Mersenne Twister random number generator MT19937, 32 bit. + */ + +enum { mt19937_N = 624, mt19937_M = 397, }; typedef struct mt19937 { - uint_fast32_t idx; - uint_fast32_t arr[mt19937_N]; /* the array for the state vector */ + uint32_t idx; + uint32_t arr[mt19937_N]; } mt19937_t; /* initializes state with a seed */ -static inline void mt19937_init(mt19937_t *state, uint_fast32_t seed) { +static inline void mt19937_init(mt19937_t *state, uint32_t seed) { state->idx = mt19937_N; state->arr[0] = seed; for (int i = 1; i < mt19937_N; ++i) { seed = state->arr[i] = 1812433253 * (seed ^ (seed >> 30)) + i; - /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ - /* In the previous versions, MSBs of the seed affect */ - /* only MSBs of the array arr[]. */ - /* 2002/01/09 modified by Makoto Matsumoto */ } } /* creates a new state from a seed */ -static inline mt19937_t mt19937_seed(uint_fast32_t seed) { +static inline mt19937_t mt19937_seed(uint32_t seed) { mt19937_t state; mt19937_init(&state, seed); return state; @@ -82,10 +39,10 @@ static inline mt19937_t mt19937_default(void) { } /* generates a random number on [0, 0xffffffff]-interval */ -static inline uint_fast32_t mt19937_rand(mt19937_t *state) { +static inline uint32_t mt19937_rand(mt19937_t *state) { enum {N = mt19937_N, M = mt19937_M}; - uint_fast32_t y, *arr = state->arr; + uint32_t y, *arr = state->arr; if (state->idx >= N) { /* generate N words at one time */ int k = 0; for (; k < N-M; ++k) { @@ -110,7 +67,9 @@ static inline uint_fast32_t mt19937_rand(mt19937_t *state) { return y; } -/* xoroshiro128** and xoshiro128** with splitmix64 seed initialization */ +/* + * xoroshiro128** with splitmix64 seed initialization + */ static inline uint64_t splitmix64(uint64_t *state) { uint64_t z = (*state += 0x9e3779b97f4a7c15); @@ -123,8 +82,6 @@ static inline uint64_t c_rotl(uint64_t x, int s) { return (x << s) | (x >> (64 - s)); } -/* xoroshiro128** */ - typedef struct xoroshiro128ss { uint64_t s[2]; } xoroshiro128ss_t; @@ -146,13 +103,14 @@ static inline uint64_t xoroshiro128ss_rand(xoroshiro128ss_t *state) { return result; } -/* http://pracrand.sourceforge.net */ +/* + * sfc64: http://pracrand.sourceforge.net + */ typedef struct sfc64 { uint64_t s[3], counter; } sfc64_t; - static inline uint64_t sfc64_rand(sfc64_t* state) { enum {LROT = 24, RSHIFT = 11, LSHIFT = 3}; uint64_t *s = state->s; @@ -160,7 +118,7 @@ static inline uint64_t sfc64_rand(sfc64_t* state) { uint64_t result = s[0] + s[1] + ++state->counter; s[0] = s[1] ^ (s[1] >> RSHIFT); s[1] = s[2] + (s[2] << LSHIFT); - s[2] = ((s[2] << LROT) | (s[2] >> (64-LROT))) + result; + s[2] = c_rotl(s[2], LROT) + result; return result; } |