Permalink
Browse files

using Mersenne Twister instead of built-in RNG

  • Loading branch information...
1 parent c95e8b2 commit dd9507318bdd2a38131fff7825d4b71e400f3e39 @ntamas committed Mar 18, 2012
Showing with 177 additions and 8 deletions.
  1. +1 −1 src/CMakeLists.txt
  2. +75 −0 src/mt.c
  3. +87 −0 src/mt.h
  4. +4 −2 src/plgen.c
  5. +3 −3 src/sampling.c
  6. +3 −1 src/sampling.h
  7. +4 −1 test/test_sampling.c
View
@@ -1,6 +1,6 @@
configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config.h.in ${CMAKE_CURRENT_BINARY_DIR}/config.h)
-set(PLFIT_CORE_SRCS lbfgs.c error.c kolmogorov.c plfit.c sampling.c zeta.c)
+set(PLFIT_CORE_SRCS lbfgs.c error.c kolmogorov.c mt.c plfit.c sampling.c zeta.c)
include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
View
@@ -0,0 +1,75 @@
+/* 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 "mt.h"
+
+void mt_init(mt_rng_t* rng) {
+ int i;
+ for (i = 0; i < 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] = (rand() << 16) + rand();
+ }
+ rng->mt_index = 0;
+}
+
+#define MT_IA 397
+#define MT_IB (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 mt_random(mt_rng_t* rng) {
+ uint32_t * b = rng->mt_buffer;
+ int idx = rng->mt_index;
+ uint32_t s;
+ int i;
+
+ if (idx == 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 < MT_LEN-1; i++) {
+ s = TWIST(b, i, i+1);
+ b[i] = b[i - MT_IB] ^ (s >> 1) ^ MAGIC(s);
+ }
+
+ s = TWIST(b, MT_LEN-1, 0);
+ b[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 mt_uniform_01(mt_rng_t* rng) {
+ return ((double)mt_random(rng)) / MT_RAND_MAX;
+}
+
View
@@ -0,0 +1,87 @@
+/* mt.h
+ *
+ * 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.
+ */
+
+#ifndef __MT_H__
+#define __MT_H__
+
+#ifdef _MSC_VER
+# define uint32_t __int32;
+#else
+# include <stdint.h>
+#endif
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+__BEGIN_DECLS
+
+#define MT_LEN 624
+
+/**
+ * \def MT_RAND_MAX
+ *
+ * The maximum random number that \c mt_random() can generate.
+ */
+#define MT_RAND_MAX 0xFFFFFFFF
+
+/**
+ * Struct that stores the internal state of a Mersenne Twister random number
+ * generator.
+ */
+typedef struct {
+ int mt_index;
+ uint32_t mt_buffer[MT_LEN];
+} mt_rng_t;
+
+/**
+ * \brief Initializes a Mersenne Twister random number generator.
+ *
+ * The random number generator is seeded with random 32-bit numbers obtained
+ * from the \em built-in random number generator using consecutive calls to
+ * \c rand().
+ *
+ * \param rng the random number generator to initialize
+ */
+void mt_init(mt_rng_t* rng);
+
+/**
+ * \brief Returns the next 32-bit random number from the given Mersenne Twister
+ * random number generator.
+ *
+ * \param rng the random number generator to use
+ * \return the next 32-bit random number from the generator
+ */
+uint32_t mt_random(mt_rng_t* rng);
+
+/**
+ * \brief Returns a uniformly distributed double from the interval [0;1)
+ * based on the next value of the given Mersenne Twister random number
+ * generator.
+ *
+ * \param rng the random number generator to use
+ * \return a uniformly distributed random number from the interval [0;1)
+ */
+double mt_uniform_01(mt_rng_t* rng);
+
+__END_DECLS
+
+#endif
+
+
View
@@ -40,6 +40,7 @@ typedef struct _cmd_options_t {
} cmd_options_t;
cmd_options_t opts;
+mt_rng_t rng;
void show_version(FILE* f) {
fprintf(f, "plgen " PLFIT_VERSION_STRING "\n");
@@ -177,7 +178,7 @@ int sample_discrete() {
/* Sampling */
while (opts.num_samples > 0) {
n = opts.num_samples > BLOCK_SIZE ? BLOCK_SIZE : opts.num_samples;
- plfit_walker_alias_sampler_sample(&sampler, samples, n);
+ plfit_walker_alias_sampler_sample(&sampler, samples, n, &rng);
for (i = 0; i < n; i++) {
printf("%ld\n", (long int)(samples[i] + opts.offset));
@@ -203,7 +204,7 @@ int sample_continuous() {
}
for (i = 0; i < opts.num_samples; i++) {
- u = rand() / ((double)RAND_MAX);
+ u = mt_uniform_01(&rng);
printf("%.8f\n", pow(u, gamma));
}
@@ -231,6 +232,7 @@ int main(int argc, char* argv[]) {
retval = 4;
} else {
srand(opts.use_seed ? opts.seed : time(0));
+ mt_init(&rng);
retval = opts.continuous ? sample_continuous() : sample_discrete();
}
}
View
@@ -115,15 +115,15 @@ void plfit_walker_alias_sampler_destroy(plfit_walker_alias_sampler_t* sampler) {
int plfit_walker_alias_sampler_sample(const plfit_walker_alias_sampler_t* sampler,
- long int *xs, size_t n) {
+ long int *xs, size_t n, mt_rng_t* rng) {
double u;
long int j;
long int *x;
x = xs;
while (n > 0) {
- u = rand() / ((double)RAND_MAX);
- j = rand() % sampler->num_bins;
+ u = mt_uniform_01(rng);
+ j = mt_random(rng) % sampler->num_bins;
*x = (u < sampler->probs[j]) ? j : sampler->indexes[j];
n--; x++;
}
View
@@ -21,6 +21,7 @@
#define __SAMPLING_H__
#include <stdlib.h>
+#include "mt.h"
#undef __BEGIN_DECLS
#undef __END_DECLS
@@ -70,10 +71,11 @@ void plfit_walker_alias_sampler_destroy(plfit_walker_alias_sampler_t* sampler);
* \param xs pointer to an array where the sampled items should be
* written
* \param n the number of samples to draw
+ * \param rng the Mersenne Twister random number generator to use
* \return error code
*/
int plfit_walker_alias_sampler_sample(const plfit_walker_alias_sampler_t* sampler,
- long int* xs, size_t n);
+ long int* xs, size_t n, mt_rng_t* rng);
__END_DECLS
View
@@ -30,6 +30,9 @@ int test_sampling() {
long int data[NUM_SAMPLES];
long int hist[NUM_ITEMS];
long int i, j, max_hist;
+ mt_rng_t rng;
+
+ mt_init(&rng);
/* Create the sampler, sample, destroy */
for (i = 0; i < NUM_ITEMS; i++) {
@@ -39,7 +42,7 @@ int test_sampling() {
if (plfit_walker_alias_sampler_init(&sampler, probs, NUM_ITEMS)) {
return 1;
}
- if (plfit_walker_alias_sampler_sample(&sampler, data, NUM_SAMPLES)) {
+ if (plfit_walker_alias_sampler_sample(&sampler, data, NUM_SAMPLES, &rng)) {
return 2;
}
plfit_walker_alias_sampler_destroy(&sampler);

0 comments on commit dd95073

Please sign in to comment.