CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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