mt19937-period.cpp

#include <cstdio>
#include <cstdint>
#include <ctime>
#include <array>

typedef uint32_t u32;
const int P = 19937, N = 624, M = 397, W = 32, R = 31;
const u32 MATRIX_A = 0x9908b0df;
const u32 LMASK = (1u << R) - 1;
const u32 UMASK = ~LMASK;

u32 mult_a(u32 v, u32 a) {
	return (v >> 1) ^ ((v & 1) ? a : 0);
}

struct State {
	u32 vector[N];
};

void phi_inverse(State *s, u32 *x, u32 a) {
	for (int i = P; i >= N; i --) {
		u32 y = (x[i] ^ x[i-N+M] ^ ((x[i-N+1] & 1) ? a : 0)) << 1;
		x[i-N+1] = (x[i-N+1] & UMASK) | (y & LMASK) | (x[i-N+1] & 1);
		x[i-N] = (y & UMASK) | (x[i-N] & LMASK);
	}
	memcpy(s->vector, &x[0], N*sizeof(u32));
}

void generate(u32 *bits, State *s, u32 a) {
	bits[0] = 0;
	u32 *vector = s->vector;
	int index = 1;
	int i;
	for (i = 1; i < N; i+=2) {
		bits[index] = vector[i] & 1;
		index += 1;
	}
	while (index <= P) {
		int k = 0;
		for (; k < N-M; k ++) {
			vector[k] = vector[k+M] ^ mult_a((vector[k] & UMASK) | (vector[k+1] & LMASK), a);
		}
		for (; k < N-1; k ++) {
			vector[k] = vector[k+M-N] ^ mult_a((vector[k] & UMASK) | (vector[k+1] & LMASK), a);
		}
		vector[N-1] = vector[M-1] ^ mult_a((vector[N-1] & UMASK) | (vector[0] & LMASK), a);
		if (N & 1) {
			i = 1 - (i & 1);
		} else {
			i = 1;
		}
		while (i < N && index <= P) {
			bits[index] = vector[i] & 1;
			index ++;
			i += 2;
		}
	}
}

bool test(u32 a) {
	std::array<u32, P+1> initial = {}, x;
	State s;
	initial[2] = 1;
	x = initial;
	for (int i = 0; i < P; i ++) {
		phi_inverse(&s, &x[0], a);
		generate(&x[0], &s, a);
	}
	return x == initial;
}

u32 nextrnd(u32 x) {
	x = x ^ (x << 13);
	x = x ^ (x >> 17);
	x = x ^ (x << 15);
	return x;
}

void search(u32 start, u32 max) {
	u32 x = start;
	for (u32 i = 0; i < max; i++) {
		if ((i % 128) == 0) printf("%d %u\n", i, time(0));
		if (test(x)) {
			printf("found %08x time=%u\n", x, time(0));
		}
		x = nextrnd(x);
	}
}

int main(void) {
	setvbuf(stdout, 0, _IONBF, 0);
	printf("time=%u\n", time(0));
#ifdef _OPENMP
	printf("OPENMP is enabled.\n");
#endif

#pragma omp parallel
#pragma omp sections
	{
#pragma omp section
	{
		search(1, 0x80000000);
	}
#pragma omp section
	{
		// 0x27070f90 = nextrnd^{0x80000000}(1)
		search(0x27070f90, 0x80000000);
	}
	}
	puts("finish");
	return 0;
}