randommt.h

#ifndef __randommt_h__
#define __randommt_h__

#define VERSION_A 0
#define VERSION_B 8
#define VERSION_C 459

#include <algorithm>
#include <cassert>
#include <cstdarg>
#include <cstdio>
#include <cstring>
#include <ctime>
#include <limits.h>
#include <stdint.h>

void initrand (uint32_t seed);
void initrand_str (const char * str);
void initrand_array (const unsigned * key, int len);
void initrand_time ();
uint32_t randhash (const char * str);
inline double rndvalue (void);
inline double rndvalue (double range);
inline double rndvalue (double lo, double hi);
inline uint32_t rndvalue (int range);
inline uint32_t rndvalue (uint32_t range);
inline int32_t rndvalue (int32_t lo, int32_t hi);
inline uint64_t rndvalue (uint64_t range);
inline int64_t rndvalue (int64_t lo, int64_t hi);
template <class T> void rndshuffle (T a, T b);

namespace gentools_random_internal
{
 const int POW = 19937, MARR = POW >> 5;
 const int ARR = MARR + 1, EXTRA = (ARR << 5) - POW;
 const int SHIFT = 397, SHIFT_LIM = ARR - SHIFT, MSHIFT = SHIFT - ARR;
 const uint32_t LO_BITS = (1U << EXTRA) - 1U, HI_BITS = 0xFFFFFFFF ^ LO_BITS;
 const uint32_t A [2] = {0U, 0x9908B0DFU};
 const int TSHR1 = 11, TSHL2 = 7, TSHL3 = 15, TSHR4 = 18;
 const uint32_t TMASK2 = 0x9D2C5680U, TMASK3 = 0xEFC60000U;
 const uint32_t RND_MULT1 = 1812433253U;
 const uint32_t RND_MULT2 = 1664525U, RND_MULT3 = 1566083941U;
 const uint32_t DEFAULT_SEED = 5489U, ARRAY_SEED = 19650218U;
 const double NORM_MULT = 1.0 / 4294967296.0;

 const int UTEST_LEN = 4;
 const uint32_t UTEST_KEY [UTEST_LEN] = {0x123, 0x234, 0x345, 0x456};
 const int UTEST_STEPS = 1000;
 const uint32_t UTEST_EXPECTED_RESULT = 0xFA0596C7U;

 uint32_t cur_state [ARR];
 int cur_index = ARR + 1;

 int64_t Time()
 {
  #ifdef __GNUC__
   int64_t res;
   asm volatile ("rdtsc" : "=A" (res));
   return res;
  #else
   int low, hi;
   __asm
   {
     rdtsc
     mov low, eax
     mov hi, edx
   }
   return ((ll)hi << 32) | low;
  #endif
 }

 inline uint32_t generate_element (uint32_t v0, uint32_t v1, uint32_t vshift)
 {
  uint32_t res, temp;

  temp = (v0 & HI_BITS) | (v1 & LO_BITS);
  res = vshift ^ (temp >> 1) ^ A[temp & 1U];

  return res;
 }

 inline void generate_array (void)
 {
  int i;

  for (i = 0; i < SHIFT_LIM; i++)
   cur_state[i] = generate_element (cur_state[i],
                                    cur_state[i + 1],
                                    cur_state[i + SHIFT]);

  for ( ; i < MARR; i++)
   cur_state[i] = generate_element (cur_state[i],
                                    cur_state[i + 1],
                                    cur_state[i + MSHIFT]);

  cur_state[i] = generate_element (cur_state[i],
                                   cur_state[0],
                                   cur_state[i + MSHIFT]);
 }

 inline uint32_t get_raw (void)
 {
  if (cur_index >= ARR)
  {
   if (cur_index > ARR)
    initrand (DEFAULT_SEED);
   assert (cur_index == ARR); // error if not initialized
   cur_index = 0;
   generate_array ();
  }

  return cur_state[cur_index++];
 }

 inline uint32_t get_tempered (void)
 {
  uint32_t val;

  val = get_raw ();
  val ^= val >> TSHR1;
  val ^= (val << TSHL2) & TMASK2;
  val ^= (val << TSHL3) & TMASK3;
  val ^= val >> TSHR4;

  return val;
 }

 inline void unit_test (void)
 {
  uint32_t res;
  int i;

  initrand_array (UTEST_KEY, UTEST_LEN);
  res = 0;
  for (i = 0; i < UTEST_STEPS; i++)
   res = res * RND_MULT1 + rndvalue (0);

  printf ("%08X\n", res);
  assert (res == UTEST_EXPECTED_RESULT);
 }

 struct random_state
 {
  uint32_t cur_state [ARR];
  int cur_index;

  random_state (void): cur_index (ARR + 1) {}
 };

 void save_state (random_state & st)
 {
  memmove (st.cur_state, cur_state, ARR * sizeof (uint32_t));
  st.cur_index = cur_index;
 }

 void restore_state (const random_state & st)
 {
  memmove (cur_state, st.cur_state, ARR * sizeof (uint32_t));
  cur_index = st.cur_index;
 }
} // namespace gentools_random_internal


void initrand (uint32_t seed)
{
 int i;

 using namespace gentools_random_internal;

 assert (seed != 0U);
 cur_state[0] = seed;
 for (i = 1; i < ARR; i++)
  cur_state[i] = RND_MULT1 * (cur_state[i - 1] ^ (cur_state[i - 1] >> 30)) + i;
 cur_index = ARR;
}

void initrand_time (){
	using namespace gentools_random_internal;
	uint32_t t = (uint32_t)Time();
	fprintf (stderr, "Using seed : %u\n", t);
	initrand (t);
}

void initrand_str (const char * str)
{
 int i, j, k, len;

 using namespace gentools_random_internal;

 initrand (ARRAY_SEED);

 len = strlen (str);
 for (i = 1, j = 0, k = std::max (ARR, len); k > 0; k--)
 {
  cur_state[i] = (cur_state[i] ^ (RND_MULT2 *
                                  (cur_state[i - 1] ^
                                   (cur_state[i - 1] >> 30)))) +
                 (uint32_t) str[j] + j;
  i++;
  if (i >= ARR)
  {
   cur_state[0] = cur_state[MARR];
   i = 1;
  }
  j++;
  if (j >= len)
   j = 0;
 }

 for (k = MARR; k > 0; k--)
 {
  cur_state[i] = (cur_state[i] ^ (RND_MULT3 *
                                  (cur_state[i - 1] ^
                                   (cur_state[i - 1] >> 30)))) - i;
  i++;
  if (i >= ARR)
  {
   cur_state[0] = cur_state[MARR];
   i = 1;
  }
 }

 cur_state[0] = 1U << 31; // only the highest bit is relevant here
 cur_index = ARR;
}

// strip all characters with ASCII codes 32 (Space) or less
void initrand_str_strip (const char * str)
{
 int old_len = strlen (str);
 char new_str [old_len + 1];
 char * new_p;
 const char * old_p;

 new_p = new_str;
 for (old_p = str; *old_p != '\0'; old_p++)
  if (*old_p > 32)
   *(new_p++) = *old_p;
 *new_p = '\0';

 initrand_str (new_str);
}

void initrand_array (const unsigned * key, int len)
{
 int i, j, k;

 using namespace gentools_random_internal;

 initrand (ARRAY_SEED);

 for (i = 1, j = 0, k = std::max (ARR, len); k > 0; k--)
 {
  cur_state[i] = (cur_state[i] ^ (RND_MULT2 *
                                  (cur_state[i - 1] ^
                                   (cur_state[i - 1] >> 30)))) +
                 key[j] + j;
  i++;
  if (i >= ARR)
  {
   cur_state[0] = cur_state[MARR];
   i = 1;
  }
  j++;
  if (j >= len)
   j = 0;
 }

 for (k = MARR; k > 0; k--)
 {
  cur_state[i] = (cur_state[i] ^ (RND_MULT3 *
                                  (cur_state[i - 1] ^
                                   (cur_state[i - 1] >> 30)))) - i;
  i++;
  if (i >= ARR)
  {
   cur_state[0] = cur_state[MARR];
   i = 1;
  }
 }

 cur_state[0] = 1U << 31; // only the highest bit is relevant here
 cur_index = ARR;
}

uint32_t randhash (const char * str)
{
 int i;
 uint32_t res;

 using namespace gentools_random_internal;

 res = 0;
 for (i = 0; str[i] != '\0'; i++)
  res = RND_MULT1 * (res ^ (res >> 30)) + (uint32_t) str[i];

 return res;
}

inline double rndvalue (void)
{
 using namespace gentools_random_internal;

 return get_tempered () * NORM_MULT;
}

inline double rndvalue (double range)
{
 return rndvalue () * (range);
}

inline double rndvalue (double lo, double hi)
{
 return rndvalue (hi - lo) + lo;
}

inline uint32_t rndvalue (uint32_t range)
{
 using namespace gentools_random_internal;

 if (!range)
  return get_tempered ();
 else
  return (uint32_t) (((uint64_t) get_tempered () *
                      range + (range >> 1)) >> 32ULL);
}

inline uint32_t rndvalue (int range)
{
 return rndvalue ((uint32_t) range);
}

inline int32_t rndvalue (int32_t lo, int32_t hi)
{
 return lo + (int32_t) rndvalue (hi - lo + 1);
}

inline uint64_t rndvalue (uint64_t range)
{
 using namespace gentools_random_internal;

 uint64_t res;
 if (!range)
 {
  res = (uint64_t) get_tempered () << 32ULL;
  res |= (uint64_t) get_tempered ();
 }
 else
 {
  uint64_t lim;
  lim = uint64_t (-1);
  if (range <= (lim >> 32ULL))
   lim >>= 32ULL;
  if (range <= (lim >> 16ULL))
   lim >>= 16ULL;
  if (range <= (lim >>  8ULL))
   lim >>=  8ULL;
  if (range <= (lim >>  4ULL))
   lim >>=  4ULL;
  if (range <= (lim >>  2ULL))
   lim >>=  2ULL;
  if (range <= (lim >>  1ULL))
   lim >>=  1ULL;
  assert (range <= lim);
  do
  {
   res = (uint64_t) get_tempered () << 32ULL;
   res |= (uint64_t) get_tempered ();
   res &= lim;
  }
  while (res >= range);
 }

 return res;
}

inline int64_t rndvalue (int64_t lo, int64_t hi)
{
 return lo + (int64_t) rndvalue ((uint64_t) (hi - lo + 1));
}

template <class T> void rndshuffle (T a, T b)
{
 int i, k;

 for (i = 0; a + i != b; i++)
 {
  k = rndvalue (i + 1);
  std::swap (a[i], a[k]);
 }
}


#endif // __randommt_h__