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
CLHEP::RandExpZiggurat::engine
HepRandomEngine & engine()
Definition: RandExpZiggurat.cc:17
CLHEP::RandExpZiggurat
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:35
CLHEP::RandExpZiggurat::ziggurat_SHR3
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:134
CLHEP::RandExpZiggurat::put
std::ostream & put(std::ostream &os) const
Definition: RandExpZiggurat.cc:73
CLHEP::RandExpZiggurat::shoot
static float shoot()
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:55
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandExpZiggurat::~RandExpZiggurat
virtual ~RandExpZiggurat()
Definition: RandExpZiggurat.cc:19
CLHEP::RandExpZiggurat::ziggurat_init
static bool ziggurat_init()
Definition: RandExpZiggurat.cc:130
CLHEP::RandExpZiggurat::name
std::string name() const
Definition: RandExpZiggurat.cc:15
CLHEP::RandExpZiggurat::kn
static unsigned long kn[128]
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:129
CLHEP::RandExpZiggurat::shootArray
static void shootArray(const int size, float *vect, float mean=1.0)
Definition: RandExpZiggurat.cc:33
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::RandExpZiggurat::ziggurat_efix
static float ziggurat_efix(unsigned long jz, HepRandomEngine *anEngine)
Definition: RandExpZiggurat.cc:112
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandExpZiggurat::fe
static float fe[256]
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:130
CLHEP::RandExpZiggurat::fn
static float fn[128]
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:130
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandExpZiggurat::fire
float fire()
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:84
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RandExpZiggurat::ziggurat_UNI
static float ziggurat_UNI(HepRandomEngine *anEngine)
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:135
CLHEP::RandExpZiggurat::ziggurat_is_init
static bool ziggurat_is_init
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:132
CLHEP::RandExpZiggurat::operator()
virtual double operator()()
Definition: RandExpZiggurat.cc:28
CLHEP::RandExpZiggurat::we
static float we[256]
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:130
CLHEP::RandExpZiggurat::RandExpZiggurat
RandExpZiggurat(HepRandomEngine &anEngine, double mean=1.0)
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:162
CLHEP::RandExpZiggurat::get
std::istream & get(std::istream &is)
Definition: RandExpZiggurat.cc:91
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::RandExpZiggurat::wn
static float wn[128]
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:130
CLHEP::RandExpZiggurat::fireArray
void fireArray(const int size, float *vect)
Definition: RandExpZiggurat.cc:53
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
CLHEP::HepRandom
Definition: Matrix/CLHEP/Random/Random.h:50
CLHEP::RandExpZiggurat::ke
static unsigned long ke[256]
Definition: Matrix/CLHEP/Random/RandExpZiggurat.h:129