/*------------------------------------------------------------*/ 

/************  Random numbers ***********/ 

/* Tausworthes generator. This file contains C-version. */ 

/* S. Tezuka and P.L Equyer: Efficient portable combined Tausworthe 
 * random number generators. ACM. Transactions on Modelling and 
 * Computer Simulation 1 (1991) 99-112.
 * Time in HP 720 (Tarzan): 1.36 microsec. 
 *
 * Contents:
 * Long random integers in [0,2^31-1]:	BRAND() (Macro)
 * Long random integers in [ia,ib] (long): irand(ia,ib) 
 * Float random numbers in (0,1):	frand()
 * Double random numbers in (0,1):	drand()
 * Set seeds to ia and ib (unsigned long): setrangen2(ia,ib) 
 * Store seeds to ia and ib (unsigned long): getseeds2(&ia,&ib) */

/*------------------------------------------------------------*/

#define BRAND CombTaus /* Basic random numbers in [1..2**31-1] */

/* Random integers in interval [a,b) (b-a < 100000) */
#define IRI(ia,ib) ((ia) + CombTaus() % ((ib) - (ia)))

/* Random integers in interval [0,z) (z < 100000) */
#define IRZ(z) (CombTaus() % (z))

extern long CombTaus(void);
extern long irand(long,long);
extern double drand(void);
extern float frand(void);
extern void setrangen2(unsigned long,unsigned long);
extern void getseeds2(unsigned long*,unsigned long*);

/*------------------------------------------------------------*/

#include <cstdio>
#include <cstdlib>

unsigned long I1 = 12345UL, I2 = 67890UL; /* Seeds. Global */ 

/* Definitions of constants */
#define Mask1 0x7fffffffUL
#define Mask2 0x1fffffffUL
#define Todoub 4.656612873e-10
#define Tofloat 4.656612e-10f

/* Basic random sequence. C-version */
long CombTaus(void)
{
unsigned long b;
b = ((I1 << 13) ^ I1) & Mask1;
I1 = ((I1 << 12) ^ (b >> 19)) & Mask1;
b = ((I2 << 2) ^ I2) & Mask2;
I2 = ((I2 << 17) ^ (b >> 12)) & Mask2;
return I1 ^ (I2 << 2);
}

/*+++ Random integer in ia..ib +++*/
long irand (long ia, long ib)
{
long range = ib-ia+1,rn;
do {rn = BRAND()/(Mask1/range); } while (rn >= range); return ia + rn;
}

/*+++ Random double in [0,1) +++*/
double drand(void)
{ return Todoub*BRAND(); }

/*+++ Random float in [0,1) +++*/
float frand(void)
{ return Tofloat*BRAND(); }

/*+++ Set seeds of the random number generator +++*/ 
void setrangen2(unsigned long seed1, unsigned long seed2) 
{ I1 = seed1 & Mask1; I2 = seed2 & Mask2; } 

/*+++ Get current seeds from the random number generator +++*/ 
void getseeds2(unsigned long *pia, unsigned long *pib) 
{ *pia = I1; *pib = I2; }

/*------------------------------------------------------------*/ 

/************  Random permutations ***********/ 

/* Swap the contents of a[i] and a[j] */
void swap(unsigned int* a, unsigned int i, unsigned j) {
  unsigned int t;

  t = a[i]; a[i] = a[j]; a[j] = t;
}

/* Permute the n numbers in array a */
void permute( unsigned int* a, unsigned int n) {
   unsigned int i;
   unsigned int j;

   setrangen2(4333, 78901234);
   for (i = n - 1; i > 0; i--) {
     j = irand(0, i);
     swap(a, i, j);
   }
}