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  }
185  RandGauss::get(is);
186  return is;
187 }
188 
189 } // namespace CLHEP
static void shootArray(const int size, float *vect, float mean=0.0, float stdDev=1.0)
std::istream & get(std::istream &is)
void fireArray(const int size, float *vect)
static float ziggurat_nfix(long hz, HepRandomEngine *anEngine)
static float ziggurat_UNI(HepRandomEngine *anEngine)
std::ostream & put(std::ostream &os) const
static float ziggurat_RNOR(HepRandomEngine *anEngine)
HepRandomEngine & engine()
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
std::string name() const
std::istream & get(std::istream &is)
Definition: RandGauss.cc:262
std::ostream & put(std::ostream &os) const
Definition: RandGauss.cc:237
shared_ptr< HepRandomEngine > localEngine
HepRandomEngine & engine()
Definition: RandGauss.cc:44