| 1 | #ifndef STSRAND_H_SEEN | 
|---|
| 2 | #define STSRAND_H_SEEN | 
|---|
| 3 |  | 
|---|
| 4 | // Classe RandomGenerator | 
|---|
| 5 | // Generateur aleatoire compatible multi-thread | 
|---|
| 6 | // | 
|---|
| 7 | // R. Ansari          (C) UPS+LAL IN2P3/CNRS  2007 | 
|---|
| 8 | // C. Magneville      (C) DAPNIA/SPP  CEA     2007 | 
|---|
| 9 |  | 
|---|
| 10 |  | 
|---|
| 11 | #include "machdefs.h" | 
|---|
| 12 | #include "objfio.h"   // Pour rendre la classe PPersist | 
|---|
| 13 |  | 
|---|
| 14 | #include "ndatablock.h" | 
|---|
| 15 |  | 
|---|
| 16 | namespace SOPHYA { | 
|---|
| 17 |  | 
|---|
| 18 | class STSRandGen : public AnyDataObj { | 
|---|
| 19 |  | 
|---|
| 20 | public: | 
|---|
| 21 | STSRandGen(size_t n = 1024, bool tsafe=true); | 
|---|
| 22 | STSRandGen(STSRandGen const & rg); | 
|---|
| 23 |  | 
|---|
| 24 | void SetBuffSize(size_t n); | 
|---|
| 25 |  | 
|---|
| 26 | virtual ~STSRandGen(); | 
|---|
| 27 |  | 
|---|
| 28 | /*! \brief Automatic initialization using the current time */ | 
|---|
| 29 | static  void AutoInit(int lp=0); | 
|---|
| 30 | /*! \brief Initialization using the specified seed  */ | 
|---|
| 31 | static  void Init(long seed_val, int lp=0); | 
|---|
| 32 | /*! \brief Initialization using the specified seed  */ | 
|---|
| 33 | static  void Init(unsigned short seed_16v[3], int lp=0); | 
|---|
| 34 | /*! \brief Getting the current seed state */ | 
|---|
| 35 | static  void GetSeed(unsigned short seed_16v[3], int lp=0); | 
|---|
| 36 |  | 
|---|
| 37 | /*! \brief Return a random number \b r with a flat distribution  0<=r<1 */ | 
|---|
| 38 | inline r_8 Flat01() {return Next();} | 
|---|
| 39 |  | 
|---|
| 40 | /*! \brief Return a random number \b r with a flat distribution  -1<=r<1 */ | 
|---|
| 41 | inline r_8 Flatpm1() {return 2.*Next()-1.;} | 
|---|
| 42 |  | 
|---|
| 43 | /*! \brief Return a random number following a gaussian distribution (sigma=1, mean=0)*/ | 
|---|
| 44 | virtual r_8 Gaussian(); | 
|---|
| 45 |  | 
|---|
| 46 | /*! \brief Return a random number following a gaussian distribution "sigma", (mean=0)*/ | 
|---|
| 47 | virtual r_8 Gaussian(double sigma)  {return sigma*Gaussian();} | 
|---|
| 48 | /*! \brief Return a random number following a gaussian distribution with mean=mu, and sigma */ | 
|---|
| 49 | virtual r_8 Gaussian(double sigma,double mu) {return sigma*Gaussian()+mu;} | 
|---|
| 50 |  | 
|---|
| 51 | /*! \brief alias for Gaussian() */ | 
|---|
| 52 | inline r_8 NorRand() {return Gaussian();} | 
|---|
| 53 |  | 
|---|
| 54 | /*! \brief Return a random number following a poisson distribution with mean mu */ | 
|---|
| 55 | virtual  uint_8 Poisson(double mu, double mumax); | 
|---|
| 56 |  | 
|---|
| 57 | /*! \brief Return a random number following a poisson distribution with mean mu */ | 
|---|
| 58 | inline uint_8 Poisson(double mu) {return Poisson(mu, -1.);} | 
|---|
| 59 |  | 
|---|
| 60 | /*! \brief Return a random number following a poisson distribution with mean mu */ | 
|---|
| 61 | inline uint_8  PoissRandLimit(double mu,double mumax=10.) {return Poisson(mu, mumax);} | 
|---|
| 62 |  | 
|---|
| 63 | /*! \brief Returns a random complex number such that real and imaginary parts are gaussians with variance sig^2 */ | 
|---|
| 64 | inline complex< r_8 > ComplexGaussRan(void) | 
|---|
| 65 | {return complex< r_8 >(NorRand(),NorRand());} | 
|---|
| 66 | inline complex< r_8 > ComplexGaussRan(double sig) | 
|---|
| 67 | {return complex< r_8 >(sig*NorRand(),sig*NorRand());} | 
|---|
| 68 | /*! \brief Returns the module of a random complex number generated by ComplexGaussRan */ | 
|---|
| 69 | inline double ModComplexGaussRan(double sig=1.) | 
|---|
| 70 | {double r=-log(1.-Flat01()); return sig*sqrt(2.*r);} | 
|---|
| 71 |  | 
|---|
| 72 |  | 
|---|
| 73 | //  Pour la gestion de persistance PPF | 
|---|
| 74 | friend class ObjFileIO<STSRandGen> ; | 
|---|
| 75 |  | 
|---|
| 76 | protected: | 
|---|
| 77 | // ---- protected data members | 
|---|
| 78 | NDataBlock<r_8>  rseq_; | 
|---|
| 79 | size_t idx_; | 
|---|
| 80 | bool fg_nothrsafe;  // if true --> NOT thread-safe | 
|---|
| 81 | // ---- protected methods | 
|---|
| 82 | void GenSeq(); | 
|---|
| 83 | inline r_8 Next() | 
|---|
| 84 | { | 
|---|
| 85 | if (rseq_.Size() == 0)  return drand48(); | 
|---|
| 86 | else { | 
|---|
| 87 | if(idx_==rseq_.Size()) GenSeq(); | 
|---|
| 88 | return(rseq_(idx_++)); | 
|---|
| 89 | } | 
|---|
| 90 | } | 
|---|
| 91 | // Non thread-safe version of Init() and GetSeed() | 
|---|
| 92 | static  void Init_P(unsigned short seed_16v[3]); | 
|---|
| 93 | static  void GetSeed_P(unsigned short seed_16v[3]); | 
|---|
| 94 |  | 
|---|
| 95 | };  // Fin de la classe STSRandGen | 
|---|
| 96 |  | 
|---|
| 97 | // Classe pour la gestion de persistance PPF :  ObjFileIO<STSRandGen> | 
|---|
| 98 |  | 
|---|
| 99 | /*! Writes the random generator object state in the POutPersist stream \b os */ | 
|---|
| 100 | inline POutPersist& operator << (POutPersist& os, STSRandGen & obj) | 
|---|
| 101 | { ObjFileIO<STSRandGen> fio(&obj);  fio.Write(os);  return(os); } | 
|---|
| 102 | /*! Reads the random generator object state from the PInPersist stream \b is */ | 
|---|
| 103 | inline PInPersist& operator >> (PInPersist& is, STSRandGen & obj) | 
|---|
| 104 | { ObjFileIO<STSRandGen> fio(&obj); is.SkipToNextObject(); fio.Read(is); return(is); } | 
|---|
| 105 |  | 
|---|
| 106 | } /* namespace SOPHYA */ | 
|---|
| 107 |  | 
|---|
| 108 | #endif | 
|---|