// Mersenne Twister random number generator
// Based on code by Makoto Matsumoto, Takuji Nishimura, Shawn Cokus, and Richard J. Wagner

// 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
// see the inventors' web page at 
//
//		http://www.math.keio.ac.jp/~matumoto/emt.html

// Reference
// M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
// Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
// Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.

// Copyright (C) 2000  Matthew Bellew
// 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
// 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

// The original code included the following notice:
//
//     Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura.
//     When you use this, send an email to: matumoto@math.keio.ac.jp
//     with an appropriate reference to your work.
//
// It would be nice to CC: rjwagner@writeme.com, cokus@math.washington.edu, and matthewb@westside.com
// when you write.



//-------------------------------------
// MMXRand.h
//
//	written by Matthew Bellew (July 2000)
//		http://www.bellew.net/rand/MMXRand.html

//
// Derived from Richard Wagner's implementation
//		http://www-personal.engin.umich.edu/~wagnerr/MersenneTwister.html
// of MersenneTwister
//		http://www.math.keio.ac.jp/~matumoto/emt.html
//



#ifndef _MMXRand_
#define _MMXRand_

class MMXTwister
	{
	protected:
		enum
			{
			N = 624,
			M = 397,
			MAGIC = 0x9908b0df
			};
		unsigned long m_state[N];
		unsigned long m_hashed[N]; // compute all at once using mmx
		long m_next;

	protected:
		static unsigned long Hash(unsigned long x)
			{
			x ^= (x >> 11);
			x ^= (x <<  7) & 0x9d2c5680U;
			x ^= (x << 15) & 0xefc60000U;
			return ( x ^ (x >> 18) );
			}
		unsigned long hiBit( unsigned long u )  const { return u & 0x80000000U; }
		unsigned long loBit( unsigned long u )  const { return u & 0x00000001U; }
		unsigned long loBits( unsigned long u ) const { return u & 0x7fffffffU; }
		unsigned long mixBits( unsigned long u, unsigned long v ) const
			{ 
			return hiBit(u) | loBits(v); 
			}
		unsigned long twist( unsigned long m, unsigned long s0, unsigned long s1 ) const
			{ 
			return m ^ (mixBits(s0,s1)>>1) ^ (loBit(s1) ? MAGIC : 0U); 
			}

		void DoTheTwist();
		void HashIt();

		void DoTheTwistC();
		void DoTheTwistMMX();
		void HashItC();
		void HashItASM();
		void HashItMMX();

		virtual void Reload()
			{
			DoTheTwistMMX();
			HashItMMX();
			m_next = N;
			}

	public:
		MMXTwister()
			{
			Seed();	// calls Reload()
			}
		void Seed(unsigned long s=0);

		unsigned long randInt()
			{
			if (--m_next < 0)
				Reload();
			return m_hashed[m_next];
			}

		double rand()
			{ return double(randInt()) * 2.3283064370807974e-10; }

		double rand( const double n )
			{ return rand() * n; }

		double randExc()
			{ return double(randInt()) * 2.3283064365386963e-10; }

		double randExc( const double n )
			{ return randExc() * n; }

		unsigned long randInt( unsigned long n )
			{ return int( randExc() * (n+1) ); }
			// Is there a faster/better way?

		unsigned long operator()()
			{
			return randInt();
			}

		unsigned long operator()(unsigned long n)
			{
			return randInt(n);
			}
	};
#endif



//-------------------------------------
// MMXRand.cpp
//
//	written by Matthew Bellew (July 2000)
//		http://www.bellew.net/rand/MMXRand.html

//
// Derived from Richard Wagner's implementation
//		http://www-personal.engin.umich.edu/~wagnerr/MersenneTwister.html
// of MersenneTwister
//		http://www.math.keio.ac.jp/~matumoto/emt.html
//


#ifndef _MMXRand_
#include "MMXRand.h"
#endif

#include "stdlib.h"
#include "time.h"



void MMXTwister::Seed(unsigned long seed)
	{
	if (seed == 0)
		{
#ifdef _DEBUG
		seed = 4357U;
#else
		seed = Hash(time(NULL) ^ clock());
#endif
		}

	unsigned long *s;
	int i;
	for( i = N, s = m_state; i-- ; s++)
		{
		*s    = seed & 0xffff0000;
		*s |= ( (seed *= 69069U)++ & 0xffff0000 ) >> 16;
		(seed *= 69069U)++;
		}
	Reload();
	}


//
// Implementation of MersenneTwister using the MMX instruction
// by Matthew Bellew (July 2000)
//
void MMXTwister::DoTheTwist()
	{
	DoTheTwistMMX();
	}


void MMXTwister::DoTheTwistC()
	{
	unsigned long *p = m_state;
	int i;
	for ( i = N - M; i--; )
		*p++ = twist(p[M], p[0], p[1]);
	for ( i = M; --i; )
		*p++ = twist( p[M-N], p[0], p[1] );

	*p = twist( p[M-N], p[0], m_state[0] );
	}


void MMXTwister::DoTheTwistMMX()
	{
	static __int64 doubleZero = 0x0000000000000000;
	static __int64 doubleOne  = 0x0000000100000001;
	static __int64 doubleMagic= 0x9908b0df9908b0df;
	static __int64 doubleLow  = 0x7fffffff7fffffff;
_asm
	{
	lea edi, DWORD PTR [ecx+4]	// skip vtable
	lea	esi, DWORD PTR [ecx+4]	// p = state

; for ( i = N - M; i--; )
;   *p++ = twist(p[M], p[0], p[1]);

	// NOTE
	// You can actually do three DWORDs at a time
	// by using the regular integer registers
	// However, it only helps performance about 10%

	mov	ecx, 113		; (624-397)/2=113 (+1)
$LOOP1:
	movq	mm1, QWORD PTR [esi]
	movq	mm2, QWORD PTR [esi+4]
	movq	mm0, mm1
	pxor	mm0, mm2
	pand	mm2, [doubleOne]
	pand	mm0, [doubleLow]
	pxor	mm0, mm1
	pcmpeqd	mm2, [doubleOne]
	psrld	mm0, 1
	pand	mm2, [doubleMagic]
	pxor	mm0, mm2
	pxor	mm0, QWORD PTR [esi+(397*4)]
	add		esi, 8
	dec		ecx
	movq	QWORD PTR [esi-8], mm0
	jne SHORT $LOOP1

	// left over
	mov	ebx, DWORD PTR [esi]
	mov	edx, DWORD PTR [esi+4]
	mov	eax, ebx
	xor	eax, edx
	and	dl, 1
	and	eax, 7ffffffeH
	xor	eax, ebx
	shr	eax, 1
	neg	dl
	sbb	edx, edx
	and	edx, 9908b0dfH
	xor	eax, edx
	xor	eax, DWORD PTR [esi+(397*4)]
	mov	DWORD PTR [esi], eax
	add	esi, 4

; for( i = M; --i; )
;   *p++ = twist( p[M-N], p[0], p[1] );

	mov	ecx, 198		; 397/2 = 198 (+1)
$LOOP2:
	movq	mm1, QWORD PTR [esi]
	movq	mm2, QWORD PTR [esi+4]
	movq	mm0, mm1
	pxor	mm0, mm2
	pand	mm2, [doubleOne]
	pand	mm0, [doubleLow]
	pxor	mm0, mm1
	pcmpeqd	mm2, [doubleOne]
	psrld	mm0, 1
	pand	mm2, [doubleMagic]
	pxor	mm0, mm2
	pxor	mm0, QWORD PTR [esi+(397-624)*4]
	add		esi, 8
	dec		ecx
	movq	QWORD PTR [esi-8], mm0
	jne SHORT $LOOP2

; *p = twist( p[M-N], p[0], state[0] );
	mov	ebx, DWORD PTR [esi]
	mov	edx, DWORD PTR [edi]
	mov	eax, ebx
	xor	eax, edx
	and	dl, 1
	and	eax, 7ffffffeH
	xor	eax, ebx
	shr	eax, 1
	neg	dl
	sbb	edx, edx
	and	edx, 9908b0dfH
	xor	eax, edx
	xor	eax, DWORD PTR [esi+(397-624)*4]
	mov	DWORD PTR [esi], eax

	emms
	}
	}

	
void MMXTwister::HashIt()
	{
	HashItMMX();
	}


void MMXTwister::HashItC()
	{
    for (int i=0 ; i<N ; i++)
		m_hashed[i] = Hash(m_state[i]);
	}


void MMXTwister::HashItASM()
	{
	_asm
	{
	lea esi, DWORD PTR [ecx+4]
	lea edi, DWORD PTR [ecx+624*4+4]
	mov	ecx, 623
$LOOP:
;		x=m_state[i]
	mov	eax, DWORD PTR [esi+ecx*4]
;		x ^= (x >> 11);
	mov	edx, eax
	shr	edx, 11
	xor	eax, edx
;		x ^= (x <<  7) & 0x9d2c5680U;
	mov	edx, eax
	shl	edx, 7
	and	edx, 9d2c5680H
	xor	eax, edx
;		x ^= (x << 15) & 0xefc60000U;
	mov	edx, eax
	shl	edx, 15
	and	edx, 0efc60000H
	xor	eax, edx
;		x = ( x ^ (x >> 18) );
	mov	edx, eax
	shr	edx, 18
	xor	eax, edx
;		m_hashed[i] = x
	mov	DWORD PTR [edi+ecx*4], eax
	sub	ecx, 1
	jns	SHORT $LOOP
	}
	}


void MMXTwister::HashItMMX()
	{
	static __int64 mask1 = 0x9d2c56809d2c5680;
	static __int64 mask2 = 0xefc60000efc60000;

	_asm
	{
	lea		esi, DWORD PTR [ecx+4]
	lea		edi, DWORD PTR [ecx+624*4+4]
	mov		ecx, 622
$LOOP:
;		x=m_state[i]
	movq	mm0, QWORD PTR [esi+ecx*4]
;		x ^= (x >> 11);
	movq	mm1, mm0
	psrld	mm1, 11		; shift right logical???
	pxor	mm0, mm1
;		x ^= (x <<  7) & 0x9d2c5680U;
	movq	mm1, mm0
	pslld	mm1, 7
	pand	mm1, QWORD PTR [mask1]
	pxor		mm0, mm1
;		x ^= (x << 15) & 0xefc60000U;
	movq	mm1, mm0
	pslld	mm1, 15
	pand	mm1, QWORD PTR [mask2]
	pxor		mm0, mm1
;		x = ( x ^ (x >> 18) );
	movq	mm1, mm0
	psrld	mm1, 18
	pxor	mm0, mm1
;		m_hashed[i] = x
	movq	QWORD PTR [edi+ecx*4], mm0
	sub		ecx, 2
	jns	SHORT $LOOP
	
	emms
	}
	}


//---------------------------------
// main.cpp
//


#include "windows.h"
#include "stdlib.h"
#include "stdio.h"


double elapsed()
	{
	static __int64 last = 0;
	static __int64 now;
	GetSystemTimeAsFileTime((FILETIME *)&now);
	__int64 e = now-last;
	last = now;
	return e / 10000000.0;
	}


class CTwister : public MMXTwister
	{
	protected:
		virtual void Reload()
			{
			DoTheTwistC();
			HashItC();
			m_next = N;
			}
	};


void main()
	{
	MMXTwister randMMX;
	CTwister   randC;
	randMMX.Seed(0x034fe934);
	randC.Seed(0x034fe934);

    for (long i=0; i<60; i++)
		{
		for (long j=0 ; j<1000*1000 ; j++)
			{
			if (randC() != randMMX())
				{
				fprintf(stderr, "error");
				exit(1);
				}
			}
		}


	elapsed();
    for (long i=0; i<60; i++)
		{
		for (long j=0 ; j<1000*1000 ; j++)
			{
			DWORD x = randMMX();
			// printf("%08x\n",x);
			}
		}
	printf("MMXTwister %3.4lf\n", elapsed());
	}