// $Id: TripleRand.cc,v 1.18 2002/05/07 16:47:30 mf Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
//                           Hep Random
//                       --- TripleRand ---
//                   class implementation file
// -----------------------------------------------------------------------
// A canopy pseudo-random number generator. Using the Tausworthe
// exclusive-or shift register, a simple Integer Coungruence generator, and
// the Hurd 288 total bit shift register, all XOR'd with each other, we
// provide an engine that should be a fairly good "mother" generator.
// =======================================================================
// Ken Smith      - Initial draft started:                    23rd Jul 1998
//                - Added conversion operators:                6th Aug 1998
// J. Marraffino  - Added some explicit casts to deal with
//                  machines where sizeof(int) != sizeof(long)  22 Aug 1998
// M. Fischler    - Modified use of the various exponents of 2
//                  to avoid per-instance space overhead and
//                  correct the rounding procedure              15 Sep 1998
//                - modify constructors so that no sequence can 
//		    ever accidentally be produced by differnt 
//		    seeds, and so that exxcept for the default
//		    constructor, the seed fully determines the 
//		    sequence.              			15 Sep 1998
// M. Fischler	  - Eliminated dependency on hepString		13 May 1999
// M. Fischler	  - Eliminated Taus() and Cong() which were	
//		    causing spurions warnings on SUN CC		27 May 1999
// M. Fischler    - Put endl at end of puts of Tausworth and	10 Apr 2001
//		    integerCong
//
// =======================================================================

#include "CLHEP/Random/TripleRand.h"
#include <string.h>

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

//********************************************************************
//   TripleRand  
//********************************************************************

// Number of instances with automatic seed selection
int TripleRand::numEngines = 0;

double TripleRand::twoToMinus_32;
double TripleRand::twoToMinus_53;
double TripleRand::nearlyTwoToMinus_54;

void TripleRand::powersOfTwo() {
  twoToMinus_32 = ldexp (1.0, -32);
  twoToMinus_53 = ldexp (1.0, -53);
  nearlyTwoToMinus_54 = ldexp (1.0, -54) - ldexp (1.0, -100);
}

TripleRand::TripleRand() 
: tausworthe (1234567 + numEngines + 175321),
  integerCong(69607 * tausworthe + 54329, numEngines),
  hurd(19781127 + integerCong)
{  
  powersOfTwo();
  theSeed = 1234567;
  ++numEngines;
}

TripleRand::TripleRand(long seed)
: tausworthe ((unsigned int)seed + 175321),
  integerCong(69607 * tausworthe + 54329, 1313),
  hurd(19781127 + integerCong)
{
  powersOfTwo();
  theSeed = seed;
}

TripleRand::TripleRand(HepStd::istream & is) 
{
  is >> *this;
}

TripleRand::TripleRand(int rowIndex, int colIndex)
: tausworthe (rowIndex + numEngines * colIndex + 175321),
  integerCong(69607 * tausworthe + 54329, 19),
  hurd(19781127 + integerCong)
{
  powersOfTwo();
  theSeed = rowIndex;
}

TripleRand::~TripleRand() { }

TripleRand::TripleRand( const TripleRand & p )
{
  *this = p;
}

TripleRand & TripleRand::operator=( const TripleRand & p ) {
  if( this != &p ) {
     tausworthe  = p.tausworthe;
     integerCong = p.integerCong;
     hurd        = p.hurd;
  }
  return *this;
}

double TripleRand::flat() {
  unsigned int ic ( integerCong );
  unsigned int t  ( tausworthe  );
  unsigned int h  ( hurd        );
  return ( (t  ^ ic ^ h) * twoToMinus_32 +          // most significant part
               (h >> 11) * twoToMinus_53 +          // fill in remaining bits
                     nearlyTwoToMinus_54	    // make sure non-zero
         );
}

void TripleRand::flatArray(const int size, double* vect) {
  for (int i = 0; i < size; ++i) {
    vect[i] = flat();
  }
}

void TripleRand::setSeed(long seed, int) {
  theSeed = seed;
  tausworthe  = Tausworthe((unsigned int)seed + numEngines + 175321);
  integerCong = IntegerCong(69607 * tausworthe + 54329, numEngines);
  hurd        = Hurd288Engine( 19781127 + integerCong );
}

void TripleRand::setSeeds(const long * seeds, int) {
  setSeed(seeds ? *seeds : 1234567, 0);
  theSeeds = seeds;
}

void TripleRand::saveStatus(const char filename[]) const {
  HepStd::ofstream outFile(filename, HepStd::ios::out);
  if (!outFile.bad()) {
    outFile << HepStd::setprecision(20) << theSeed << " ";
    tausworthe.put ( outFile );
    integerCong.put( outFile);
    outFile << ConstHurd() << HepStd::endl;
  }
}

void TripleRand::restoreStatus(const char filename[]) {
  HepStd::ifstream inFile(filename, HepStd::ios::in);
  if (!inFile.bad()) {
    inFile >> theSeed;
    tausworthe.get ( inFile );
    integerCong.get( inFile );
    inFile >> Hurd();
  }
}

void TripleRand::showStatus() const {
  HepStd::cout << HepStd::setprecision(20) << HepStd::endl;
  HepStd::cout <<         "-------- TripleRand engine status ---------" 
								<< HepStd::endl;
  HepStd::cout << "Initial seed           = " << theSeed << HepStd::endl;
  HepStd::cout << "Tausworthe generator   = " << HepStd::endl;
  tausworthe.put( HepStd::cout ); 
  HepStd::cout << "IntegerCong generator  = " << HepStd::endl;
  integerCong.put( HepStd::cout );
  HepStd::cout << "Hurd288Engine generator= " << HepStd::endl << ConstHurd();
  HepStd::cout << HepStd::endl << "-----------------------------------------" 
								<< HepStd::endl;
}

TripleRand::operator float() {
  return (float)
    ( ( integerCong ^ tausworthe ^ hurd ) * twoToMinus_32 
					+ nearlyTwoToMinus_54 );
					// make sure non-zero!
}

TripleRand::operator unsigned int() {
  return integerCong ^ tausworthe ^ hurd;
}

Hurd288Engine & TripleRand::Hurd() 	       { return hurd; }

const Hurd288Engine & TripleRand::ConstHurd() const 
						{ return hurd; }

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

  os << " " << beginMarker << " ";
  os << HepStd::setprecision(20) << e.theSeed << " ";
  e.tausworthe.put( os );
  e.integerCong.put( os );
  os << e.ConstHurd();
  os << " " <<  endMarker  << " ";
  return os;
}

HepStd::istream & operator>>(HepStd::istream & is, TripleRand & e) {
  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,"TripleRand-begin")) {
    is.clear(HepStd::ios::badbit | is.rdstate());
    HepStd::cerr << "\nInput mispositioned or"
         << "\nTripleRand state description missing or"
         << "\nwrong engine type found." << HepStd::endl;
    return is;
  }
  is >> e.theSeed;
  e.tausworthe.get( is );
  e.integerCong.get( is );
  is >> e.Hurd();
  is >> HepStd::ws;
  is.width(MarkerLen);  
  is >> endMarker;
  if (strcmp(endMarker,"TripleRand-end")) {
    is.clear(HepStd::ios::badbit | is.rdstate());
    HepStd::cerr << "\nTripleRand state description incomplete."
         << "\nInput stream is probably mispositioned now." << HepStd::endl;
    return is;
  }
  return is;
}

//********************************************************************
//   Tausworthe
//********************************************************************

TripleRand::Tausworthe::Tausworthe() {
  words[0] = 1234567;
  for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
    words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
  } 
}

TripleRand::Tausworthe::Tausworthe(unsigned int seed) {
  words[0] = seed;
  for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
    words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
  }
}

TripleRand::Tausworthe::operator unsigned int() {

// Mathematically:  Consider a sequence of bits b[n].  Repeatedly form
// b[0]' = b[127] ^ b[97]; b[n]' = b[n-1].  This sequence will have a very
// long period (2**127-1 according to Tausworthe's work).

// The actual method used relies on the fact that the operations needed to
// form bit 0 for up to 96 iterations never depend on the results of the
// previous ones.  So you can actually compute many bits at once.  In fact
// you can compute 32 at once -- despite 127 - 97 < 32 -- but 24 was used in
// the method used in Canopy, where they wanted only single-precision float
// randoms.  I will do 32 here.

// When you do it this way, this looks disturbingly like the dread lagged XOR
// Fibonacci.  And indeed, it is a lagged Fibonacii, F(4,3, op) with the op
// being the XOR of a combination of shifts of the two numbers.  Although
// Tausworthe asserted excellent properties, I would be scared to death.
// However, the shifting and bit swapping really does randomize this in a
// serious way.

// Statements have been made to the effect that shift register sequences fail
// the parking lot test because they achieve randomness by multiple foldings,
// and this produces a characteristic pattern.  We observe that in this
// specific algorithm, which has a fairly long lever arm, the foldings become
// effectively random.  This is evidenced by the fact that the generator
// does pass the Diehard tests, including the parking lot test.

// To avoid shuffling of variables in memory, you either have to use circular
// pointers (and those give you ifs, which are also costly) or compute at least
// a few iterations at once.  We do the latter.  Although there is a possible
// trade of room for more speed, by computing and saving 256 instead of 128
// bits at once, I will stop at this level of optimization.

// To remind:  Each (32-bit) step takes the XOR of bits [127-96] with bits
// [95-64] and places it in bits [0-31].  But in the first step, we designate
// word0 as bits [0-31], in the second step, word 1 (since the bits it holds
// will no longer be needed), then word 2, then word 3.  After this, the
// stream contains 128 random bits which we will use as 4 valid 32-bit
// random numbers.

// Thus at the start of the first step, word[0] contains the newest (used)
// 32-bit random, and word[3] the oldest.   After four steps, word[0] again
// contains the newest (now unused) random, and word[3] the oldest.
// Bit 0 of word[0] is logically the newest bit, and bit 31 of word[3]
// the oldest.

  if (wordIndex <= 0) {
    for (wordIndex = 0; wordIndex < 4; ++wordIndex) {
      words[wordIndex] = ( (words[(wordIndex+1) & 3] << 1 ) |
                                   (words[wordIndex] >> 31)   )
                       ^ ( (words[(wordIndex+1) & 3] << 31) |
                                   (words[wordIndex] >>  1)   );
    }
  }
  return words[--wordIndex] & 0xffffffff;
}

void TripleRand::Tausworthe::put( HepStd::ostream & os ) const {
  char beginMarker[] = "Tausworthe-begin";
  char endMarker[]   = "Tausworthe-end";

  os << " " << beginMarker << " ";
  os << HepStd::setprecision(20);
  for (int i = 0; i < 4; ++i) {
    os << words[i] << " ";
  }
  os << wordIndex;
  os << " " <<  endMarker  << " ";
  os << HepStd::endl;
}

void TripleRand::Tausworthe::get( HepStd::istream & is ) {
  char beginMarker [MarkerLen];
  char endMarker   [MarkerLen];

  is >> HepStd::ws;
  is.width(MarkerLen);  
  is >> beginMarker;
  if (strcmp(beginMarker,"Tausworthe-begin")) {
    is.clear(HepStd::ios::badbit | is.rdstate());
    HepStd::cerr << "\nInput mispositioned or"
              << "\nTausworthe state description missing or"
              << "\nwrong engine type found." << HepStd::endl;
  }
  for (int i = 0; i < 4; ++i) {
    is >> words[i];
  }
  is >> wordIndex;
  is >> HepStd::ws;
  is.width(MarkerLen);  
  is >> endMarker;
  if (strcmp(endMarker,"Tausworthe-end")) {
    is.clear(HepStd::ios::badbit | is.rdstate());
    HepStd::cerr << "\nTausworthe state description incomplete."
              << "\nInput stream is probably mispositioned now." << HepStd::endl;
  }
}

//********************************************************************
//   IntegerCong
//********************************************************************

TripleRand::IntegerCong::IntegerCong() 
: state((unsigned int)3758656018U),
  multiplier(66565),
  addend(12341) 
{
}

TripleRand::IntegerCong::IntegerCong(unsigned int seed, int streamNumber)
: state(seed),
  multiplier(65536 + 1024 + 5 + (8 * 1017 * streamNumber)),
  addend(12341)
{
  // As to the multiplier, the following comment was made:
  // We want our multipliers larger than 2^16, and equal to
  // 1 mod 4 (for max. period), but not equal to 1 mod 8 
  // (for max. potency -- the smaller and higher dimension the 
  // stripes, the better)

  // All of these will have fairly long periods.  Depending on the choice
  // of stream number, some of these will be quite decent when considered
  // as independent random engines, while others will be poor.  Thus these
  // should not be used as stand-alone engines; but when combined with a
  // generator known to be good, they allow creation of millions of good
  // independent streams, without fear of two streams accidentally hitting
  // nearby places in the good random sequence.
}

TripleRand::IntegerCong::operator unsigned int() {
  return state = (state * multiplier + addend) & 0xffffffff;
}

void TripleRand::IntegerCong::put( HepStd::ostream & os ) const {
  char beginMarker[] = "IntegerCong-begin";
  char endMarker[]   = "IntegerCong-end";

  os << " " << beginMarker << " ";
  os << HepStd::setprecision(20);
  os << state << " " << multiplier << " " << addend;
  os << " " <<  endMarker  << " ";
  os << HepStd::endl;
}

void TripleRand::IntegerCong::get( HepStd::istream & is ) {
  char beginMarker [MarkerLen];
  char endMarker   [MarkerLen];

  is >> HepStd::ws;
  is.width(MarkerLen);  
  is >> beginMarker;
  if (strcmp(beginMarker,"IntegerCong-begin")) {
    is.clear(HepStd::ios::badbit | is.rdstate());
    HepStd::cerr << "\nInput mispositioned or"
              << "\nIntegerCong state description missing or"
              << "\nwrong engine type found." << HepStd::endl;
  }
  is >> state >> multiplier >> addend;
  is >> HepStd::ws;
  is.width(MarkerLen);  
  is >> endMarker;
  if (strcmp(endMarker,"IntegerCong-end")) {
    is.clear(HepStd::ios::badbit | is.rdstate());
    HepStd::cerr << "\nIntegerCong state description incomplete."
              << "\nInput stream is probably mispositioned now." << HepStd::endl;
  }
}