1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandExpZiggurat ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9
10 // =======================================================================
11
12 #include "CLHEP/Random/defs.h"
13 #include "CLHEP/Random/DoubConv.hh"
14
15 #include "CLHEP/Random/RandExpZiggurat.h"
16
17 #include <iostream>
18 #include <cmath> // for std::log()
19
20 namespace CLHEP {
21
22 CLHEP_THREAD_LOCAL unsigned long RandExpZiggurat::kn[128], RandExpZiggurat::ke[256];
23 CLHEP_THREAD_LOCAL float RandExpZiggurat::wn[128],RandExpZiggurat::fn[128],RandExpZiggurat::we[256],RandExpZiggurat::fe[256];
24 CLHEP_THREAD_LOCAL bool RandExpZiggurat::ziggurat_is_init = false;
25
name() const26 std::string RandExpZiggurat::name() const {return "RandExpZiggurat";}
27
engine()28 HepRandomEngine & RandExpZiggurat::engine() {return *localEngine;}
29
~RandExpZiggurat()30 RandExpZiggurat::~RandExpZiggurat() {
31 }
32
RandExpZiggurat(const RandExpZiggurat & right)33 RandExpZiggurat::RandExpZiggurat(const RandExpZiggurat& right) : HepRandom(right),defaultMean(right.defaultMean)
34 {
35 }
36
operator ()()37 double RandExpZiggurat::operator()()
38 {
39 return fire( defaultMean );
40 }
41
shootArray(const int size,float * vect,float mean)42 void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
43 {
44 for (int i=0; i<size; ++i) vect[i] = shoot(mean);
45 }
46
shootArray(const int size,double * vect,double mean)47 void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
48 {
49 for (int i=0; i<size; ++i) vect[i] = shoot(mean);
50 }
51
shootArray(HepRandomEngine * anEngine,const int size,float * vect,float mean)52 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, float* vect, float mean )
53 {
54 for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
55 }
56
shootArray(HepRandomEngine * anEngine,const int size,double * vect,double mean)57 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, double* vect, double mean )
58 {
59 for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
60 }
61
fireArray(const int size,float * vect)62 void RandExpZiggurat::fireArray( const int size, float* vect)
63 {
64 for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
65 }
66
fireArray(const int size,double * vect)67 void RandExpZiggurat::fireArray( const int size, double* vect)
68 {
69 for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
70 }
71
fireArray(const int size,float * vect,float mean)72 void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
73 {
74 for (int i=0; i<size; ++i) vect[i] = fire( mean );
75 }
76
fireArray(const int size,double * vect,double mean)77 void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
78 {
79 for (int i=0; i<size; ++i) vect[i] = fire( mean );
80 }
81
put(std::ostream & os) const82 std::ostream & RandExpZiggurat::put ( std::ostream & os ) const {
83 int pr=os.precision(20);
84 std::vector<unsigned long> t(2);
85 os << " " << name() << "\n";
86 os << "Uvec" << "\n";
87 t = DoubConv::dto2longs(defaultMean);
88 os << defaultMean << " " << t[0] << " " << t[1] << "\n";
89 os.precision(pr);
90 return os;
91 #ifdef REMOVED
92 int pr=os.precision(20);
93 os << " " << name() << "\n";
94 os << defaultMean << "\n";
95 os.precision(pr);
96 return os;
97 #endif
98 }
99
get(std::istream & is)100 std::istream & RandExpZiggurat::get ( std::istream & is ) {
101 std::string inName;
102 is >> inName;
103 if (inName != name()) {
104 is.clear(std::ios::badbit | is.rdstate());
105 std::cerr << "Mismatch when expecting to read state of a "
106 << name() << " distribution\n"
107 << "Name found was " << inName
108 << "\nistream is left in the badbit state\n";
109 return is;
110 }
111 if (possibleKeywordInput(is, "Uvec", defaultMean)) {
112 std::vector<unsigned long> t(2);
113 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
114 return is;
115 }
116 // is >> defaultMean encompassed by possibleKeywordInput
117 return is;
118 }
119
120
ziggurat_efix(unsigned long jz,HepRandomEngine * anEngine)121 float RandExpZiggurat::ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine)
122 {
123 if(!ziggurat_is_init) ziggurat_init();
124
125 unsigned long iz=jz&255;
126
127 float x;
128 for(;;)
129 {
130 if(iz==0) return (7.69711-std::log(ziggurat_UNI(anEngine))); /* iz==0 */
131 x=jz*we[iz];
132 if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1]-fe[iz]) < std::exp(-x) ) return (x);
133
134 /* initiate, try to exit for(;;) loop */
135 jz=ziggurat_SHR3(anEngine);
136 iz=(jz&255);
137 if(jz<ke[iz]) return (jz*we[iz]);
138 }
139 }
140
ziggurat_init()141 bool RandExpZiggurat::ziggurat_init()
142 {
143 const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
144 double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
145 double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
146 int i;
147
148 /* Set up tables for RNOR */
149 q=vn/std::exp(-.5*dn*dn);
150 kn[0]=(unsigned long)((dn/q)*rzm1);
151 kn[1]=0;
152
153 wn[0]=q/rzm1;
154 wn[127]=dn/rzm1;
155
156 fn[0]=1.;
157 fn[127]=std::exp(-.5*dn*dn);
158
159 for(i=126;i>=1;i--) {
160 dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
161 kn[i+1]=(unsigned long)((dn/tn)*rzm1);
162 tn=dn;
163 fn[i]=std::exp(-.5*dn*dn);
164 wn[i]=dn/rzm1;
165 }
166
167 /* Set up tables for REXP */
168 q = ve/std::exp(-de);
169 ke[0]=(unsigned long)((de/q)*rzm2);
170 ke[1]=0;
171
172 we[0]=q/rzm2;
173 we[255]=de/rzm2;
174
175 fe[0]=1.;
176 fe[255]=std::exp(-de);
177
178 for(i=254;i>=1;i--) {
179 de=-std::log(ve/de+std::exp(-de));
180 ke[i+1]= (unsigned long)((de/te)*rzm2);
181 te=de;
182 fe[i]=std::exp(-de);
183 we[i]=de/rzm2;
184 }
185 ziggurat_is_init=true;
186 return true;
187 }
188
189 } // namespace CLHEP
190