| 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 
 | 
|---|