1)ISTORIC

Marin Mersenne s-a nascut la 8 septembrie 1588 in Oize in Maine,Franta si a murit la 1 septembrie 1648 la Paris.
Mersenne a studiat teologia la Sorbonne din 1609 pina in 1611.In 1912 s-a intors la Paris unde a devenit preot la Palatul Regal.
Mersenne Marin(1588 – 1648)a fost calugar,filozof si matematician francez renumit pentru interesul sau in formularea si generarea numerelor prime bazate pe ceea ce astazi sunt cunoscute sub numele de numere Mersenne. Totusi Mersenne nu a fost numai matematician,el a scris despre teoria muzicii si alte subiecte ,a lucrat cu Euclid,Archimedes si alti matematicieni.
El a studiat numerele prime si a incercat sa gaseasca o formula prin care sa reprezinte toate aceste numere.
Marin Mersenne a studiat numerele prime si a incercat sa gaseasca o formula care sa determine toate numerele prime.umere.Este usor sa demonstram ca daca numarul n = 2^p – 1 atunci n trebuie sa fie prim.In 1644 Mersenne afirma in lucrarea “Cugetari fizico – matematice “ ca numarele de forma 2^p – 1 sunt prime pentru n = 2,3,5,7,13,17,19,31,67,127,257.
Multa vreme s-a crezut ca numerele de forma 2^n- 1 sunt numere prime pentru orice n,dar in 1536 Hudalricus Regius a aratat ca 2^11 -1 =2047 nu este numar prim(este egal cu 23*89).In 1603 Pietro Cataldi a verificat corect ca 2^17-1 si 2^19 – 1 sunt ambele numere prime,si de asemenea sunt prime si 23,29,31 si 37.In 1640 Fermat demonstreaza ca Cataldi a gresit in legatura cu numerele 23 si 37,acestea nefiind prime;apoi Euler in 1738 demonstreaza ca acest Cataldi a gresit si in privinta numarului 29.Ceva mai tirziu Euler arata ca in privinta numarului 31 Cataldi a facut o afirmatie corecta,acesta fiind prim.
Mersenne a prefatat lucrarea”Cugetari fizico-matematice”(1644) cu afirmatia prin care sustinea ca numerele de forma 2^n-1 sunt prime pentru n=2,3,5,7,13,17,19,31,67,127 si 257.
Daca un numar de forma 2^n-1(n-numar prim) este prim spunem ca este un numar prim Mersenne.
Era evident ca Mersenne nu a reusit sa testeze toate aceste numere .In 1750,cind Euler a verificat urmatorul numar din lista numerelor Mersenne,2^31-1,a descoperit ca si acest numar este prim.in 1876,Lucas verifica numarul 2^127-1 si ajunge la concluzia ca si acest numar este prim.Sapte ani mai tirziu Pervouchine arata ca numarul 2^61-1 este prim,deci Mersenne a omis acest numar.in jurul anului 1900 Powers demonstreaza ca Mersenne a omis de asemenea alte doua numere:pe 2^89 - 1 si 2^107 - 1 .In sfirsit,in 1947,rangul numerelor Mersenne ,<=258,a fost complet verificat si s-a determinat lista corecta a numerelor Mersenne:
N=2,3,5,7,13,17,19,31,61,89,107 si 127.

2)Definitia Numerelor Mersenne

Fie n un numar intrg prim.Numarul 2^n – 1 se numeste numar Mersenne prim.

Numerele Mersenne sunt numere de forma Mn=2^n – 1,unde n este un numar prim.
Primele numere Mersenne prime sunt 3,7,31,127,8191,131071,… corespunzind indicilor n = 2,3,5,7,13,17,… .
Numerele Mersenne prime au fost studiate pentru proprietatea remarcabila pe care aceste numere o au si anume:fiecare numar prim corespunde exact unui numar patrat perfect.

S-a presupus ca exista o ifinitate de numere Mersenne prime.

Tabela cu numerele Mersenne prime:

Fie M(p) = 2p-1 si P(p) = 2p-1(2p-1). Lista tuturor numerelor Mersenne prime p,cunoscute pentru fiecare M(p) (therefore P(p) este un numar perfect) :
## p
(exponent) digits
in Mp digits
in Pp year discoverer notes
1 2 1 1 ---- ----
2 3 1 2 ---- ----
3 5 2 3 ---- ----
4 7 3 4 ---- ----
5 13 4 8 1456 Anonymous
6 17 6 10 1588 Cataldi
7 19 6 12 1588 Cataldi
8 31 10 19 1772 Euler
9 61 19 37 1883 Pervushin
10 89 27 54 1911 Powers
11 107 33 65 1914 Powers note12 127 39 77 1876 Lucas
13 521 157 314 1952 Robinson
14 607 183 366 1952 Robinson
15 1279 386 770 1952 Robinson
16 2203 664 1327 1952 Robinson
17 2281 687 1373 1952 Robinson
18 3217 969 1937 1957 Riesel
19 4253 1281 2561 1961 Hurwitz
20 4423 1332 2663 1961 Hurwitz
21 9689 2917 5834 1963 Gillies

22 9941 2993 5985 1963 Gillies
23 11213 3376 6751 1963 Gillies
24 19937 6002 12003 1971 Tuckerman [Tuckerman71]25 21701 6533 13066 1978 Noll & Nickel
[NN80]26 23209 6987 13973 1979 Noll "
27 44497 13395 26790 1979 Nelson & Slowinski
[Slowinski79]28 86243 25962 51924 1982 Slowinski [Ewing83]29 110503 33265 66530 1988 Colquitt & Welsh [CW91]30 132049 39751 79502 1983 Slowinski
31 216091 65050 130100 1985 Slowinski
32 756839 227832 455663 1992 Slowinski & Gage
[Peterson92]33 859433 258716 517430 1994 Slowinski & Gage
34 1257787 378632 757263 1996 Slowinski & Gage (web page)35 1398269 420921 841842 1996 Armengaud, Woltman,
et. al. (GIMPS)
(web page)36 2976221 895932 1791864 1997 Spence, Woltman,
et. al. (GIMPS) (web page)37 3021377 909526 1819050 1998 Clarkson, Woltman, Kurowski
et. al. (GIMPS, PrimeNet)
(web page)38 6972593 2098960 4197919 1999 Hajratwala, Woltman, Kurowski
et. al. (GIMPS, PrimeNet) (web page)?? 13466917 4053946 8107892 2001 Cameron, Woltman, Kurowski
et. al. (GIMPS, PrimeNet) (web page)?? 20996011 6320430 12640858 2003 Shafer, Woltman, Kurowski
et. al. (GIMPS, PrimeNet) (web page)?? 24036583 7235733 14471465 2004 Findley, Woltman, Kurowski
et. al. (GIMPS, PrimeNet) (web page)

3)Numere perfecte si citeva teoreme

Multe culture antice au fost preocupate de relatiile care le are un numar cu suma si divizorii sai,de multe ori dinduli-se interpretari mistice.Pe noi ne preocupa gasirea unei singure relatii de acest fel(relatia dintre un numar dat ,suma si divizorii acestui numar).

Definitie:

Un numar intreg pozitiv n se numeste numar perfect daca acesta este egal cu suma tuturor divizorilor sai proprii,mai putin el insusi(excuzindu-l pe n din suma).

De exemplu, 6 este primul numar perfect pentru ca 6=1+2+3.Urmatorul numar perfect este 28=1+2+4+7+14.Urmatoarele doua numere prime sunt 496 si 8128.Daca privim cu atentie aceste patru numere observam ca 2*3,4*7,16*31,64*127 .Numerele au aceeasi forma :2^n-1(2^n-1)(pentru n=2,3,5 si 7).In fiecare caz 2^n – 1 este un numar prim Mersenne.

Teorema1:
Un numar k este patrat perfect daca si numai daca are forma: 2^n-1(2^n-1) si 2^n-1 sunt numere prime.

Teorema2:
Daca 2^n-1 este prim,atunci n este de asemenea prim.

Deci cautarea numerelor Mersenne implica de asemenea cautarea numerelor care sunt patrate perfecte!!!

Daca dorim putem sa folosim si reprezentarea binara a numerelor.Pentru primele patru numere patrate perfecte reprezentarea binara arata astfel:
6=110
28 =11100
496 = 111110000
8128 = 1111111000000

Reprezentarea in forma binara a numerelor este o consecinta a primei teoreme enuntate mai sus.
Cind dorim sa verificam daca un numar este numar prim Mersenne,de obicei,cautam mai intai divizori mai mici.Urmatoarea teorema a lui Euler si Fermat este folositoare in acest sens:

Teorema3:
Presupunem ca p si q sunt doua numere prime.Daca q divide pe Mp=2^p – 1, atunci
Q=+/-1( mod 8) si q = 2kp + 1
pentru k intreg.

Teorema4:
Sa presupunem ca p=3 (mod 4 ) ar fi prim.2p+1 este de asemenea primdaca si numai daca 2p+1 il divide pe Mp.

Teorema5:
Daca adunam cifrele oricarui numar perfect(cu exceptia lui 6),apoi suma cifrelor numarului rezutat,si repetam acest process pina cind ne ramine o singura cifra,acea cifra va fi unu.

Legatura dintre numerele Mersenne si numerele perfecte:

Exista o strinsa relatie intre numerele Mersenne si numerele perfecte.
Numerele perfecte sunt acele numere care sunt egale cu suma divizorilor lor proprii.
Ex:1+ 2 + 3 = 6
1 + 2 + 4 + 7 + 14 = 28
acest lucru a fost demonstrate de Euler :
cind 2^n – 1 este un numar prim,
numarul 2^(n-1) * ( 2^n – 1 ) este patrat perfect.

Numerele Mersenne prime cunoscute :

Mersenne Primes Numere perfecte

2
2 - 1 3 6

3
2 - 1 7 28

5
2 - 1 31 496

7
2 - 1 127 8,128

13
2 - 1 8,191 33,550,336

17
2 - 1 131,071 8,589,869,056

19
2 - 1 524,287 137,438,691,328

31
2 - 1 2,147,483,647

61
2 - 1 2,305,843,009,213,693,952

89
2 - 1 618,970,019,642,690,137,449,562,112

107
2 - 1 162,259,276,829,213,363,391,578,010,288,127

127
2 - 1 170,141,183,460,469,231,731,687,303,715,884,105,727

521
2 - 1

607
2 - 1

1,279
2 - 1

2,203
2 - 1

2,281
2 - 1

3,217
2 - 1

4,253
2 - 1

4,423
2 - 1

9,689
2 - 1

9,941
2 - 1

11,213
2 - 1

19,937
2 - 1

21,701
2 - 1

23,209
2 - 1

44,497
2 - 1

86,243
2 - 1

110,503
2 - 1

132,049
2 - 1

216,091
2 - 1

756,839
2 - 1

859,433
2 - 1

1,257,787
2 - 1

1,398,269
2 - 1

2,976,221
2 - 1

3,021,377
2 - 1

6,972,593
2 - 1

13,466,917
2 - 1

20,996,011
2 - 1

24,036,583
2 - 1

25,964,951
2 - 1

In 1985 cel mai mare numar prim cunoscut era
2^(216,091) – 1.

4) Testele Lucas-Lehmer(Testul de primalitate al unui numar prim)
Generarea numerelor Mersenne prime

Numerele Mersenne prime(asadar si numerele perfecte) au avut la baza urmatoarea teorema:

Testul Lucas-Lehmer :

Pentru p prim,numarul Mersenne 2^p – 1 este prim daca si numai daca 2^p – 1 divide pe S
(p – 1) unde S(n + 1) = S(n)^2 – 2, si S(1) = 4.

(Putem incepe si cu S(1) = 10 )
In pseudocod acest test se scrie astfel:

Lucas_Lehmer_Test(p):
s := 4;
for i from 3 to p do s := s^2 – 2 mod 2^p – 1;
if s = = 0 then
2^p este numar prim
else
2^p – 1 este numar compus

bazele acestei teorii au fost puse de Lucas in 1870 si a fost concretizata intr-un test de catre Lehmer in 1930.
Acest test este ideal pentru calculatoarele binare pentru ca impartirea prin 2^n – 1( in cod binary) poate fi facuta folosind numai rotirea si adunarea.

5)Ipoteze si probleme nerezolvate referitoare la numerele Mersenne

Sunt infinite multe numere Mersenne compuse?
Euler a aratat:

Teorema:
Daca k>1 si p = 4k + 1 este prim,atunci 2*p + 1 este pri daca si numai daca 2^p = 1(mod 2*p + 1).
Deci,daca p = 4k + 3 si 2*p + 1 sunt prime atunci numarul Mersenne 2^p – 1 este compus.

Noua ipoteza referitoare la numere Mersenne:
Bateman,Selfridge si Wagstaff au presupus urmatoarele:

Sa presupunem ca p este un numar natural.Daca sunt indeplinite doua din conditiile de mai jos atunci este indeplinita si a treia conditie:
1. p = 2^k+/-1 sau p = 4^k +/-3
2. 2^p – 1 este prim (evident este un numar Mersenne prim)
3. (2^p + 1) /3 este prim.

Exista mai multe numere Mersenne prime duble?

Daca n = Mp este prim,atunci si Mn este tot prim;sa numin acest numar MMp(un numar “dublu Mersenn”).Intr – adevar fiecare din primele patru numere sunt prime:

MM2 = 2^3 – 1 = 7
MM3 = 2^7 – 1 = 127
MM5 = 2^31 – 1 = 2147483647

MM7 = 2^127 – 1 = 170141183460469231731687303715884105727

6)Cifru Mersenne

a)Algoritm de determinare a cifrului Mersenne:

Acest algoritm este format din doua parti:
-partea recurenta(schimba starea bitilor 19,937 cu cea a succesorilor)
- procedura de extragere a multiplilor
Specificarea algoritmului (schimbarea registrilor):

for i = 0 to 623
temp = first bit of a(i) followed by last 31 bits of a(i+1) ;
a(i) = temp shifted right one bit xor
X'9908B0DF' if temp is odd xor
a(i+397) ;
next i
Implementarea algoritmului in C arata astfel:

for i = 0 to 226
temp = first bit of a(i) followed by last 31 bits of a(i+1) ;
a(i) = temp shifted right one bit xor
X'9908B0DF' if temp is odd xor
a(i+397) ;
next i
for i = 227 to 622
temp = first bit of a(i) followed by last 31 bits of a(i+1) ;
a(i) = temp shifted right one bit xor
X'9908B0DF' if temp is odd xor
a(i-227) ;
next i
temp = first bit of a(623) followed by last 31 bits of a(0) ;
a(i) = temp shifted right one bit xor
X'9908B0DF' if temp is odd xor
a(396) ;
Cifru Mersenne parametrizeaza familia generatorilor ,parametrii aratind astfel:
L=19937
W=32
M=397
A = X'9908B0DF'
Unde W reprezinta lungimea in biti a sirurilor iar L este lungimea in biti a registrului pe care dorim sa-l mutam,care trebuie sa fie exponent pentru numerele Mersenne prime parametrii suplimentari putind fi dervati din urmatoarele:

N = floor(L/W) + 1
R = L mod W
Astfel,dupa ce L a fost impartit la W,R este restul si N este un coeficient intreg.Aceasta ecuatie se aplica numai in cazul in care R este diferit de zero,care este intotdeauna adevarata dar in cazul trivial pentru un numar prim N,N trebuie sa fie prim pentru a fi un exponent Mersene valid.

Daca aplicam parametrizarile de mai sus obtinem urmatorul cod:

for i = 0 to N-M-1
temp = first R bits of x(i) followed by last W-R bits of x(i+1) ;
x(i) = temp shifted right one bit xor
A if temp is odd xor
x(i+M) ;
next i
for i = N-M to N-2
temp = first R bits of x(i) followed by last W-R bits of x(i+1) ;
x(i) = temp shifted right one bit xor
A if temp is odd xor
x(i+M-N) ;
next i
temp = first R bits of x(N-1) followed by last W-R bits of x(0) ;
x(N-1) = temp shifted right one bit xor
A if temp is odd xor
x(M-1) ;

Algoritmul pentru determinarea cifrului Mersenne poate fi ilustrat ca un fel de feedback liniar cu schimbarea registrului:

Un alt algoritm pentru determinarea cifrului Mersenne este urmatorul:

-partea de specificarea a algoritmului
for i = 0 to 350
temp = first 13 bits of a(i) followed by last 19 bits of a(i+1) ;
a(i) = temp shifted right one bit xor
X'E4BD75F5' [or X'CCAB8EE7'] if temp is odd xor
a(i+175) ;
next i

-implementarea partii de specificare
for i = 0 to 175
temp = first 13 bits of a(i) followed by last 19 bits of a(i+1) ;
a(i) = temp shifted right one bit xor
X'E4BD75F5' [or X'CCAB8EE7'] if temp is odd xor
a(i+175) ;
next i
for i = 176 to 349
temp = first 13 bits of a(i) followed by last 19 bits of a(i+1) ;
a(i) = temp shifted right one bit xor
X'E4BD75F5' [or X'CCAB8EE7'] if temp is odd xor
a(i-176) ;
next i
temp = first 13 bits of a(350) followed by last 19 bits of a(0) ;
a(i) = temp shifted right one bit xor
X'E4BD75F5' [or X'CCAB8EE7'] if temp is odd xor
a(174) ;

- functia principala

temp = selected array element ;
temp = temp xor
temp shifted right 11 bits ;
temp = temp xor
((temp shifted left 7 bits) and X'655E5280' [or X'31B6AB00'] ) ;
temp = temp xor
((temp shifted left 15 bits) and X'FFD58000' [or X'FFE50000']) ;
value = temp xor
temp shifted right 17 bits ;

b)Algoritm pentru determinarea cifrului Mersenne:

Acest algoritm a fost proiectat tinind cont de numerosi generatori.Este un generator rapid si evita imultirile si impartirile,acesta beneficiind si de pipelines.
Algoritmul lui Mersenne a fost implementat in C++.Este un algoritm liber si portabil.

Descrierea algoritmului:
-crearea generatorului cu MTRand r;
-generarea numerelor intregi
- Optiuni
-forma implicita
-numar intreg unic
-sir de orice lungime
-este usor de schimbat declaratiile
-exemplu de program
-teste de validare a rezultatelor si performantelor

// The Mersenne Twister is an algorithm for generating random numbers. It
// was designed with consideration of the flaws in various other generators.
// The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
// are far greater. The generator is also fast; it avoids multiplication and
// division, and it benefits from caches and pipelines. For more information

#ifndef MERSENNETWISTER_H
#define MERSENNETWISTER_H

// Not thread safe (unless auto-initialization is avoided and each thread has
// its own MTRand object)

#include <iostream>
#include <limits.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

class MTRand {
// Data
public:
typedef unsigned long uint32; // unsigned integer type, at least 32 bits

enum { N = 624 }; // length of state vector
enum { SAVE = N + 1 }; // length of array for save()

protected:
enum { M = 397 }; // period parameter

uint32 state[N]; // internal state
uint32 *pNext; // next value to get from state
int left; // number of values left before reload needed

//Methods
public:
MTRand( const uint32& oneSeed ); // initialize with a simple uint32
MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or an array
MTRand(); // auto-initialize with /dev/urandom or time() and clock()

// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
// values together, otherwise the generator state can be learned after
// reading 624 consecutive values.

// Access to 32-bit random numbers
double rand(); // real number in [0,1]
double rand( const double& n ); // real number in [0,n]
double randExc(); // real number in [0,1)
double randExc( const double& n ); // real number in [0,n)
double randDblExc(); // real number in (0,1)
double randDblExc( const double& n ); // real number in (0,n)
uint32 randInt(); // integer in [0,2^32-1]
uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32
double operator()() { return rand(); } // same as rand()

// Access to 53-bit random numbers (capacity of IEEE double precision)
double rand53(); // real number in [0,1)

// Access to nonuniform random number distributions
double randNorm( const double& mean = 0.0, const double& variance = 0.0 );

// Re-seeding functions with same behavior as initializers
void seed( const uint32 oneSeed );
void seed( uint32 *const bigSeed, const uint32 seedLength = N );
void seed();

// Saving and loading generator state
void save( uint32* saveArray ) const; // to array of size SAVE
void load( uint32 *const loadArray ); // from such array
friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
friend std::istream& operator>>( std::istream& is, MTRand& mtrand );

protected:
void initialize( const uint32 oneSeed );
void reload();
uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
uint32 mixBits( const uint32& u, const uint32& v ) const
{ return hiBit(u) | loBits(v); }
uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
{ return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
static uint32 hash( time_t t, clock_t c );
};

inline MTRand::MTRand( const uint32& oneSeed )
{ seed(oneSeed); }

inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
{ seed(bigSeed,seedLength); }

inline MTRand::MTRand()
{ seed(); }

inline double MTRand::rand()
{ return double(randInt()) * (1.0/4294967295.0); }

inline double MTRand::rand( const double& n )
{ return rand() * n; }

inline double MTRand::randExc()
{ return double(randInt()) * (1.0/4294967296.0); }

inline double MTRand::randExc( const double& n )
{ return randExc() * n; }

inline double MTRand::randDblExc()
{ return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }

inline double MTRand::randDblExc( const double& n )
{ return randDblExc() * n; }

inline double MTRand::rand53()
{
uint32 a = randInt() >> 5, b = randInt() >> 6;
return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada
}

inline double MTRand::randNorm( const double& mean, const double& variance )
{
// Return a real number from a normal (Gaussian) distribution with given
// mean and variance by Box-Muller method
double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
return mean + r * cos(phi);
}

inline MTRand::uint32 MTRand::randInt()
{
// Pull a 32-bit integer from the generator state
// Every other access function simply transforms the numbers extracted here

if( left == 0 ) reload();
--left;

register uint32 s1;
s1 = *pNext++;
s1 ^= (s1 >> 11);
s1 ^= (s1 << 7) & 0x9d2c5680UL;
s1 ^= (s1 << 15) & 0xefc60000UL;
return ( s1 ^ (s1 >> 18) );
}

inline MTRand::uint32 MTRand::randInt( const uint32& n )
{
// Find which bits are used in n
// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
uint32 used = n;
used |= used >> 1;
used |= used >> 2;
used |= used >> 4;
used |= used >> 8;
used |= used >> 16;

// Draw numbers until one is found in [0,n]
uint32 i;
do
i = randInt() & used; // toss unused bits to shorten search
while( i > n );
return i;
}

inline void MTRand::seed( const uint32 oneSeed )
{
// Seed the generator with a simple uint32
initialize(oneSeed);
reload();
}

inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
{
// Seed the generator with an array of uint32's
// There are 2^19937-1 possible initial states. This function allows
// all of those to be accessed by providing at least 19937 bits (with a
// default seed length of N = 624 uint32's). Any bits above the lower 32
// in each element are discarded.
// Just call seed() if you want to get array from /dev/urandom
initialize(19650218UL);
register int i = 1;
register uint32 j = 0;
register int k = ( N > seedLength ? N : seedLength );
for( ; k; --k )
{
state[i] =
state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
state[i] &= 0xffffffffUL;
++i; ++j;
if( i >= N ) { state[0] = state[N-1]; i = 1; }
if( j >= seedLength ) j = 0;
}
for( k = N - 1; k; --k )
{
state[i] =
state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
state[i] -= i;
state[i] &= 0xffffffffUL;
++i;
if( i >= N ) { state[0] = state[N-1]; i = 1; }
}
state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
reload();
}

inline void MTRand::seed()
{
// Seed the generator with an array from /dev/urandom if available
// Otherwise use a hash of time() and clock() values

// First try getting an array from /dev/urandom
FILE* urandom = fopen( "/dev/urandom", "rb" );
if( urandom )
{
uint32 bigSeed[N];
register uint32 *s = bigSeed;
register int i = N;
register bool success = true;
while( success && i-- )
success = fread( s++, sizeof(uint32), 1, urandom );
fclose(urandom);
if( success ) { seed( bigSeed, N ); return; }
}

// Was not successful, so use time() and clock() instead
seed( hash( time(NULL), clock() ) );
}

inline void MTRand::initialize( const uint32 seed )
{
// Initialize generator state with seed
// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
// In previous versions, most significant bits (MSBs) of the seed affect
// only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
register uint32 *s = state;
register uint32 *r = state;
register int i = 1;
*s++ = seed & 0xffffffffUL;
for( ; i < N; ++i )
{
*s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
r++;
}
}

inline void MTRand::reload()
{
// Generate N new values in state
// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
register uint32 *p = state;
register int i;
for( i = N - M; i--; ++p )
*p = twist( p[M], p[0], p[1] );
for( i = M; --i; ++p )
*p = twist( p[M-N], p[0], p[1] );
*p = twist( p[M-N], p[0], state[0] );

left = N, pNext = state;
}

inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
{
// Get a uint32 from t and c
// Better than uint32(x) in case x is floating point in [0,1]
// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)

static uint32 differ = 0; // guarantee time-based seeds will change

uint32 h1 = 0;
unsigned char *p = (unsigned char *) &t;
for( size_t i = 0; i < sizeof(t); ++i )
{
h1 *= UCHAR_MAX + 2U;
h1 += p[i];
}
uint32 h2 = 0;
p = (unsigned char *) &c;
for( size_t j = 0; j < sizeof(c); ++j )
{
h2 *= UCHAR_MAX + 2U;
h2 += p[j];
}
return ( h1 + differ++ ) ^ h2;
}

inline void MTRand::save( uint32* saveArray ) const
{
register uint32 *sa = saveArray;
register const uint32 *s = state;
register int i = N;
for( ; i--; *sa++ = *s++ ) {}
*sa = left;
}

inline void MTRand::load( uint32 *const loadArray )
{
register uint32 *s = state;
register uint32 *la = loadArray;
register int i = N;
for( ; i--; *s++ = *la++ ) {}
left = *la;
pNext = &state[N-left];
}

inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand )
{
register const MTRand::uint32 *s = mtrand.state;
register int i = mtrand.N;
for( ; i--; os << *s++ << "\t" ) {}
return os << mtrand.left;
}

inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
{
register MTRand::uint32 *s = mtrand.state;
register int i = mtrand.N;
for( ; i--; is >> *s++ ) {}
is >> mtrand.left;
mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
return is;
}

#endif // MERSENNETWISTER_H

// Change log:
//
// v0.1 - First release on 15 May 2000
// - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
// - Translated from C to C++
// - Made completely ANSI compliant
// - Designed convenient interface for initialization, seeding, and
// obtaining numbers in default or user-defined ranges
// - Added automatic seeding from /dev/urandom or time() and clock()
// - Provided functions for saving and loading generator state
//
// v0.2 - Fixed bug which reloaded generator one step too late
//
// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
//
// v0.4 - Removed trailing newline in saved generator format to be consistent
// with output format of built-in types
//
// v0.5 - Improved portability by replacing static const int's with enum's and
// clarifying return values in seed(); suggested by Eric Heimburg
// - Removed MAXINT constant; use 0xffffffffUL instead
//
// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
// - Changed integer [0,n] generator to give better uniformity
//
// v0.7 - Fixed operator precedence ambiguity in reload()
// - Added access for real numbers in (0,1) and (0,n)
//
// v0.8 - Included time.h header to properly support time_t and clock_t
//
// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
// - Allowed for seeding with arrays of any length
// - Added access for real numbers in [0,1) with 53-bit resolution
// - Added access for real numbers from normal (Gaussian) distributions
// - Increased overall speed by optimizing twist()
// - Doubled speed of integer [0,n] generation
// - Fixed out-of-range number generation on 64-bit machines
// - Improved portability by substituting literal constants for long enum's
// - Changed license from GNU LGPL to BSD