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

Random/CLHEP/Random/RandExpZiggurat.h
Go to the documentation of this file.
1/*
2Code adapted from:
3http://www.jstatsoft.org/v05/i08/
4
5
6Original disclaimer:
7
8The ziggurat method for RNOR and REXP
9Combine the code below with the main program in which you want
10normal or exponential variates. Then use of RNOR in any expression
11will provide a standard normal variate with mean zero, variance 1,
12while use of REXP in any expression will provide an exponential variate
13with density exp(-x),x>0.
14Before using RNOR or REXP in your main, insert a command such as
15zigset(86947731 );
16with 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.)
18For details of the method, see Marsaglia and Tsang, "The ziggurat method
19for 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
29namespace CLHEP {
30
35class RandExpZiggurat : public HepRandom {
36
37public:
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
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();
112protected:
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
143private:
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
157using namespace CLHEP;
158#endif
159
160namespace CLHEP {
161
162inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine), deleteEngine(false), defaultMean(mean)
163{
165}
166
167inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), deleteEngine(true), defaultMean(mean)
168{
170}
171
172} // namespace CLHEP
173
174#endif
virtual double flat()=0
static HepRandomEngine * getTheEngine()
Definition Random.cc:166
static float ziggurat_REXP(HepRandomEngine *anEngine)
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
HepRandomEngine & engine()
static float ziggurat_efix(unsigned long jz, HepRandomEngine *anEngine)
static void shootArray(const int size, double *vect, double mean=1.0)
static bool ziggurat_init()
std::istream & get(std::istream &is)
static void shootArray(HepRandomEngine *anEngine, const int size, float *vect, float mean=1.0)
static float shoot(HepRandomEngine *anEngine)
static void shootArray(const int size, float *vect, float mean=1.0)
void fireArray(const int size, float *vect, float mean)
RandExpZiggurat(HepRandomEngine &anEngine, double mean=1.0)
static void shootArray(HepRandomEngine *anEngine, const int size, double *vect, double mean=1.0)
void fireArray(const int size, double *vect, double mean)
static float shoot(HepRandomEngine *anEngine, float mean)
std::string name() const
void fireArray(const int size, double *vect)
RandExpZiggurat(HepRandomEngine *anEngine, double mean=1.0)
static float ziggurat_efix(unsigned long jz, HepRandomEngine *anEngine)
void fireArray(const int size, float *vect)
static float ziggurat_UNI(HepRandomEngine *anEngine)
std::ostream & put(std::ostream &os) const
virtual double operator()()