https://github.com/ntamas/plfit
Tip revision: 651adf2bdb8d6b468f16c930898b8a5c2597a253 authored by Tamas Nepusz on 12 March 2024, 13:46:47 UTC
chore: bumped version to 0.9.6
chore: bumped version to 0.9.6
Tip revision: 651adf2
mt.c
/* mt.c
*
* Mersenne Twister random number generator, based on the implementation of
* Michael Brundage (which has been placed in the public domain).
*
* Author: Tamas Nepusz (original by Michael Brundage)
*
* See the following URL for the original implementation:
* http://www.qbrundage.com/michaelb/pubs/essays/random_number_generation.html
*
* This file has been placed in the public domain.
*/
#include <stdlib.h>
#include "plfit_mt.h"
static uint16_t get_random_uint16() {
return rand();
}
void plfit_mt_init(plfit_mt_rng_t* rng) {
plfit_mt_init_from_rng(rng, 0);
}
void plfit_mt_init_from_rng(plfit_mt_rng_t* rng, plfit_mt_rng_t* seeder) {
int i;
if (seeder == 0) {
for (i = 0; i < PLFIT_MT_LEN; i++) {
/* RAND_MAX is guaranteed to be at least 32767, so we can use two
* calls to rand() to produce a random 32-bit number */
rng->mt_buffer[i] = (((uint32_t) get_random_uint16()) << 16) + get_random_uint16();
}
} else {
for (i = 0; i < PLFIT_MT_LEN; i++) {
rng->mt_buffer[i] = plfit_mt_random(seeder);
}
}
rng->mt_index = 0;
}
#define MT_IA 397
#define MT_IB (PLFIT_MT_LEN - MT_IA)
#define UPPER_MASK 0x80000000
#define LOWER_MASK 0x7FFFFFFF
#define MATRIX_A 0x9908B0DF
#define TWIST(b,i,j) ((b)[i] & UPPER_MASK) | ((b)[j] & LOWER_MASK)
#define MAGIC(s) (((s)&1)*MATRIX_A)
uint32_t plfit_mt_random(plfit_mt_rng_t* rng) {
uint32_t * b = rng->mt_buffer;
int idx = rng->mt_index;
uint32_t s;
int i;
if (idx == PLFIT_MT_LEN * sizeof(uint32_t)) {
idx = 0;
i = 0;
for (; i < MT_IB; i++) {
s = TWIST(b, i, i+1);
b[i] = b[i + MT_IA] ^ (s >> 1) ^ MAGIC(s);
}
for (; i < PLFIT_MT_LEN-1; i++) {
s = TWIST(b, i, i+1);
b[i] = b[i - MT_IB] ^ (s >> 1) ^ MAGIC(s);
}
s = TWIST(b, PLFIT_MT_LEN-1, 0);
b[PLFIT_MT_LEN-1] = b[MT_IA-1] ^ (s >> 1) ^ MAGIC(s);
}
rng->mt_index = idx + sizeof(uint32_t);
return *(uint32_t *)((unsigned char *)b + idx);
/*
Matsumoto and Nishimura additionally confound the bits returned to the caller
but this doesn't increase the randomness, and slows down the generator by
as much as 25%. So I omit these operations here.
r ^= (r >> 11);
r ^= (r << 7) & 0x9D2C5680;
r ^= (r << 15) & 0xEFC60000;
r ^= (r >> 18);
*/
}
double plfit_mt_uniform_01(plfit_mt_rng_t* rng) {
return ((double)plfit_mt_random(rng)) / PLFIT_MT_RAND_MAX;
}