RandGaussZiggurat.h
Go to the documentation of this file.
1 /*
2 Code adapted from:
3 http://www.jstatsoft.org/v05/i08/
4 
5 
6 Original disclaimer:
7 
8 The ziggurat method for RNOR and REXP
9 Combine the code below with the main program in which you want
10 normal or exponential variates. Then use of RNOR in any expression
11 will provide a standard normal variate with mean zero, variance 1,
12 while use of REXP in any expression will provide an exponential variate
13 with density exp(-x),x>0.
14 Before using RNOR or REXP in your main, insert a command such as
15 zigset(86947731 );
16 with your own choice of seed value>0, rather than 86947731.
17 (If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
18 For details of the method, see Marsaglia and Tsang, "The ziggurat method
19 for generating random variables", Journ. Statistical Software.
20 
21 */
22 
23 
24 #ifndef RandGaussZiggurat_h
25 #define RandGaussZiggurat_h 1
26 
27 #include "cmath"
28 #include "CLHEP/Random/RandGauss.h"
30 
31 namespace CLHEP {
32 
33 /**
34  * @author ATLAS
35  * @ingroup random
36  */
37 class RandGaussZiggurat : public RandGauss {
38 
39 public:
40 
41  inline RandGaussZiggurat ( HepRandomEngine& anEngine, double mean=0.0, double stdDev=1.0 );
42  inline RandGaussZiggurat ( HepRandomEngine* anEngine, double mean=0.0, double stdDev=1.0 );
43 
44  // Destructor
45  virtual ~RandGaussZiggurat();
46 
47  // Static methods to shoot random values using the static generator
48 
49  static inline float shoot() {return ziggurat_RNOR(HepRandom::getTheEngine());};
50  static inline float shoot( float mean, float stdDev ) {return shoot()*stdDev + mean;};
51 
52  static void shootArray ( const int size, float* vect, float mean=0.0, float stdDev=1.0 );
53  static void shootArray ( const int size, double* vect, double mean=0.0, double stdDev=1.0 );
54 
55  // Static methods to shoot random values using a given engine
56  // by-passing the static generator.
57 
58  static inline float shoot( HepRandomEngine* anotherEngine ) {return ziggurat_RNOR(anotherEngine);};
59  static inline float shoot( HepRandomEngine* anotherEngine, float mean, float stdDev ) {return shoot(anotherEngine)*stdDev + mean;};
60 
61  static void shootArray ( HepRandomEngine* anotherEngine, const int size, float* vect, float mean=0.0, float stdDev=1.0 );
62  static void shootArray ( HepRandomEngine* anotherEngine, const int size, double* vect, double mean=0.0, double stdDev=1.0 );
63 
64  // Instance methods using the localEngine to instead of the static
65  // generator, and the default mean and stdDev established at construction
66 
67  inline float fire() {return ziggurat_RNOR(localEngine.get()) * defaultStdDev + defaultMean;};
68 
69  inline float fire ( float mean, float stdDev ) {return ziggurat_RNOR(localEngine.get()) * stdDev + mean;};
70 
71  void fireArray ( const int size, float* vect);
72  void fireArray ( const int size, double* vect);
73  void fireArray ( const int size, float* vect, float mean, float stdDev );
74  void fireArray ( const int size, double* vect, double mean, double stdDev );
75 
76  virtual double operator()();
77  virtual double operator()( double mean, double stdDev );
78 
79  // Save and restore to/from streams
80 
81  std::ostream & put ( std::ostream & os ) const;
82  std::istream & get ( std::istream & is );
83 
84  std::string name() const;
86 
87  static std::string distributionName() {return "RandGaussZiggurat";}
88  // Provides the name of this distribution class
89 
90  static bool ziggurat_init();
91 protected:
92  //////////////////////////
93  // Ziggurat Original code:
94  //////////////////////////
95  //static unsigned long jz,jsr=123456789;
96  //
97  //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
98  //#define UNI (.5 + (signed) SHR3*.2328306e-9)
99  //#define IUNI SHR3
100  //
101  //static long hz;
102  //static unsigned long iz, kn[128], ke[256];
103  //static float wn[128],fn[128], we[256],fe[256];
104  //
105  //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
106  //#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix())
107 
108  static CLHEP_THREAD_LOCAL unsigned long kn[128], ke[256];
109  static CLHEP_THREAD_LOCAL float wn[128],fn[128], we[256],fe[256];
110 
112 
113  static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
114  static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
115  static inline float ziggurat_RNOR(HepRandomEngine* anEngine) {
116  if(!ziggurat_is_init) ziggurat_init();
117  long hz=(signed)ziggurat_SHR3(anEngine);
118  unsigned long iz=hz&127;
119  return ((unsigned long)std::abs(hz)<kn[iz]) ? hz*wn[iz] : ziggurat_nfix(hz,anEngine);
120  };
121  static float ziggurat_nfix(long hz,HepRandomEngine* anEngine);
122 
123 private:
124 
125  // Private copy constructor. Defining it here disallows use.
127 
128  // All the engine info, and the default mean and sigma, are in the RandGauss
129  // base class.
130 
131 };
132 
133 } // namespace CLHEP
134 
135 namespace CLHEP {
136 
137 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine & anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
138 {
139 }
140 
141 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine * anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
142 {
143 }
144 
145 } // namespace CLHEP
146 
147 #endif
RandGaussZiggurat(HepRandomEngine &anEngine, double mean=0.0, double stdDev=1.0)
static CLHEP_THREAD_LOCAL bool ziggurat_is_init
static CLHEP_THREAD_LOCAL unsigned long ke[256]
static CLHEP_THREAD_LOCAL float fe[256]
std::shared_ptr< HepRandomEngine > localEngine
Definition: RandGauss.h:155
double defaultMean
Definition: RandGauss.h:152
static CLHEP_THREAD_LOCAL float fn[128]
static CLHEP_THREAD_LOCAL float wn[128]
std::string name() const
static float ziggurat_UNI(HepRandomEngine *anEngine)
static float shoot(HepRandomEngine *anotherEngine, float mean, float stdDev)
virtual double operator()()
float abs(float number)
Definition: d0nt_math.hpp:39
void fireArray(const int size, float *vect)
static float ziggurat_RNOR(HepRandomEngine *anEngine)
virtual double flat()=0
static float shoot(HepRandomEngine *anotherEngine)
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
std::ostream & put(std::ostream &os) const
static bool ziggurat_init()
Float_t d
Definition: plot.C:236
HepRandomEngine & engine()
static void shootArray(const int size, float *vect, float mean=0.0, float stdDev=1.0)
#define CLHEP_THREAD_LOCAL
Definition: thread_local.h:25
static CLHEP_THREAD_LOCAL float we[256]
float fire(float mean, float stdDev)
static float shoot(float mean, float stdDev)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static float ziggurat_nfix(long hz, HepRandomEngine *anEngine)
static CLHEP_THREAD_LOCAL unsigned long kn[128]
static std::string distributionName()
double defaultStdDev
Definition: RandGauss.h:153
static HepRandomEngine * getTheEngine()