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 
9 namespace CLHEP {
10 
12 unsigned long RandExpZiggurat::kn[128], RandExpZiggurat::ke[256];
14 
15 std::string RandExpZiggurat::name() const {return "RandExpZiggurat";}
16 
17 HepRandomEngine & RandExpZiggurat::engine() {return *localEngine;}
18 
20  if ( deleteEngine ) delete localEngine;
21 }
22 
23 RandExpZiggurat::RandExpZiggurat(const RandExpZiggurat& right) : HepRandom(right),defaultMean(right.defaultMean)
24 {
26 }
27 
29 {
30  return fire( defaultMean );
31 }
32 
33 void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
34 {
35  for (int i=0; i<size; ++i) vect[i] = shoot(mean);
36 }
37 
38 void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
39 {
40  for (int i=0; i<size; ++i) vect[i] = shoot(mean);
41 }
42 
43 void 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 
48 void 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 
53 void RandExpZiggurat::fireArray( const int size, float* vect)
54 {
55  for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
56 }
57 
58 void RandExpZiggurat::fireArray( const int size, double* vect)
59 {
60  for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
61 }
62 
63 void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
64 {
65  for (int i=0; i<size; ++i) vect[i] = fire( mean );
66 }
67 
68 void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
69 {
70  for (int i=0; i<size; ++i) vect[i] = fire( mean );
71 }
72 
73 std::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 
91 std::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 
112 float 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)