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