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
9 CLHEP_THREAD_LOCAL unsigned long RandGaussZiggurat::kn[128], RandGaussZiggurat::ke[256];
10 CLHEP_THREAD_LOCAL float RandGaussZiggurat::wn[128],RandGaussZiggurat::fn[128],RandGaussZiggurat::we[256],RandGaussZiggurat::fe[256];
11 CLHEP_THREAD_LOCAL bool RandGaussZiggurat::ziggurat_is_init = false;
12
engine()13 HepRandomEngine & RandGaussZiggurat::engine() {return RandGauss::engine();}
14
~RandGaussZiggurat()15 RandGaussZiggurat::~RandGaussZiggurat() {
16 }
17
name() const18 std::string RandGaussZiggurat::name() const
19 {
20 return "RandGaussZiggurat";
21 }
22
ziggurat_init()23 bool RandGaussZiggurat::ziggurat_init()
24 {
25 const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
26 double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
27 double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
28 int i;
29
30 /* Set up tables for RNOR */
31 q=vn/exp(-.5*dn*dn);
32 kn[0]=(unsigned long)((dn/q)*rzm1);
33 kn[1]=0;
34
35 wn[0]=q/rzm1;
36 wn[127]=dn/rzm1;
37
38 fn[0]=1.;
39 fn[127]=exp(-.5*dn*dn);
40
41 for(i=126;i>=1;i--) {
42 dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
43 kn[i+1]=(unsigned long)((dn/tn)*rzm1);
44 tn=dn;
45 fn[i]=exp(-.5*dn*dn);
46 wn[i]=dn/rzm1;
47 }
48
49 /* Set up tables for REXP */
50 q = ve/exp(-de);
51 ke[0]=(unsigned long)((de/q)*rzm2);
52 ke[1]=0;
53
54 we[0]=q/rzm2;
55 we[255]=de/rzm2;
56
57 fe[0]=1.;
58 fe[255]=exp(-de);
59
60 for(i=254;i>=1;i--) {
61 de=-log(ve/de+exp(-de));
62 ke[i+1]= (unsigned long)((de/te)*rzm2);
63 te=de;
64 fe[i]=exp(-de);
65 we[i]=de/rzm2;
66 }
67 ziggurat_is_init=true;
68
69 //std::cout<<"Done RandGaussZiggurat::ziggurat_init()"<<std::endl;
70
71 return true;
72 }
73
ziggurat_nfix(long hz,HepRandomEngine * anEngine)74 float RandGaussZiggurat::ziggurat_nfix(long hz,HepRandomEngine* anEngine)
75 {
76 if(!ziggurat_is_init) ziggurat_init();
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)std::abs(hz)<kn[iz]) return (hz*wn[iz]);
99 }
100 }
101
operator ()()102 double RandGaussZiggurat::operator()() {
103 return ziggurat_RNOR(localEngine.get()) * defaultStdDev + defaultMean;
104 }
105
operator ()(double mean,double stdDev)106 double RandGaussZiggurat::operator()( double mean, double stdDev ) {
107 return ziggurat_RNOR(localEngine.get()) * stdDev + mean;
108 }
109
shootArray(const int size,float * vect,float mean,float stdDev)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
shootArray(const int size,double * vect,double mean,double stdDev)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
shootArray(HepRandomEngine * anEngine,const int size,float * vect,float mean,float stdDev)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
shootArray(HepRandomEngine * anEngine,const int size,double * vect,double mean,double stdDev)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
fireArray(const int size,float * vect)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
fireArray(const int size,double * vect)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
fireArray(const int size,float * vect,float mean,float stdDev)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
fireArray(const int size,double * vect,double mean,double stdDev)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
put(std::ostream & os) const166 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
get(std::istream & is)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
190