103 lines
2.4 KiB
C
103 lines
2.4 KiB
C
/**
|
|
* twister.c
|
|
*
|
|
* Creation Date: 09/05/2019
|
|
*
|
|
* Authors:
|
|
* Leonardo Vencovsky (https://github.com/LeoVen)
|
|
*
|
|
* Mersenne Twister PRNG utility
|
|
*
|
|
*/
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <stdint.h>
|
|
|
|
#define NN 312
|
|
#define MM 156
|
|
#define MATRIX_A UINT64_C(0xB5026F5AA96619E9)
|
|
#define UM UINT64_C(0xFFFFFFFF80000000) // Most significant 33 bits
|
|
#define LM UINT64_C(0x7FFFFFFF) // Least significant 31 bits
|
|
|
|
// State vector
|
|
typedef struct mt_state_s
|
|
{
|
|
uint64_t mt[NN];
|
|
int mti;
|
|
} mt_state, *mt_state_ptr;
|
|
|
|
// Initializes mt[NN] with a seed
|
|
void twist_init(mt_state_ptr st, uint64_t seed);
|
|
|
|
// Generates a random uint64_t in [0, 2^64-1] interval
|
|
uint64_t twist_uint64(mt_state_ptr st);
|
|
|
|
// Generates a double in [0, 1) interval
|
|
double twist_real(mt_state_ptr st);
|
|
|
|
//
|
|
// Implementation
|
|
//
|
|
|
|
// Initialize a state mt[NN] with a seed
|
|
void twist_init(mt_state_ptr st, uint64_t seed)
|
|
{
|
|
st->mti = NN + 1;
|
|
|
|
st->mt[0] = seed;
|
|
|
|
for (st->mti = 1; st->mti < NN; st->mti++)
|
|
{
|
|
st->mt[st->mti] = (UINT64_C(6364136223846793005) * (st->mt[st->mti - 1] ^ (st->mt[st->mti - 1] >> 62)) + st->mti);
|
|
}
|
|
}
|
|
|
|
// Generates a random uint64_t in [0, 2^64-1] interval
|
|
uint64_t twist_uint64(mt_state_ptr st)
|
|
{
|
|
int i;
|
|
uint64_t x;
|
|
static uint64_t mag01[2] = {UINT64_C(0), MATRIX_A};
|
|
|
|
if (st->mti >= NN)
|
|
{
|
|
// Generate NN words at one time if twist_init() has not been called
|
|
// A default initial seed is used
|
|
if (st->mti == NN + 1)
|
|
twist_init(st, UINT64_C(5489));
|
|
|
|
for (i = 0; i < NN - MM; i++)
|
|
{
|
|
x = (st->mt[i] & UM) | (st->mt[i + 1] & LM);
|
|
st->mt[i] = st->mt[i + MM] ^ (x >> 1) ^ mag01[(int)(x & UINT64_C(1))];
|
|
}
|
|
|
|
for (; i < NN - 1; i++)
|
|
{
|
|
x = (st->mt[i] & UM) | (st->mt[i + 1] & LM);
|
|
st->mt[i] = st->mt[i + (MM - NN)] ^ (x >> 1) ^ mag01[(int)(x & UINT64_C(1))];
|
|
}
|
|
|
|
x = (st->mt[NN - 1] & UM) | (st->mt[0] & LM);
|
|
|
|
st->mt[NN - 1] = st->mt[MM - 1] ^ (x >> 1) ^ mag01[(int)(x & UINT64_C(1))];
|
|
|
|
st->mti = 0;
|
|
}
|
|
|
|
x = st->mt[st->mti++];
|
|
|
|
x ^= (x >> 29) & UINT64_C(0x5555555555555555);
|
|
x ^= (x << 17) & UINT64_C(0x71D67FFFEDA60000);
|
|
x ^= (x << 37) & UINT64_C(0xFFF7EEE000000000);
|
|
x ^= (x >> 43);
|
|
|
|
return x;
|
|
}
|
|
|
|
// Generates a double in [0, 1) interval
|
|
double twist_real(mt_state_ptr st)
|
|
{
|
|
return (twist_uint64(st) >> 11) * (1.0 / 9007199254740992.0);
|
|
}
|