// $Id: RandEngine.cc,v 1.28 2004/08/27 19:08:48 marafino Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
//                             HEP Random
//                         --- RandEngine ---
//                      class implementation file
// -----------------------------------------------------------------------
// This file is part of Geant4 (simulation toolkit for HEP).

// =======================================================================
// Gabriele Cosmo - Created: 5th September 1995
//                - Minor corrections: 31st October 1996
//                - Added methods for engine status: 19th November 1996
//                - mx is initialised to RAND_MAX: 2nd April 1997
//                - Fixed bug in setSeeds(): 15th September 1997
//                - Private copy constructor and operator=: 26th Feb 1998
// J.Marraffino   - Added stream operators and related constructor.
//                  Added automatic seed selection from seed table and
//                  engine counter. Removed RAND_MAX and replaced by
//                  pow(0.5,32.). Flat() returns now 32 bits values
//                  obtained by concatenation: 15th Feb 1998
// Ken Smith      - Added conversion operators:  6th Aug 1998
// J. Marraffino  - Remove dependence on hepString class  13 May 1999
// M. Fischler    - Rapaired bug that in flat() that relied on rand() to
//                  deliver 15-bit results.  This bug was causing flat()
//		    on Linux systems to yield randoms with mean of 5/8(!)
//		  - Modified algorithm such that on systems with 31-bit rand()
//		    it will call rand() only once instead of twice. Feb 2004
// M. Fischler    - Modified the general-case template for RandEngineBuilder
//		    such that when RAND_MAX is an unexpected value the routine
//		    will still deliver a sensible flat() random. Feb 2004
// M. Fischler	  - Tweaked a line of code intended to test whehter RANDMAX
//                  delivers 32-bits or more.  The original code emitted a
//                  warning, as the test makes little sense working with 
//                  32-bit integers.  Aug 2004
//
// =======================================================================

#include "CLHEP/Random/RandEngine.h"
#include "CLHEP/Random/Random.h"
#include <string.h>
#include <math.h>

#ifdef NOTDEF
#undef  RAND_MAX
#define RAND_MAX 9382956
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/RandFlat.h"
MTwistEngine * fakeFlat = new MTwistEngine;
RandFlat rflat (fakeFlat, 0, RAND_MAX+1);
int rand() { return (int)rflat(); }
#endif

static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 

// number of instances with automatic seed selection
int RandEngine::numEngines = 0;

// Maximum index into the seed table
int RandEngine::maxIndex = 215;

RandEngine::RandEngine(long seed) 
: mantissa_bit_32( pow(0.5,32.) )
{
   setSeed(seed,0);
   setSeeds(&theSeed,0);
   seq = 0;
}

RandEngine::RandEngine(int rowIndex, int colIndex)
: mantissa_bit_32( pow(0.5,32.) )
{
   long seeds[2];
   long seed;

   int cycle = abs(int(rowIndex/maxIndex));
   int row = abs(int(rowIndex%maxIndex));
   int col = abs(int(colIndex%2));
   long mask = ((cycle & 0x000007ff) << 20 );
   HepRandom::getTheTableSeeds( seeds, row );
   seed = (seeds[col])^mask;
   setSeed(seed,0);
   setSeeds(&theSeed,0);
   seq = 0;
}

RandEngine::RandEngine()
: mantissa_bit_32( pow(0.5,32.) )
{
   long seeds[2];
   long seed;
   int cycle,curIndex;

   cycle = abs(int(numEngines/maxIndex));
   curIndex = abs(int(numEngines%maxIndex));
   numEngines += 1;
   long mask = ((cycle & 0x007fffff) << 8);
   HepRandom::getTheTableSeeds( seeds, curIndex );
   seed = seeds[0]^mask;
   setSeed(seed,0);
   setSeeds(&theSeed,0);
   seq = 0;
}

RandEngine::RandEngine(HepStd::istream& is)
: mantissa_bit_32( pow(0.5,32.) )
{
   is >> *this;
}

RandEngine::~RandEngine() {}

RandEngine::RandEngine(const RandEngine &p)
: mantissa_bit_32( pow(0.5,32.) )
{
  // Assignment and copy of RandEngine objects may provoke
  // undesired behavior in a single thread environment.
  
  HepStd::cerr << "!!! WARNING !!! - Illegal operation." << HepStd::endl;
  HepStd::cerr << "- Copy constructor and operator= are NOT allowed on "
	    << "RandEngine objects -" << HepStd::endl;
  *this = p;
}

RandEngine & RandEngine::operator = (const RandEngine &p)
{
  // Assignment and copy of RandEngine objects may provoke
  // undesired behavior in a single thread environment.

  HepStd::cerr << "!!! WARNING !!! - Illegal operation." << HepStd::endl;
  HepStd::cerr << "- Copy constructor and operator= are NOT allowed on "
	    << "RandEngine objects -" << HepStd::endl;
  *this = p;
  return *this;
}

void RandEngine::setSeed(long seed, int)
{
   theSeed = seed;
   srand( int(seed) );
   seq = 0;
}

void RandEngine::setSeeds(const long* seeds, int)
{
  setSeed(seeds ? *seeds : 19780503L, 0);
  theSeeds = seeds;
}

void RandEngine::saveStatus( const char filename[] ) const
{
   HepStd::ofstream outFile( filename, HepStd::ios::out ) ;

   if (!outFile.bad()) {
     outFile << theSeed << HepStd::endl;
     outFile << seq << HepStd::endl;
   }
}

void RandEngine::restoreStatus( const char filename[] )
{
   // The only way to restore the status of RandEngine is to
   // keep track of the number of shooted random sequences, reset
   // the engine and re-shoot them again. The Rand algorithm does
   // not provide any way of getting its internal status.

   HepStd::ifstream inFile( filename, HepStd::ios::in);
   long count;

   if (!inFile.bad() && !inFile.eof()) {
     inFile >> theSeed;
     inFile >> count;
     setSeed(theSeed,0);
     for (int i=0; i<count; ++i) flat();
   }
}

void RandEngine::showStatus() const
{
   HepStd::cout << HepStd::endl;
   HepStd::cout << "---------- Rand engine status ----------" << HepStd::endl;
   HepStd::cout << " Initial seed  = " << theSeed << HepStd::endl;
   HepStd::cout << " Shooted sequences = " << seq << HepStd::endl;
   HepStd::cout << "----------------------------------------" << HepStd::endl;
}

// Here we set up such that **at compile time**, the compiler decides based on
// RAND_MAX how to form a random double with 32 leading random bits, using
// one or two calls to rand().  Some of the lowest order bits of 32 are allowed to
// be as weak as mere XORs of some higher bits, but not to be always fixed.
//
// The method decision is made at compile time, rather than using a run-time
// if on the value of RAND_MAX.  Although it is possible to cope with arbitrary
// values of RAND_MAX of the form 2**N-1, with the same efficiency, the 
// template techniques needed would look mysterious and perhaps over-stress 
// some older compilers.  We therefore only treat RAND_MAX = 2**15-1 (as on 
// most older systems) and 2**31-1 (as on the later Linux systems) as special
// and super-efficient cases.  We detect any different value, and use an
// algorithm which is correct even if RAND_MAX is not one less than a power 
// of 2.

template <int> struct RandEngineBuilder {     // RAND_MAX any arbitrary value
  static unsigned int thirtyTwoRandomBits() {
  
  static bool prepared = false;
  static unsigned int iT;
  static unsigned int iK;
  static unsigned int iS;
  static int iN;
  static double fS;
  static double fT;
  
  if ( (RAND_MAX >> 31) > 0 ) 
  {
    // Here, we know that integer arithmetic is 64 bits.
    if ( !prepared ) {
      iS = RAND_MAX + 1;
      iK = 1;
//    int StoK = S;
      int StoK = iS;
      int HALF_RAND_MAX = RAND_MAX >> 1; // Tweak to avoid warning on 32-bit 
      if ( (HALF_RAND_MAX >> 31) == 0) { // shift:  Use HALF_RAND_MAX
        iK = 2;
//	StoK = S*S;
	StoK = iS*iS;
      }
      int m;
      for ( m = 0; m < 64; ++m ) {
        StoK >>= 1;
	if (StoK == 0) break; 
      }   
      iT = 1 << m;
      prepared = true;
    }
    int v = 0;
    do { 
      for ( int i = 0; i < iK; ++i ) {
        v = iS*v+rand();
      }
    } while (v < iT);   
    return v & 0xFFFFFFFF;
    
  } 
  else if ( (RAND_MAX >> 26) == 0 ) 
  {
    // Here, we know it is safe to work in terms of doubles without loss 
    // of precision, but we have no idea how many randoms we will need to
    // generate 32 bits.
    if ( !prepared ) {
      fS = RAND_MAX + 1;
      double twoTo32 = ldexp(1.0,32);
      double StoK = fS;
      for ( iK = 1; StoK < twoTo32; StoK *= fS, iK++ ) { } 
      int m;
      fT = 1.0;
      for ( m = 0; m < 64; ++m ) {
        StoK *= .5;
	if (StoK < 1.0) break; 
	fT *= 2.0;
      }   
      prepared = true;
    }
    double v = 0;
    do { 
      for ( int i = 0; i < iK; ++i ) {
        v = fS*v+rand();
      }
    } while (v < fT);    
    return ((unsigned int)v) & 0xFFFFFFFF;
    
  } 
  else 
  {
    // Here, we know that 16 random bits are available from each of 
    // two random numbers.
    if ( !prepared ) {
      iS = RAND_MAX + 1;
      int SshiftN = iS;
      for (iN = 0; SshiftN > 1; SshiftN >>= 1, iN++) { }
      iN -= 17;
    prepared = true;
    }
    unsigned int x1, x2;
    do { x1 = rand(); } while (x1 < (1<<16) );
    do { x2 = rand(); } while (x2 < (1<<16) );
    return x1 | (x2 << 16);    
  }
  
  }
};

template <> struct RandEngineBuilder<2147483647> { // RAND_MAX = 2**31 - 1
  inline static unsigned int thirtyTwoRandomBits() {
    unsigned int x = rand() << 1;	// bits 31-1
    x ^= ( (x>>23) ^ (x>>7) ) ^1;	// bit 0 (weakly pseudo-random)
    return x & 0xFFFFFFFF;		// mask in case int is 64 bits 
    }
};

template <> struct RandEngineBuilder<32767> { // RAND_MAX = 2**15 - 1
  inline static unsigned int thirtyTwoRandomBits() {
    unsigned int x = rand() << 17;	// bits 31-17
    x ^= rand() << 2;			// bits 16-2
    x ^= ( (x>>23) ^ (x>>7) ) ^3;	// bits  1-0 (weakly pseudo-random)
    return x & 0xFFFFFFFF;		// mask in case int is 64 bits 
    }
};

double RandEngine::flat()
{
  double r;
  do { r = RandEngineBuilder<RAND_MAX>::thirtyTwoRandomBits(); 
     } while ( r == 0 );
  return r/4294967296.0;
}

void RandEngine::flatArray(const int size, double* vect)
{
   int i;

   for (i=0; i<size; ++i)
     vect[i]=flat();
}

RandEngine::operator unsigned int() {
  return RandEngineBuilder<RAND_MAX>::thirtyTwoRandomBits();
}

HepStd::ostream & operator << ( HepStd::ostream& os, const RandEngine& e ) 
{
     char beginMarker[] = "RandEngine-begin";
     char endMarker[]   = "RandEngine-end";

     os << " " << beginMarker << " ";
     os << e.theSeed << " " << e.seq << " ";
     os << endMarker << " ";
     return os;
}

HepStd::istream & operator >> ( HepStd::istream& is, RandEngine& e )
{
   // The only way to restore the status of RandEngine is to
   // keep track of the number of shooted random sequences, reset
   // the engine and re-shoot them again. The Rand algorithm does
   // not provide any way of getting its internal status.

  long count;
  char beginMarker [MarkerLen];
  char endMarker   [MarkerLen];

  is >> HepStd::ws;
  is.width(MarkerLen);  // causes the next read to the char* to be <=
			// that many bytes, INCLUDING A TERMINATION \0 
			// (Stroustrup, section 21.3.2)
  is >> beginMarker;
  if (strcmp(beginMarker,"RandEngine-begin")) {
     is.clear(HepStd::ios::badbit | is.rdstate());
     HepStd::cerr << "\nInput stream mispositioned or"
	       << "\nRandEngine state description missing or"
	       << "\nwrong engine type found." << HepStd::endl;
     return is;
  }
  is >> e.theSeed;
  is >> count;
  is >> HepStd::ws;
  is.width(MarkerLen);  
  is >> endMarker;
  if (strcmp(endMarker,"RandEngine-end")) {
     is.clear(HepStd::ios::badbit | is.rdstate());
     HepStd::cerr << "\nRandEngine state description incomplete."
	       << "\nInput stream is probably mispositioned now." << HepStd::endl;
     return is;
   }

   e.setSeed(e.theSeed,0);
   for (int i=0; i<count; ++i)
      { e.flat(); }	// { dummy = flat(); }
   return is;
}