CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

Random/Random/RandExpZiggurat.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 #ifndef RandExpZiggurat_h
24 #define RandExpZiggurat_h 1
25 
26 #include "CLHEP/Random/defs.h"
27 #include "CLHEP/Random/Random.h"
28 
29 namespace CLHEP {
30 
35 class RandExpZiggurat : public HepRandom {
36 
37 public:
38 
39  inline RandExpZiggurat ( HepRandomEngine& anEngine, double mean=1.0 );
40  inline RandExpZiggurat ( HepRandomEngine* anEngine, double mean=1.0 );
41  // These constructors should be used to instantiate a RandExpZiggurat
42  // distribution object defining a local engine for it.
43  // The static generator will be skipped using the non-static methods
44  // defined below.
45  // If the engine is passed by pointer the corresponding engine object
46  // will be deleted by the RandExpZiggurat destructor.
47  // If the engine is passed by reference the corresponding engine object
48  // will not be deleted by the RandExpZiggurat destructor.
49 
50  virtual ~RandExpZiggurat();
51  // Destructor
52 
53  // Static methods to shoot random values using the static generator
54 
55  static float shoot() {return shoot(HepRandom::getTheEngine());};
56  static float shoot( float mean ) {return shoot(HepRandom::getTheEngine(),mean);};
57 
58  /* ENGINE IS INTRINSIC FLOAT
59  static double shoot() {return shoot(HepRandom::getTheEngine());};
60  static double shoot( double mean ) {return shoot(HepRandom::getTheEngine(),mean);};
61  */
62 
63  static void shootArray ( const int size, float* vect, float mean=1.0 );
64  static void shootArray ( const int size, double* vect, double mean=1.0 );
65 
66  // Static methods to shoot random values using a given engine
67  // by-passing the static generator.
68 
69  static inline float shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
70  static inline float shoot( HepRandomEngine* anEngine, float mean ) {return shoot(anEngine)*mean;};
71 
72  /* ENGINE IS INTRINSIC FLOAT
73  static inline double shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
74 
75  static inline double shoot( HepRandomEngine* anEngine, double mean ) {return shoot(anEngine)*mean;};
76  */
77 
78  static void shootArray ( HepRandomEngine* anEngine, const int size, float* vect, float mean=1.0 );
79  static void shootArray ( HepRandomEngine* anEngine, const int size, double* vect, double mean=1.0 );
80 
81  // Methods using the localEngine to shoot random values, by-passing
82  // the static generator.
83 
84  inline float fire() {return fire(defaultMean);};
85  inline float fire( float mean ) {return ziggurat_REXP(localEngine)*mean;};
86 
87  /* ENGINE IS INTRINSIC FLOAT
88  inline double fire() {return fire(defaultMean);};
89  inline double fire( double mean ) {return ziggurat_REXP(localEngine)*mean;};
90  */
91 
92  void fireArray ( const int size, float* vect );
93  void fireArray ( const int size, double* vect );
94  void fireArray ( const int size, float* vect, float mean );
95  void fireArray ( const int size, double* vect, double mean );
96 
97  virtual double operator()();
98  inline float operator()( float mean ) {return fire( mean );};
99 
100  // Save and restore to/from streams
101 
102  std::ostream & put ( std::ostream & os ) const;
103  std::istream & get ( std::istream & is );
104 
105  std::string name() const;
107 
108  static std::string distributionName() {return "RandExpZiggurat";}
109  // Provides the name of this distribution class
110 
111  static bool ziggurat_init();
112 protected:
114  // Ziggurat Original code:
116  //static unsigned long jz,jsr=123456789;
117  //
118  //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
119  //#define UNI (.5 + (signed) SHR3*.2328306e-9)
120  //#define IUNI SHR3
121  //
122  //static long hz;
123  //static unsigned long iz, kn[128], ke[256];
124  //static float wn[128],fn[128], we[256],fe[256];
125  //
126  //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
127  //#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix())
128 
129  static unsigned long kn[128], ke[256];
130  static float wn[128],fn[128], we[256],fe[256];
131 
132  static bool ziggurat_is_init;
133 
134  static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
135  static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
136  static inline float ziggurat_REXP(HepRandomEngine* anEngine) {
137  unsigned long jz=ziggurat_SHR3(anEngine);
138  unsigned long iz=jz&255;
139  return (jz<ke[iz]) ? jz*we[iz] : ziggurat_efix(jz,anEngine);
140  };
141  static float ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine);
142 
143 private:
144 
145  // Private copy constructor. Defining it here disallows use.
147 
148  HepRandomEngine* localEngine;
149  bool deleteEngine;
150  double defaultMean;
151 };
152 
153 } // namespace CLHEP
154 
155 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
156 // backwards compatibility will be enabled ONLY in CLHEP 1.9
157 using namespace CLHEP;
158 #endif
159 
160 namespace CLHEP {
161 
162 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine), deleteEngine(false), defaultMean(mean)
163 {
165 }
166 
167 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), deleteEngine(true), defaultMean(mean)
168 {
170 }
171 
172 } // namespace CLHEP
173 
174 #endif
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
static void shootArray(const int size, float *vect, float mean=1.0)
static float ziggurat_UNI(HepRandomEngine *anEngine)
virtual double operator()()
void fireArray(const int size, float *vect)
static float shoot(float mean)
RandExpZiggurat(HepRandomEngine &anEngine, double mean=1.0)
virtual double flat()=0
std::string name() const
HepRandomEngine & engine()
static float ziggurat_efix(unsigned long jz, HepRandomEngine *anEngine)
static float shoot(HepRandomEngine *anEngine)
std::ostream & put(std::ostream &os) const
static float shoot(HepRandomEngine *anEngine, float mean)
static float ziggurat_REXP(HepRandomEngine *anEngine)