// Taken from
// http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.c

#define _XOPEN_SOURCE 700

#include <assert.h>
#include <inttypes.h>
#include <stdbool.h>
#include <stdio.h>


#define N 624U
static const uint16_t M = 397;
static const uint32_t MATRIX_A   = 0x9908b0dfUL;
static const uint32_t UPPER_MASK = 0x80000000UL;
static const uint32_t LOWER_MASK = 0x7fffffffUL;

static uint32_t mt[N];
static uint16_t mti;
static bool initialized = false;

static void
init(uint32_t s) {
	initialized = true;
	mt[0] = s & 0xffffffffUL;
	for (mti = 1; mti < N; mti++) {
		mt[mti] = (1812433253 * (mt[mti - 1] ^ (mt[mti - 1] >> 30))) + mti;
	}
}

/// Generates a random number on the [ 0, 0xffffffff ]-interval.
static uint32_t
generate(void) {
	assert(initialized == true);
	static uint32_t mag01[2] = { 0, MATRIX_A };

	if (mti >= N) {
		uint16_t kk;

		for (kk = 0; kk < N - M; kk++) {
			const uint32_t y0 =
				(mt[kk]     & UPPER_MASK) |
				(mt[kk + 1] & LOWER_MASK);
			mt[kk] = mt[kk + M] ^ (y0 >> 1) ^ mag01[y0 & 1U];
		}

		for (; kk < N - 1; kk++) {
			const uint32_t y1 =
				(mt[kk]     & UPPER_MASK) |
				(mt[kk + 1] & LOWER_MASK);
			mt[kk] = mt[kk + (M - N)] ^ (y1 >> 1) ^ mag01[y1 & 1U];
		}

		const uint32_t y2 =
			(mt[N - 1] & UPPER_MASK) |
			(mt[0]     & LOWER_MASK);
		mt[N - 1] = mt[M - 1] ^ (y2 >> 1) ^ mag01[y2 & 1U];
		mti = 0;
	}

	uint32_t y = mt[mti++];

	y ^= (y >> 11);
	y ^= (y << 7)  & 0x9d2c5680UL;
	y ^= (y << 15) & 0xefc60000UL;
	y ^= (y >> 18);

	return y;
}


int
main(void) {
	init(2222);

	for (uint16_t i = 0; i < 10000; i++) {
		assert(
			printf("%" PRIu32 "\n", generate()) > 0
		);
	}

	return 0;
}