We wish it's quadratic. It's actually exponential in complexity. However there
are tricks that you can do, but you can't be sure that your result is optimal.
@Georges: I do not think we complicated this problem. I did however push this conversation into something collaborative, competitive and exciting at the same time.
@ZeRaW: Yes, I do have about 10 hours of free time a day, be it a work day or not.
Anyway, here's the 64-bit Mersene Twister ported to C#. I ported it based on a 32 bit implementation in C#, and matched it with the C implementation I am using. You're welcome.
main.cs
using System;
namespace MerseneTwisterTester
{
class Program
{
static void Main(string[] args)
{
mt19937.sgenrand(255);
for(int i = 0; i < 10; i++)
Console.WriteLine(mt19937.genrand());
}
}
}
mt19937.cs
using System;
public class mt19937
{
/* Period parameters */
private const int N = 312;
private const int M = 156;
private const ulong MATRIX_A = 0xB5026F5AA96619E9; /* constant vector a */
private const ulong UPPER_MASK = 0xFFFFFFFF80000000; /* most significant w-r bits */
private const ulong LOWER_MASK = 0x7FFFFFFF; /* least significant r bits */
/* Tempering parameters */
private const ulong TEMPERING_MASK_B = 0x9d2c5680;
private const ulong TEMPERING_MASK_C = 0xefc60000;
private static ulong TEMPERING_SHIFT_U(ulong y) { return (y >> 29); }
private static ulong TEMPERING_SHIFT_S(ulong y) { return (y << 17); }
private static ulong TEMPERING_SHIFT_T(ulong y) { return (y << 37); }
private static ulong TEMPERING_SHIFT_L(ulong y) { return (y >> 43); }
//static unsigned long mt[N]; /* the array for the state vector */
private static ulong[] mt = new ulong[N];
// static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
private static uint mti = N + 1; /* mti==N+1 means mt[N] is not initialized */
/* initializing the array with a NONZERO seed */
public static void sgenrand(ulong seed)
{
/* setting initial seeds to mt[N] using */
/* the generator Line 25 of Table 1 in */
/* [KNUTH 1981, The Art of Computer Programming */
/* Vol. 2 (2nd Ed.), pp102] */
mt[0] = seed;
for(mti = 1; mti < N; mti++) {
mt[mti] = (6364136223846793005 * (mt[mti-1] ^ (mt[mti-1] >> 62)) + mti);
}
}
private static ulong[/* 2 */] mag01 = { 0x0, MATRIX_A };
/* generating reals */
/* unsigned long */
/* for integer generation */
public static ulong genrand()
{
ulong y;
/* mag01[x] = x * MATRIX_A for x=0,1 */
if(mti >= N) /* generate N words at one time */ {
short kk;
if(mti == N + 1) /* if sgenrand() has not been called, */ {
sgenrand(5489); /* a default initial seed is used */
}
for(kk = 0; kk < N - M; kk++) {
y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
}
for(; kk < N - 1; kk++) {
y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
}
y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1];
mti = 0;
}
y = mt[mti++];
y ^= (y >> 29) & 0x5555555555555555;
y ^= (y << 17) & 0x71D67FFFEDA60000;
y ^= (y << 37) & 0xFFF7EEE000000000;
y ^= (y >> 43);
return y;
}
}
For the rare C/C++ developer, here we go:
main.cpp
#include <iostream>
using namespace std;
extern "C" { // when using C, just include the file.
#include "mt64.h"
}
int main(){
init_genrand64(255);
for(int i = 0; i < 10; i++)
cout << genrand64_int64() << endl;
return 0;
}
Header file, summary for user and for compiler (yes, we fluff our compiler that much)
mt64.h
/* initializes mt[NN] with a seed */
void init_genrand64(unsigned long long seed);
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
void init_by_array64(unsigned long long init_key[],
unsigned long long key_length);
/* generates a random number on [0, 2^64-1]-interval */
unsigned long long genrand64_int64(void);
/* generates a random number on [0, 2^63-1]-interval */
long long genrand64_int63(void);
/* generates a random number on [0,1]-real-interval */
double genrand64_real1(void);
/* generates a random number on [0,1)-real-interval */
double genrand64_real2(void);
/* generates a random number on (0,1)-real-interval */
double genrand64_real3(void);
and finally:
mt19937-64.c. Yes, I mix C and C++ in source code without fear, bitch.
#include <stdio.h>
#include "mt64.h"
#define NN 312
#define MM 156
#define MATRIX_A 0xB5026F5AA96619E9ULL
#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */
#define LM 0x7FFFFFFFULL /* Least significant 31 bits */
/* The array for the state vector */
static unsigned long long mt[NN];
/* mti==NN+1 means mt[NN] is not initialized */
static int mti=NN+1;
/* initializes mt[NN] with a seed */
void init_genrand64(unsigned long long seed)
{
mt[0] = seed;
for (mti=1; mti<NN; mti++)
mt[mti] = (6364136223846793005ULL * (mt[mti-1] ^ (mt[mti-1] >> 62)) + mti);
}
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
void init_by_array64(unsigned long long init_key[],
unsigned long long key_length)
{
unsigned long long i, j, k;
init_genrand64(19650218ULL);
i=1; j=0;
k = (NN>key_length ? NN : key_length);
for (; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 3935559000370003845ULL))
+ init_key[j] + j; /* non linear */
i++; j++;
if (i>=NN) { mt[0] = mt[NN-1]; i=1; }
if (j>=key_length) j=0;
}
for (k=NN-1; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 2862933555777941757ULL))
- i; /* non linear */
i++;
if (i>=NN) { mt[0] = mt[NN-1]; i=1; }
}
mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */
}
/* generates a random number on [0, 2^64-1]-interval */
unsigned long long genrand64_int64(void)
{
int i;
unsigned long long x;
static unsigned long long mag01[2]={0ULL, MATRIX_A};
if (mti >= NN) { /* generate NN words at one time */
/* if init_genrand64() has not been called, */
/* a default initial seed is used */
if (mti == NN+1)
init_genrand64(5489ULL);
for (i=0;i<NN-MM;i++) {
x = (mt[i]&UM)|(mt[i+1]&LM);
mt[i] = mt[i+MM] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
}
for (;i<NN-1;i++) {
x = (mt[i]&UM)|(mt[i+1]&LM);
mt[i] = mt[i+(MM-NN)] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
}
x = (mt[NN-1]&UM)|(mt[0]&LM);
mt[NN-1] = mt[MM-1] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
mti = 0;
}
x = mt[mti++];
x ^= (x >> 29) & 0x5555555555555555ULL;
x ^= (x << 17) & 0x71D67FFFEDA60000ULL;
x ^= (x << 37) & 0xFFF7EEE000000000ULL;
x ^= (x >> 43);
return x;
}
/* generates a random number on [0, 2^63-1]-interval */
long long genrand64_int63(void)
{
return (long long)(genrand64_int64() >> 1);
}
/* generates a random number on [0,1]-real-interval */
double genrand64_real1(void)
{
return (genrand64_int64() >> 11) * (1.0/9007199254740991.0);
}
/* generates a random number on [0,1)-real-interval */
double genrand64_real2(void)
{
return (genrand64_int64() >> 11) * (1.0/9007199254740992.0);
}
/* generates a random number on (0,1)-real-interval */
double genrand64_real3(void)
{
return ((genrand64_int64() >> 12) + 0.5) * (1.0/4503599627370496.0);
}
No matter what you use, this output will signify your success.
3220586997909315655
3303451203970382242
11896436706893466529
8960318650144385956
4679212705770455613
15567843309247195414
6961994097256010468
10807484256991480663
11890420171946432686
15474158341220671739
So far I am taking the random number and masking it with 0xFFFF (though the output above is the raw 64bit numbers). 0xFFFF is even larger than the largest number returned by the C/C++ stdlib random function (0x7FFF) (which is pretty cute now that I compare).
Gladiators, begin!