From 61b30811c138dd3b011b0d0e5c7c265c05fa0893 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Daniel=20Kochma=C5=84ski?= Date: Tue, 22 Sep 2015 22:37:48 +0200 Subject: [PATCH] random: provide 64 bit implementation MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit This change pulls dependency on C99 types Signed-off-by: Daniel Kochmański --- src/c/num_rand.d | 161 +++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 120 insertions(+), 41 deletions(-) diff --git a/src/c/num_rand.d b/src/c/num_rand.d index ab13ea1..2ec7ffe 100644 --- a/src/c/num_rand.d +++ b/src/c/num_rand.d @@ -36,64 +36,120 @@ n*/ * Mersenne-Twister random number generator */ -/* Period parameters */ -#define MT_N 624 -#define MT_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 */ +#if ECL_FIXNUM_BITS > 32 +/* + * 64 bit version + */ -/* INV: for 64 bit implementation modify accordingly _hash_equal in - hash.d file for case t_random */ -#define ulong unsigned long +/* Period parameters */ +#define MT_N 312 +#define MT_M 156 +#define MATRIX_A 0xB5026F5AA96619E9ULL /* constant vector a */ +#define UPPER_MASK 0xFFFFFFFF80000000ULL /* most significant 33 bits */ +#define LOWER_MASK 0x7FFFFFFFULL /* least significant 31 bits */ +#define ulong uint64_t -cl_object +static cl_object init_genrand(ulong seed) { - cl_object array = - ecl_alloc_simple_vector - ((MT_N + 1) * sizeof(ulong), - ecl_symbol_to_elttype(@'ext::byte8')); - ulong *mt = (ulong*)(array->vector.self.b8); - int j = 0; - mt[0] = seed & 0xffffffffUL; - for (j=1; j < MT_N; j++) - mt[j] = (1812433253UL * (mt[j-1] ^ (mt[j-1] >> 30)) + j) - & 0xffffffffUL; + cl_object array = ecl_alloc_simple_vector((MT_N + 1), ecl_aet_b64); + ulong *mt = array->vector.self.b64; + mt[0] = seed; + for (int j=1; j> 62)) + j); mt[MT_N] = MT_N+1; return array; } -cl_object -init_random_state(void) +static ulong +generate_int64(cl_object state) { - ulong seed; -#if !defined(ECL_MS_WINDOWS_HOST) - /* fopen() might read full 4kB blocks and discard - * a lot of entropy, so use open() */ - int file_handler = open("/dev/urandom", O_RDONLY); - if (file_handler != -1) { - read(file_handler, &seed, sizeof(ulong)); - close(file_handler); - } else -#endif - { - /* cant get urandom, use crappy source */ - /* and/or fill rest of area */ - seed = (rand() + time(0)); + static mag01[2]={0x0UL, MATRIX_A}; + ulong y; + ulong *mt = state->vector.self.b64; + + if (mt[MT_N] >= MT_N) { + /* refresh data */ + int kk; + for (kk=0; kk < (MT_N - MT_M); kk++) { + y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); + mt[kk] = mt[kk + MT_M] ^ (y >> 1) ^ mag01[y & 0x1ULL]; + } + for (; kk < (MT_N - 1); kk++) { + y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); + mt[kk] = mt[kk+(MT_M-MT_N)] ^ (y >> 1) ^ mag01[y & 0x1ULL]; + } + y = (mt[MT_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK); + mt[MT_N-1] = mt[MT_M-1] ^ (y >> 1) ^ mag01[y & 0x1ULL]; + mt[MT_N] = 0; } + /* get random 64 bit num */ + y = mt[mt[MT_N]++]; + /* Tempering */ + y ^= (y >> 29) & 0x5555555555555555ULL; + y ^= (y << 17) & 0x71D67FFFEDA60000ULL; + y ^= (y << 37) & 0xFFF7EEE000000000ULL; + y ^= (y >> 43); - return init_genrand(seed); + return y; +} + +static double +generate_double(cl_object state) +{ + return (generate_int64(state) >> 11) * (1.0 / 9007199254740991.0); +} + +static mp_limb_t +generate_limb(cl_object state) +{ +#if GMP_LIMB_BITS <= 32 + return generate_int64(state); +#else +# if GMP_LIMB_BITS <= 64 + return generate_int64(state); +# else +# if GMP_LIMB_BITS <= 128 + mp_limb_t high = generate_int64(state); + return (high << 64) | generate_int64(state); +# endif +# endif +#endif } +#else +/* + * 32 bit version + */ -ulong +/* Period parameters */ +#define MT_N 624 +#define MT_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 */ +#define ulong uint32_t + +static cl_object +init_genrand(ulong seed) +{ + cl_object array = ecl_alloc_simple_vector((MT_N + 1), ecl_aet_b32); + ulong *mt = array->vector.self.b32; + mt[0] = seed; + for (int j=1; j < MT_N; j++) + mt[j] = (1812433253UL * (mt[j-1] ^ (mt[j-1] >> 30)) + j); + + mt[MT_N] = MT_N+1; + return array; +} + +static ulong generate_int32(cl_object state) { - static ulong mag01[2]={0x0UL, MATRIX_A}; + static mag01[2]={0x0UL, MATRIX_A}; ulong y; - ulong *mt = (ulong*)state->vector.self.b8; - if (mt[MT_N] >= MT_N){ + ulong *mt = state->vector.self.b32; + if (mt[MT_N] >= MT_N) { /* refresh data */ int kk; for (kk=0; kk < (MT_N - MT_M); kk++) { @@ -144,6 +200,29 @@ generate_limb(cl_object state) # endif #endif } +#endif + +cl_object +init_random_state(void) +{ + ulong seed; +#if !defined(ECL_MS_WINDOWS_HOST) + /* fopen() might read full 4kB blocks and discard + * a lot of entropy, so use open() */ + int file_handler = open("/dev/urandom", O_RDONLY); + if (file_handler != -1) { + read(file_handler, &seed, sizeof(ulong)); + close(file_handler); + } else +#endif + { + /* cant get urandom, use crappy source */ + /* and/or fill rest of area */ + seed = (rand() + time(0)); + } + + return init_genrand(seed); +} static cl_object random_integer(cl_object limit, cl_object state) -- 2.9.0