CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
RandGaussZiggurat.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
24#ifndef RandGaussZiggurat_h
25#define RandGaussZiggurat_h 1
26
27#include "cmath"
28#include "CLHEP/Random/defs.h"
31
32namespace CLHEP {
33
39
40public:
41
42 inline RandGaussZiggurat ( HepRandomEngine& anEngine, double mean=0.0, double stdDev=1.0 );
43 inline RandGaussZiggurat ( HepRandomEngine* anEngine, double mean=0.0, double stdDev=1.0 );
44
45 // Destructor
47
48 // Static methods to shoot random values using the static generator
49
50 static inline float shoot() {return ziggurat_RNOR(HepRandom::getTheEngine());};
51 static inline float shoot( float mean, float stdDev ) {return shoot()*stdDev + mean;};
52
53 static void shootArray ( const int size, float* vect, float mean=0.0, float stdDev=1.0 );
54 static void shootArray ( const int size, double* vect, double mean=0.0, double stdDev=1.0 );
55
56 // Static methods to shoot random values using a given engine
57 // by-passing the static generator.
58
59 static inline float shoot( HepRandomEngine* anotherEngine ) {return ziggurat_RNOR(anotherEngine);};
60 static inline float shoot( HepRandomEngine* anotherEngine, float mean, float stdDev ) {return shoot(anotherEngine)*stdDev + mean;};
61
62 static void shootArray ( HepRandomEngine* anotherEngine, const int size, float* vect, float mean=0.0, float stdDev=1.0 );
63 static void shootArray ( HepRandomEngine* anotherEngine, const int size, double* vect, double mean=0.0, double stdDev=1.0 );
64
65 // Instance methods using the localEngine to instead of the static
66 // generator, and the default mean and stdDev established at construction
67
68 inline float fire() {return ziggurat_RNOR(localEngine.get()) * defaultStdDev + defaultMean;};
69
70 inline float fire ( float mean, float stdDev ) {return ziggurat_RNOR(localEngine.get()) * stdDev + mean;};
71
72 void fireArray ( const int size, float* vect);
73 void fireArray ( const int size, double* vect);
74 void fireArray ( const int size, float* vect, float mean, float stdDev );
75 void fireArray ( const int size, double* vect, double mean, double stdDev );
76
77 virtual double operator()();
78 virtual double operator()( double mean, double stdDev );
79
80 // Save and restore to/from streams
81
82 std::ostream & put ( std::ostream & os ) const;
83 std::istream & get ( std::istream & is );
84
85 std::string name() const;
87
88 static std::string distributionName() {return "RandGaussZiggurat";}
89 // Provides the name of this distribution class
90
91 static bool ziggurat_init();
92protected:
94 // Ziggurat Original code:
96 //static unsigned long jz,jsr=123456789;
97 //
98 //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
99 //#define UNI (.5 + (signed) SHR3*.2328306e-9)
100 //#define IUNI SHR3
101 //
102 //static long hz;
103 //static unsigned long iz, kn[128], ke[256];
104 //static float wn[128],fn[128], we[256],fe[256];
105 //
106 //#define RNOR (hz=SHR3, iz=hz&127, (std::fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
107 //#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix())
108
109 static CLHEP_THREAD_LOCAL unsigned long kn[128], ke[256];
110 static CLHEP_THREAD_LOCAL float wn[128],fn[128], we[256],fe[256];
111
113
114 static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
115 static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
116 static inline float ziggurat_RNOR(HepRandomEngine* anEngine) {
118 long hz=(signed)ziggurat_SHR3(anEngine);
119 unsigned long iz=hz&127;
120 return ((unsigned long)std::abs(hz)<kn[iz]) ? hz*wn[iz] : ziggurat_nfix(hz,anEngine);
121 };
122 static float ziggurat_nfix(long hz,HepRandomEngine* anEngine);
123
124private:
125
126 // Private copy constructor. Defining it here disallows use.
128
129 // All the engine info, and the default mean and sigma, are in the RandGauss
130 // base class.
131
132};
133
134} // namespace CLHEP
135
136#ifdef ENABLE_BACKWARDS_COMPATIBILITY
137// backwards compatibility will be enabled ONLY in CLHEP 1.9
138using namespace CLHEP;
139#endif
140
141namespace CLHEP {
142
143RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine & anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
144{
145}
146
147RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine * anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
148{
149}
150
151} // namespace CLHEP
152
153#endif
virtual double flat()=0
static HepRandomEngine * getTheEngine()
static float shoot(float mean, float stdDev)
void fireArray(const int size, double *vect, double mean, double stdDev)
std::istream & get(std::istream &is)
static CLHEP_THREAD_LOCAL float fn[128]
void fireArray(const int size, float *vect)
static CLHEP_THREAD_LOCAL unsigned long ke[256]
static void shootArray(HepRandomEngine *anotherEngine, const int size, double *vect, double mean=0.0, double stdDev=1.0)
static float ziggurat_UNI(HepRandomEngine *anEngine)
static bool ziggurat_init()
void fireArray(const int size, double *vect)
std::ostream & put(std::ostream &os) const
static void shootArray(HepRandomEngine *anotherEngine, const int size, float *vect, float mean=0.0, float stdDev=1.0)
static CLHEP_THREAD_LOCAL float fe[256]
static CLHEP_THREAD_LOCAL unsigned long kn[128]
static float ziggurat_RNOR(HepRandomEngine *anEngine)
virtual double operator()(double mean, double stdDev)
float fire(float mean, float stdDev)
void fireArray(const int size, float *vect, float mean, float stdDev)
static CLHEP_THREAD_LOCAL float we[256]
HepRandomEngine & engine()
static float shoot(HepRandomEngine *anotherEngine, float mean, float stdDev)
RandGaussZiggurat(HepRandomEngine &anEngine, double mean=0.0, double stdDev=1.0)
static float shoot(HepRandomEngine *anotherEngine)
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
std::string name() const
virtual double operator()()
static void shootArray(const int size, double *vect, double mean=0.0, double stdDev=1.0)
static CLHEP_THREAD_LOCAL bool ziggurat_is_init
static void shootArray(const int size, float *vect, float mean=0.0, float stdDev=1.0)
static std::string distributionName()
static float ziggurat_nfix(long hz, HepRandomEngine *anEngine)
static CLHEP_THREAD_LOCAL float wn[128]
double defaultStdDev
Definition RandGauss.h:154
RandGauss(HepRandomEngine &anEngine, double mean=0.0, double stdDev=1.0)
Definition RandGauss.icc:17
std::shared_ptr< HepRandomEngine > localEngine
Definition RandGauss.h:156
#define CLHEP_THREAD_LOCAL