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
7namespace CLHEP {
8
10//bool RandGaussZiggurat::ziggurat_is_init=false;
11unsigned long RandGaussZiggurat::kn[128], RandGaussZiggurat::ke[256];
13
15
18
19std::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 }
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
105
106double RandGaussZiggurat::operator()( double mean, double stdDev ) {
107 return ziggurat_RNOR(localEngine.get()) * stdDev + mean;
108}
109
110void 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
117void 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
124void 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
131void 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
138void 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
145void 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
152void 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
159void 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
166std::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
174std::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::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