source: Sophya/trunk/SophyaLib/BaseTools/randr48.cc@ 3611

#include "machdefs.h"
#include <math.h> 
#include <stdlib.h>
#include "thsafeop.h"
#include "fiondblock.h"

#include "randr48.h"

namespace SOPHYA {

DR48RandGen::DR48RandGen(long int  seed)
{
  srand48(seed);
}

DR48RandGen::~DR48RandGen()
{
}

void DR48RandGen::SetSeed(long int  seed)
{
  srand48(seed);
}

void DR48RandGen::SetSeed(uint_2 seed[3])
{
  seed48(seed);
}
void DR48RandGen::GetSeed(uint_2 seed[3])
{
  uint_2 *p, seed_dummy[3] = {0,0,0};
  p = seed48(seed_dummy);
  memcpy(seed,p,3*sizeof(uint_2));
  // on re-initialise a ce qui etait avant
  seed48(seed);
}

r_8 DR48RandGen::Next()
{
  return drand48();
}

void DR48RandGen::AutoInit(int lp)
{
  vector<uint_2> seed;
  GenerateSeedVector(0,seed,lp);
  uint_2 s[3] = {seed[0],seed[1],seed[2]};
  SetSeed(s);
}

//------------------------------------------------------------
//------------------------------------------------------------
//------------------------------------------------------------
/*!
   \class ThSDR48RandGen
   \ingroup BaseTools
   \brief Random number generator

   This class is a thread-safe random number generator.
   Its PPF handler can be used to save the complete state of the class and the underlying
   random number generator used.

   \sa SOPHYA::ObjFileIO<ThSDR48RandGen>

*/

// Objet statique global pour gestion de lock entre threads
static ThSafeOp* ths_rand = NULL;

ThSDR48RandGen::ThSDR48RandGen(size_t n, bool tsafe) 
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  if (tsafe) {   // thread-safe
    fg_nothrsafe = false;
    if (n < 1) n = 1024;
    rseq_.ReSize(n, false);
    idx_ = n;
  }
  else {   // NOT thread-safe
    fg_nothrsafe = true;
    idx_ = 1;
  }
}

ThSDR48RandGen::ThSDR48RandGen(ThSDR48RandGen const & rg)
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  if (!rg.fg_nothrsafe) { // thread-safe
    fg_nothrsafe = false;
    rseq_.ReSize(rg.rseq_.Size(), false);
    idx_ = rseq_.Size();
  }
  else {   // NOT thread-safe
    fg_nothrsafe = true;
    idx_ = 1;
  }
}


ThSDR48RandGen::~ThSDR48RandGen(void) 
{
  // rien a faire 
}

void ThSDR48RandGen::SetBuffSize(size_t n)
// redimensionnement du buffer
{
  if(fg_nothrsafe) return;
  if (n < 1) n = 1024;
  rseq_.ReSize(n, false);
  idx_ = n;
}


void ThSDR48RandGen::SetSeed(long int seed)
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  ths_rand->lock();
  DR48RandGen::SetSeed(seed);
  ths_rand->unlock();
  return;
}

void ThSDR48RandGen::SetSeed(uint_2 seed[3])
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  ths_rand->lock();
  SetSeed_P(seed);
  ths_rand->unlock();
}

void ThSDR48RandGen::GetSeed(uint_2 seed[3])
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  ths_rand->lock();
  GetSeed_P(seed);
  ths_rand->unlock();
  return;
}

void ThSDR48RandGen::SetSeed_P(uint_2 seed[3])
{
  DR48RandGen::SetSeed(seed);
}

void ThSDR48RandGen::GetSeed_P(uint_2 seed[3])
{
  DR48RandGen::GetSeed(seed);
}

void ThSDR48RandGen::AutoInit(int lp)
{
  vector<uint_2> seed;
  GenerateSeedVector(0,seed,lp);
  uint_2 s[3] = {seed[0],seed[1],seed[2]};
  SetSeed(s);
}


void ThSDR48RandGen::GenSeq(void)
{
  ths_rand->lock();
  for(size_t k=0; k<rseq_.Size(); k++)  rseq_(k) = drand48();
  ths_rand->unlock();
  idx_ = 0;
}
 
r_8 ThSDR48RandGen::Next() 
{
  if (rseq_.Size() == 0)  return drand48(); 
  else { 
    if(idx_==rseq_.Size()) GenSeq(); 
    return(rseq_(idx_++));
  }
}

//----------------------------------------------------------
// Classe pour la gestion de persistance
// ObjFileIO<DR48RandGen>
//----------------------------------------------------------

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<DR48RandGen>::WriteSelf(POutPersist& s) const
{
  if (dobj == NULL)
    throw NullPtrError("ObjFileIO<DR48RandGen>::WriteSelf() dobj=NULL");
  uint_2 seed[3];
  dobj->GetSeed(seed);
  s.Put(seed,3);
  return;
}

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<DR48RandGen>::ReadSelf(PInPersist& s)
{
  uint_2 seed[3];
  s.Get(seed,3);
  if(dobj == NULL) dobj = new DR48RandGen();
  dobj->SetSeed(seed);
  return;
}

//----------------------------------------------------------
// Classe pour la gestion de persistance
// ObjFileIO<ThSDR48RandGen>
//----------------------------------------------------------

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<ThSDR48RandGen>::WriteSelf(POutPersist& s) const
{
  if (dobj == NULL)
    throw NullPtrError("ObjFileIO<ThSDR48RandGen>::WriteSelf() dobj=NULL");
  ths_rand->lock();  // thread-safety
  uint_4 itab[6];
  //itab : [0]: version, [1,2,3] = srand48 state/seed  , [4,5] = reserved for future use
  itab[0] = 1;
  // On recupere et on ecrit ds le PPF l'etat du generateur aleatoire
  uint_2 seed_16v[3];
  dobj->GetSeed_P(seed_16v);
  for(int i=0; i<3; i++)  itab[i+1] = seed_16v[i];
  itab[4] = 0;
  s.Put(itab, 6);
  uint_8 sz = dobj->rseq_.Size();
  s.Put(sz);  // Taille du tableau intermediaire
  uint_8 ix = dobj->idx_;
  s.Put(ix);  // valeur de l'index 

  if (dobj->rseq_.Size() > 0)  s << dobj->rseq_;  // On ecrit le tableau (NDataBlock) si necessaire
  ths_rand->unlock(); // thread-safety
  return;
}

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<ThSDR48RandGen>::ReadSelf(PInPersist& s)
{
  uint_4 itab[6];
  //itab : [0]: version, [1,2,3] = srand48 state/seed  , [4] = reserved for future use
  s.Get(itab, 6);
  uint_8 sz,ix;
  s.Get(sz);  // Taille du tableau intermediaire
  s.Get(ix);  // Taille du tableau intermediaire

  if (dobj == NULL) dobj = new ThSDR48RandGen(sz, (sz>0)?true:false);
  dobj->idx_ = ix;
  if (sz > 0) { 
    s >> dobj->rseq_;  // On lit le tableau (NDataBlock) si necessaire
    dobj->fg_nothrsafe = false;
  }
  else {  // Objet lu est NON thread-safe, taille_tableau rseq_ = 0
    dobj->fg_nothrsafe = true;
    if (dobj->rseq_.Size() > 0)  dobj->rseq_.Dealloc(); 
  }
  // On initialise l'etat du generateur aleatoire avec les valeurs lues
  uint_2 seed_16v[3];
  //NON ? pourquoi faire GetSeed ? : dobj->GetSeed_P(seed_16v);
  for(int i=0; i<3; i++)  seed_16v[i] = itab[i+1];
  dobj->SetSeed(seed_16v);
  return;
}

// ---------------------------------------------------------
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<DR48RandGen>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ObjFileIO<DR48RandGen>;
#endif

// ---------------------------------------------------------
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<ThSDR48RandGen>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ObjFileIO<ThSDR48RandGen>;
#endif
// ---------------------------------------------------------


}  /* namespace SOPHYA */