/*<<<***************************************************************************
The MIT License (MIT)

Copyright(C) 2013 - 2015 Marc Mirti Mancinelli, Christophe Lardau and Philippe Pamart (xdcube.org)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
////////////////////////////////////////////////////////////////////////////////
#Title: xbRandom
#Project: xdcBase
////////////////////////////////////////////////////////////////////////////////
#Author: Marc MIRTI MANCINELLI
#History:
<02/06/2013: creation.
12/02/2014: this notice inserted.
***************************************************************************>>>*/
#include "stdafx.h"

#include "xbRandom.h"

xbRanrotWGenerator g_xbRand;

/************************* RANROTW.CPP ****************** AgF 1999-03-03 *
*  Random Number generator 'RANROT' type W                               *
*  This version is used when a resolution higher that 32 bits is desired.*
*                                                                        *
*  This is a lagged-Fibonacci type of random number generator with       *
*  rotation of bits.  The algorithm is:                                  *
*  Z[n] = (Y[n-j] + (Y[n-k] rotl r1)) modulo 2^(b/2)                     *
*  Y[n] = (Z[n-j] + (Z[n-k] rotl r2)) modulo 2^(b/2)                     *
*  X[n] = Y[n] + Z[n]*2^(b/2)                                            *
*                                                                        *
*  The last k values of Y and Z are stored in a circular buffer named    *
*  randbuffer.                                                           *
*  The code includes a self-test facility which will detect any          *
*  repetition of previous states.                                        *
*  The function uses a fast method for conversion to floating point.     *
*  This method relies on floating point numbers being stored in the      *
*  standard 64-bit IEEE format or the 80-bit double format.         *
*                                                                        *
*  The theory of the RANROT type of generators and the reason for the    *
*  self-test are described at www.agner.org/random/theory                *
*                                                                        *
* ©2002 A. Fog. GNU General Public License www.gnu.org/copyleft/gpl.html *
*************************************************************************/

// If your system doesn't have a rotate function for 32 bits integers,
// then use the definition below. If your system has the _lrotl function 
// then remove this.
// ulong _lrotl (ulong x, long r) {
//   return (x << r) | (x >> (sizeof(x)*8-r));}


ulong xbRanrotWGenerator::BRandom() 
{
	// generate next random number
	ulong y, z;
	// generate next number
	z = _lrotl(randbuffer[p1][0], R1) + randbuffer[p2][0];
	y = _lrotl(randbuffer[p1][1], R2) + randbuffer[p2][1];
	randbuffer[p1][0] = y; randbuffer[p1][1] = z;
	// rotate list pointers
	if (--p1 < 0) p1 = KK - 1;
	if (--p2 < 0) p2 = KK - 1;
	// perform self-test
	if (randbuffer[p1][0] == randbufcopy[0][0] &&
		memcmp(randbuffer, randbufcopy[KK-p1], 2*KK*sizeof(ulong)) == 0) 
	{
		// self-test failed
		if ((p2 + KK - p1) % KK != JJ) 
			// note: the way of printing error messages depends on system
			// In Windows you may use FatalAppExit
			xbThrowEsys("Random number generator not initialized");
		/*else 
			xbThrowEsys("Random number generator returned to initial state");*/
	}
	
	randbits[0] = y;
	randbits[1] = z;
	
	return y;
}


double xbRanrotWGenerator::RandomDouble() 
{
	if (Architecture == archNotInit)
	{
		SYSTEMTIME systi;
		::GetSystemTime(&systi);
		ulong rinit = ulong(uint64(systi.wMilliseconds )+ uint64(systi.wHour)*3600*1000+uint64(systi.wMinute)*60*1000+uint64(systi.wSecond)*1000
			+::GetTickCount());
		RandomInit(rinit);
	}

	// returns a random number between 0 and 1.
	ulong z = BRandom();  // generate 64 random bits

	switch (Architecture) 
	{
	case archEXTENDEDPRECISIONLITTLEENDIAN:
		// 80 bits floats = 63 bits resolution
		randbits[1] = z | 0x80000000;  break;
	case archLITTLEENDIAN:
		// 64 bits floats = 52 bits resolution
		randbits[1] = (z & 0x000FFFFF) | 0x3FF00000;  break;
	case archBIGENDIAN:
		// 64 bits floats = 52 bits resolution
		randbits[0] = (randbits[0] & 0x000FFFFF) | 0x3FF00000;  break;
	case archNONIEEE: default:
		// not a recognized floating point format. 32 bits resolution
		return (double)z * (1./((double)(ulong)(-1L)+1.));
	}
	return randp1 - 1.0;
}


long xbRanrotWGenerator::RandomInt32(long min, long max) 
{
	// get integer random number in desired interval
	long iinterval = max - min + 1;
//	if (iinterval <= 0) return 0x80000000;  // error
	long i = long(iinterval * RandomDouble());      // truncate
	if (i >= iinterval) i = iinterval-1;
	else if(i<0)
		i=0;
	return min + i;
}

int64 xbRanrotWGenerator::RandomInt64(int64 min, int64 max) 
{
	// get integer random number in desired interval
	int64 iinterval = max - min + 1;
//	if (iinterval <= 0) return 0x8000000000000000;  // error
	int64 i = int64(iinterval * RandomDouble());      // truncate
	if (i >= iinterval) i = iinterval-1;
	else if(i<0)
		i=0;
	return min + i;
}

double xbRanrotWGenerator::RandomDouble(double min, double max) 
{
	// get integer random number in desired interval
	double iinterval = max - min + 1;
//	if (iinterval <= 0) return 0x8000000000000000;  // error
	double i = double(iinterval * RandomDouble());      // truncate
	if (i >= iinterval) i = iinterval-1;
	else if(i<0)
		i=0;
	return min + i;
}

void xbRanrotWGenerator::RandomInit (ulong seed) 
{
	// this function initializes the random number generator.
	long i, j;
	
	// make random numbers and put them into the buffer
	for (i=0; i<KK; i++) 
	{
		for (j=0; j<2; j++) 
		{
			seed = seed * 2891336453UL + 1;
			randbuffer[i][j] = seed;
		}
	}
	
	// set exponent of randp1
	randbits[2] = 0; randp1 = 1.0;
	// initialize pointers to circular buffer
	p1 = 0;  p2 = JJ;
	// store state for self-test
	memcpy (randbufcopy, randbuffer, 2*KK*sizeof(ulong));
	memcpy (randbufcopy[KK], randbuffer, 2*KK*sizeof(ulong));
	// randomize some more
	
	for (i=0; i<31; i++) BRandom();

	if (Architecture == archNotInit)
	{
		// detect computer architecture
		randbits[2] = 0; randp1 = 1.0;
		if (randbits[2] == 0x3FFF) Architecture = archEXTENDEDPRECISIONLITTLEENDIAN;
		else if (randbits[1] == 0x3FF00000) Architecture = archLITTLEENDIAN;
		else if (randbits[0] == 0x3FF00000) Architecture = archBIGENDIAN;
		else Architecture = archNONIEEE;
	}
}

