/*	"cswrand.c"
 *
 *	entries:
 *		real random(seed)  - real uniform in 0 .. 1
 *	 	int irandm(seed)   - integer uniform in 0 .. 2**31-1
 *		void initseed(see) - use system clock to initialize seed vector
 *		int seed[] ;  integer seed vector of 2 elements
 */

/*
 *	A random number generator called as a function by
 *	random (iseed)	or	irandm (iseed)
 *	The parameter is called by reference, not by value.
 *	The routine(s) are callable directly as Fortran functions.
 *	From C, I expect you have to pass a pointer to a 2-element vector
 *	The first call gives a real uniform in 0 .. 1.
 *	The second gives an integer uniform in 0 .. 2**31-1
 *	Both update iseed in exactly the same way.
 *	iseed must be a 2-element integer vector.
 *	The initial value of the second element may be anything.
 *
 *	The period of the random sequence is 2**32 * (2**32-1)
 *	The table mt(0:127) is defined by mt(i) = 69069 ** (128-i)
 */
#include "rand.h"
int seed[2] ;

/*
 *  use system clock to initialize the seed vector for random number generator
 */
void GenSeed(int seed[])
{
	long clock ;

	clock = time((long *) 0) ;
	seed[0] = (int) clock ;
	clock = time((long *) 0) ;
	seed[1] = (int) clock ;
}


static int mt[] = {
	      902906369,
	     2030498053,
	     -473499623,
	     1640834941,
	      723406961,
	     1993558325,
	     -257162999,
	    -1627724755,
	      913952737,
	      278845029,
	     1327502073,
	    -1261253155,
	      981676113,
	    -1785280363,
	     1700077033,
	      366908557,
	    -1514479167,
	     -682799163,
	      141955545,
	     -830150595,
	      317871153,
	     1542036469,
	     -946413879,
	    -1950779155,
	      985397153,
	      626515237,
	      530871481,
	      783087261,
	    -1512358895,
	     1031357269,
	    -2007710807,
	    -1652747955,
	    -1867214463,
	      928251525,
	     1243003801,
	    -2132510467,
	     1874683889,
	     -717013323,
	      218254473,
	    -1628774995,
	    -2064896159,
	       69678053,
	      281568889,
	    -2104168611,
	     -165128239,
	     1536495125,
	      -39650967,
	      546594317,
	     -725987007,
	     1392966981,
	     1044706649,
	      687331773,
	    -2051306575,
	     1544302965,
	     -758494647,
	    -1243934099,
	      -75073759,
	      293132965,
	    -1935153095,
	      118929437,
	      807830417,
	    -1416222507,
	    -1550074071,
	      -84903219,
	     1355292929,
	     -380482555,
	    -1818444007,
	     -204797315,
	      170442609,
	    -1636797387,
	      868931593,
	     -623503571,
	     1711722209,
	      381210981,
	     -161547783,
	     -272740131,
	    -1450066095,
	     2116588437,
	     1100682473,
	      358442893,
	    -1529216831,
	     2116152005,
	     -776333095,
	     1265240893,
	     -482278607,
	     1067190005,
	      333444553,
	       86502381,
	      753481377,
	       39000101,
	     1779014585,
	      219658653,
	     -920253679,
	     2029538901,
	     1207761577,
	    -1515772851,
	     -236195711,
	      442620293,
	      423166617,
	    -1763648515,
	     -398436623,
	    -1749358155,
	     -538598519,
	     -652439379,
	      430550625,
	    -1481396507,
	     2093206905,
	    -1934691747,
	     -962631983,
	     1454463253,
	    -1877118871,
	     -291917555,
	    -1711673279,
	      201201733,
	     -474645415,
	      -96764739,
	    -1587365199,
	     1945705589,
	     1303896393,
	     1744831853,
	      381957665,
	     2135332261,
	      -55996615,
	    -1190135011,
	     1790562961,
	    -1493191723,
	      475559465,
	          69069 
		} ;
#define mask  593970775     			/* = hex 23674657 */
#define scale 4.6566128730773926e-010 	/* = 2.0 ** -31 */

float random (int is[])
{ 
	int it, nit, leh ;
	float retval ;

	it = is[0] ;
	leh = is[1] ;
	if (it <= 0) 
		it = (it + it) ^ mask ;
	else
		it = it + it ;

	nit = it - 1 ;
	/* 
	 *	to ensure all-ones pattern omitted
	 */
	leh = leh * mt[nit & 127] + nit ;
	is[0] = it ;
	is[1] = leh ;
	if (leh < 0) leh = ~ leh ;
	
	retval = scale * (float) (leh | 1) ;
	return retval ;
}
 

static int imt[] = {
	      902906369,
	     2030498053,
	     -473499623,
	     1640834941,
	      723406961,
	     1993558325,
	     -257162999,
	    -1627724755,
	      913952737,
	      278845029,
	     1327502073,
	    -1261253155,
	      981676113,
	    -1785280363,
	     1700077033,
	      366908557,
	    -1514479167,
	     -682799163,
	      141955545,
	     -830150595,
	      317871153,
	     1542036469,
	     -946413879,
	    -1950779155,
	      985397153,
	      626515237,
	      530871481,
	      783087261,
	    -1512358895,
	     1031357269,
	    -2007710807,
	    -1652747955,
	    -1867214463,
	      928251525,
	     1243003801,
	    -2132510467,
	     1874683889,
	     -717013323,
	      218254473,
	    -1628774995,
	    -2064896159,
	       69678053,
	      281568889,
	    -2104168611,
	     -165128239,
	     1536495125,
	      -39650967,
	      546594317,
	     -725987007,
	     1392966981,
	     1044706649,
	      687331773,
	    -2051306575,
	     1544302965,
	     -758494647,
	    -1243934099,
	      -75073759,
	      293132965,
	    -1935153095,
	      118929437,
	      807830417,
	    -1416222507,
	    -1550074071,
	      -84903219,
	     1355292929,
	     -380482555,
	    -1818444007,
	     -204797315,
	      170442609,
	    -1636797387,
	      868931593,
	     -623503571,
	     1711722209,
	      381210981,
	     -161547783,
	     -272740131,
	    -1450066095,
	     2116588437,
	     1100682473,
	      358442893,
	    -1529216831,
	     2116152005,
	     -776333095,
	     1265240893,
	     -482278607,
	     1067190005,
	      333444553,
	       86502381,
	      753481377,
	       39000101,
	     1779014585,
	      219658653,
	     -920253679,
	     2029538901,
	     1207761577,
	    -1515772851,
	     -236195711,
	      442620293,
	      423166617,
	    -1763648515,
	     -398436623,
	    -1749358155,
	     -538598519,
	     -652439379,
	      430550625,
	    -1481396507,
	     2093206905,
	    -1934691747,
	     -962631983,
	     1454463253,
	    -1877118871,
	     -291917555,
	    -1711673279,
	      201201733,
	     -474645415,
	      -96764739,
	    -1587365199,
	     1945705589,
	     1303896393,
	     1744831853,
	      381957665,
	     2135332261,
	      -55996615,
	    -1190135011,
	     1790562961,
	    -1493191723,
	      475559465,
	          69069
		} ;

/*
 *	A random number generator called as a function by
 *	random (iseed)	or	irandm (iseed)
 *	The parameter is called by reference, not by value.
 *	The routine(s) are callable directly as Fortran functions.
 *	From C, I expect you have to pass a pointer to a 2-element vector
 *	The first call gives a real uniform in 0 .. 1.
 *	The second gives an integer uniform in 0 .. 2**31-1
 *	Both update iseed in exactly the same way.
 *	iseed must be a 2-element integer vector.
 *	The initial value of the second element may be anything.
 *
 *	The period of the random sequence is 2**32 * (2**32-1)
 *	The table mt(0:127) is defined by mt(i) = 69069 ** (128-i)
 */

/* #define  mask  593970775		*/
/* parameter (scale = 2.0 ** -31)*/

int irandm(int is[])
{
	int itemp, nit, it, leh ;

	it = is[0] ;
	leh = is[1] ;
	if (it <= 0) 
		it = (it + it) ^ mask ;
	else
		it = it + it ;
	nit = it - 1 ;
	/* leh = leh * imt[nit & 127] + nit ;*/
	itemp = nit & 127 ;
	leh = leh * imt[itemp] + nit ;
	is[0] = it ;
	is[1] = leh ;
	/*
 	 *	Want a positive result but dont want to lose top bit.
 	 */
	if (leh < 0) leh = ~ leh ;
	return leh ;
}