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

RandGaussZiggurat.cc
Go to the documentation of this file.
1 #include "CLHEP/Random/defs.h"
2 #include "CLHEP/Random/RandGaussZiggurat.h"
3 #include "CLHEP/Units/PhysicalConstants.h"
4 #include <iostream>
5 #include <cmath> // for log()
6 
7 namespace CLHEP {
8 
10 //bool RandGaussZiggurat::ziggurat_is_init=false;
11 unsigned long RandGaussZiggurat::kn[128], RandGaussZiggurat::ke[256];
13 
15 
17 }
18 
19 std::string RandGaussZiggurat::name() const
20 {
21  return "RandGaussZiggurat";
22 }
23 
25 {
26  const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
27  double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
28  double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
29  int i;
30 
31 /* Set up tables for RNOR */
32  q=vn/exp(-.5*dn*dn);
33  kn[0]=(unsigned long)((dn/q)*rzm1);
34  kn[1]=0;
35 
36  wn[0]=q/rzm1;
37  wn[127]=dn/rzm1;
38 
39  fn[0]=1.;
40  fn[127]=exp(-.5*dn*dn);
41 
42  for(i=126;i>=1;i--) {
43  dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
44  kn[i+1]=(unsigned long)((dn/tn)*rzm1);
45  tn=dn;
46  fn[i]=exp(-.5*dn*dn);
47  wn[i]=dn/rzm1;
48  }
49 
50 /* Set up tables for REXP */
51  q = ve/exp(-de);
52  ke[0]=(unsigned long)((de/q)*rzm2);
53  ke[1]=0;
54 
55  we[0]=q/rzm2;
56  we[255]=de/rzm2;
57 
58  fe[0]=1.;
59  fe[255]=exp(-de);
60 
61  for(i=254;i>=1;i--) {
62  de=-log(ve/de+exp(-de));
63  ke[i+1]= (unsigned long)((de/te)*rzm2);
64  te=de;
65  fe[i]=exp(-de);
66  we[i]=de/rzm2;
67  }
68  ziggurat_is_init=true;
69 
70  //std::cout<<"Done RandGaussZiggurat::ziggurat_init()"<<std::endl;
71 
72  return true;
73 }
74 
76 {
77  const float r = 3.442620f; /* The start of the right tail */
78  float x, y;
79  unsigned long iz=hz&127;
80  for(;;)
81  {
82  x=hz*wn[iz]; /* iz==0, handles the base strip */
83  if(iz==0) {
84  do {
85  /* change to (1.0 - UNI) as argument to log(), because CLHEP generates [0,1),
86  while the original UNI generates (0,1] */
87  x=-log(1.0 - ziggurat_UNI(anEngine))*0.2904764; /* .2904764 is 1/r */
88  y=-log(1.0 - ziggurat_UNI(anEngine));
89  } while(y+y<x*x);
90  return (hz>0)? r+x : -r-x;
91  }
92  /* iz>0, handle the wedges of other strips */
93  if( fn[iz]+(1.0 - ziggurat_UNI(anEngine))*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) return x;
94 
95  /* initiate, try to exit for(;;) for loop*/
96  hz=(signed)ziggurat_SHR3(anEngine);
97  iz=hz&127;
98  if((unsigned long)abs(hz)<kn[iz]) return (hz*wn[iz]);
99  }
100 }
101 
104 }
105 
106 double RandGaussZiggurat::operator()( double mean, double stdDev ) {
107  return ziggurat_RNOR(localEngine.get()) * stdDev + mean;
108 }
109 
110 void RandGaussZiggurat::shootArray( const int size, float* vect, float mean, float stdDev )
111 {
112  for (int i=0; i<size; ++i) {
113  vect[i] = shoot(mean,stdDev);
114  }
115 }
116 
117 void RandGaussZiggurat::shootArray( const int size, double* vect, double mean, double stdDev )
118 {
119  for (int i=0; i<size; ++i) {
120  vect[i] = shoot(mean,stdDev);
121  }
122 }
123 
124 void RandGaussZiggurat::shootArray( HepRandomEngine* anEngine, const int size, float* vect, float mean, float stdDev )
125 {
126  for (int i=0; i<size; ++i) {
127  vect[i] = shoot(anEngine,mean,stdDev);
128  }
129 }
130 
131 void RandGaussZiggurat::shootArray( HepRandomEngine* anEngine, const int size, double* vect, double mean, double stdDev )
132 {
133  for (int i=0; i<size; ++i) {
134  vect[i] = shoot(anEngine,mean,stdDev);
135  }
136 }
137 
138 void RandGaussZiggurat::fireArray( const int size, float* vect)
139 {
140  for (int i=0; i<size; ++i) {
141  vect[i] = fire( defaultMean, defaultStdDev );
142  }
143 }
144 
145 void RandGaussZiggurat::fireArray( const int size, double* vect)
146 {
147  for (int i=0; i<size; ++i) {
148  vect[i] = fire( defaultMean, defaultStdDev );
149  }
150 }
151 
152 void RandGaussZiggurat::fireArray( const int size, float* vect, float mean, float stdDev )
153 {
154  for (int i=0; i<size; ++i) {
155  vect[i] = fire( mean, stdDev );
156  }
157 }
158 
159 void RandGaussZiggurat::fireArray( const int size, double* vect, double mean, double stdDev )
160 {
161  for (int i=0; i<size; ++i) {
162  vect[i] = fire( mean, stdDev );
163  }
164 }
165 
166 std::ostream & RandGaussZiggurat::put ( std::ostream & os ) const {
167  int pr=os.precision(20);
168  os << " " << name() << "\n";
169  RandGauss::put(os);
170  os.precision(pr);
171  return os;
172 }
173 
174 std::istream & RandGaussZiggurat::get ( std::istream & is ) {
175  std::string inName;
176  is >> inName;
177  if (inName != name()) {
178  is.clear(std::ios::badbit | is.rdstate());
179  std::cerr << "Mismatch when expecting to read state of a "
180  << name() << " distribution\n"
181  << "Name found was " << inName
182  << "\nistream is left in the badbit state\n";
183  return is;
184  }
186  return is;
187 }
188 
189 } // namespace CLHEP
CLHEP::RandGaussZiggurat::ziggurat_init
static bool ziggurat_init()
Definition: RandGaussZiggurat.cc:24
CLHEP::RandGaussZiggurat::get
std::istream & get(std::istream &is)
Definition: RandGaussZiggurat.cc:174
CLHEP::RandGauss::engine
HepRandomEngine & engine()
Definition: RandGauss.cc:44
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandGaussZiggurat::fn
static float fn[128]
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:109
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::RandGaussZiggurat::fe
static float fe[256]
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:109
CLHEP::RandGaussZiggurat::fire
float fire()
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:67
CLHEP::RandGaussZiggurat::ziggurat_RNOR
static float ziggurat_RNOR(HepRandomEngine *anEngine)
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:115
CLHEP::RandGauss::get
std::istream & get(std::istream &is)
Definition: RandGauss.cc:262
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandGaussZiggurat::shootArray
static void shootArray(const int size, float *vect, float mean=0.0, float stdDev=1.0)
Definition: RandGaussZiggurat.cc:110
CLHEP::RandGauss::defaultStdDev
double defaultStdDev
Definition: Matrix/CLHEP/Random/RandGauss.h:153
CLHEP::RandGaussZiggurat::operator()
virtual double operator()()
Definition: RandGaussZiggurat.cc:102
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandGauss::defaultMean
double defaultMean
Definition: Matrix/CLHEP/Random/RandGauss.h:152
CLHEP::RandGaussZiggurat::ziggurat_nfix
static float ziggurat_nfix(long hz, HepRandomEngine *anEngine)
Definition: RandGaussZiggurat.cc:75
CLHEP::RandGaussZiggurat::engine
HepRandomEngine & engine()
Definition: RandGaussZiggurat.cc:14
CLHEP::RandGaussZiggurat::fireArray
void fireArray(const int size, float *vect)
Definition: RandGaussZiggurat.cc:138
CLHEP::RandGaussZiggurat::wn
static float wn[128]
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:109
CLHEP::RandGauss::put
std::ostream & put(std::ostream &os) const
Definition: RandGauss.cc:237
CLHEP::RandGaussZiggurat::ziggurat_is_init
static bool ziggurat_is_init
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:111
CLHEP::RandGaussZiggurat::name
std::string name() const
Definition: RandGaussZiggurat.cc:19
CLHEP::RandGaussZiggurat::put
std::ostream & put(std::ostream &os) const
Definition: RandGaussZiggurat.cc:166
CLHEP::RandGaussZiggurat::kn
static unsigned long kn[128]
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:108
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::RandGaussZiggurat::~RandGaussZiggurat
virtual ~RandGaussZiggurat()
Definition: RandGaussZiggurat.cc:16
CLHEP::RandGaussZiggurat::we
static float we[256]
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:109
CLHEP::RandGauss::localEngine
shared_ptr< HepRandomEngine > localEngine
Definition: Matrix/CLHEP/Random/RandGauss.h:155
CLHEP::RandGaussZiggurat::ziggurat_UNI
static float ziggurat_UNI(HepRandomEngine *anEngine)
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:114
CLHEP::RandGaussZiggurat::ke
static unsigned long ke[256]
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:108
CLHEP::RandGaussZiggurat::shoot
static float shoot()
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:49
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
CLHEP::RandGaussZiggurat::ziggurat_SHR3
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
Definition: Matrix/CLHEP/Random/RandGaussZiggurat.h:113