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

RandGamma.cc
Go to the documentation of this file.
1// $Id: RandGamma.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2// -*- C++ -*-
3//
4// -----------------------------------------------------------------------
5// HEP Random
6// --- RandGamma ---
7// class implementation file
8// -----------------------------------------------------------------------
9
10// =======================================================================
11// John Marraffino - Created: 12th May 1998
12// M Fischler - put and get to/from streams 12/13/04
13// M Fischler - put/get to/from streams uses pairs of ulongs when
14// + storing doubles avoid problems with precision
15// 4/14/05
16// =======================================================================
17
18#include "CLHEP/Random/defs.h"
19#include "CLHEP/Random/RandGamma.h"
20#include "CLHEP/Random/DoubConv.hh"
21#include <cmath> // for std::log()
22
23namespace CLHEP {
24
25std::string RandGamma::name() const {return "RandGamma";}
26HepRandomEngine & RandGamma::engine() {return *localEngine;}
27
29}
30
31double RandGamma::shoot( HepRandomEngine *anEngine, double k,
32 double lambda ) {
33 return genGamma( anEngine, k, lambda );
34}
35
36double RandGamma::shoot( double k, double lambda ) {
38 return genGamma( anEngine, k, lambda );
39}
40
41double RandGamma::fire( double k, double lambda ) {
42 return genGamma( localEngine.get(), k, lambda );
43}
44
45void RandGamma::shootArray( const int size, double* vect,
46 double k, double lambda )
47{
48 for( double* v = vect; v != vect + size; ++v )
49 *v = shoot(k,lambda);
50}
51
53 const int size, double* vect,
54 double k, double lambda )
55{
56 for( double* v = vect; v != vect + size; ++v )
57 *v = shoot(anEngine,k,lambda);
58}
59
60void RandGamma::fireArray( const int size, double* vect)
61{
62 for( double* v = vect; v != vect + size; ++v )
63 *v = fire(defaultK,defaultLambda);
64}
65
66void RandGamma::fireArray( const int size, double* vect,
67 double k, double lambda )
68{
69 for( double* v = vect; v != vect + size; ++v )
70 *v = fire(k,lambda);
71}
72
73double RandGamma::genGamma( HepRandomEngine *anEngine,
74 double a, double lambda ) {
75/*************************************************************************
76 * Gamma Distribution - Rejection algorithm gs combined with *
77 * Acceptance complement method gd *
78 *************************************************************************/
79
80static double aa = -1.0, aaa = -1.0, b, c, d, e, r, s, si, ss, q0,
81 q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875,
82 q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332,
83 q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320,
84 a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867,
85 a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581,
86 a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866,
87 e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848,
88 e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826,
89 e7 = 0.000247453;
90
91double gds,p,q,t,sign_u,u,v,w,x;
92double v1,v2,v12;
93
94// Check for invalid input values
95
96 if( a <= 0.0 ) return (-1.0);
97 if( lambda <= 0.0 ) return (-1.0);
98
99 if (a < 1.0)
100 { // CASE A: Acceptance rejection algorithm gs
101 b = 1.0 + 0.36788794412 * a; // Step 1
102 for(;;)
103 {
104 p = b * anEngine->flat();
105 if (p <= 1.0)
106 { // Step 2. Case gds <= 1
107 gds = std::exp(std::log(p) / a);
108 if (std::log(anEngine->flat()) <= -gds) return(gds/lambda);
109 }
110 else
111 { // Step 3. Case gds > 1
112 gds = - std::log ((b - p) / a);
113 if (std::log(anEngine->flat()) <= ((a - 1.0) * std::log(gds))) return(gds/lambda);
114 }
115 }
116 }
117 else
118 { // CASE B: Acceptance complement algorithm gd
119 if (a != aa)
120 { // Step 1. Preparations
121 aa = a;
122 ss = a - 0.5;
123 s = std::sqrt(ss);
124 d = 5.656854249 - 12.0 * s;
125 }
126 // Step 2. Normal deviate
127 do {
128 v1 = 2.0 * anEngine->flat() - 1.0;
129 v2 = 2.0 * anEngine->flat() - 1.0;
130 v12 = v1*v1 + v2*v2;
131 } while ( v12 > 1.0 );
132 t = v1*std::sqrt(-2.0*std::log(v12)/v12);
133 x = s + 0.5 * t;
134 gds = x * x;
135 if (t >= 0.0) return(gds/lambda); // Immediate acceptance
136
137 u = anEngine->flat(); // Step 3. Uniform random number
138 if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance
139
140 if (a != aaa)
141 { // Step 4. Set-up for hat case
142 aaa = a;
143 r = 1.0 / a;
144 q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
145 r + q3) * r + q2) * r + q1) * r;
146 if (a > 3.686)
147 {
148 if (a > 13.022)
149 {
150 b = 1.77;
151 si = 0.75;
152 c = 0.1515 / s;
153 }
154 else
155 {
156 b = 1.654 + 0.0076 * ss;
157 si = 1.68 / s + 0.275;
158 c = 0.062 / s + 0.024;
159 }
160 }
161 else
162 {
163 b = 0.463 + s - 0.178 * ss;
164 si = 1.235;
165 c = 0.195 / s - 0.079 + 0.016 * s;
166 }
167 }
168 if (x > 0.0) // Step 5. Calculation of q
169 {
170 v = t / (s + s); // Step 6.
171 if (std::fabs(v) > 0.25)
172 {
173 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v);
174 }
175 else
176 {
177 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
178 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
179 } // Step 7. Quotient acceptance
180 if (std::log(1.0 - u) <= q) return(gds/lambda);
181 }
182
183 for(;;)
184 { // Step 8. Double exponential deviate t
185 do
186 {
187 e = -std::log(anEngine->flat());
188 u = anEngine->flat();
189 u = u + u - 1.0;
190 sign_u = (u > 0)? 1.0 : -1.0;
191 t = b + (e * si) * sign_u;
192 }
193 while (t <= -0.71874483771719); // Step 9. Rejection of t
194 v = t / (s + s); // Step 10. New q(t)
195 if (std::fabs(v) > 0.25)
196 {
197 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v);
198 }
199 else
200 {
201 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
202 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
203 }
204 if (q <= 0.0) continue; // Step 11.
205 if (q > 0.5)
206 {
207 w = std::exp(q) - 1.0;
208 }
209 else
210 {
211 w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
212 q + e1) * q;
213 } // Step 12. Hat acceptance
214 if ( c * u * sign_u <= w * std::exp(e - 0.5 * t * t))
215 {
216 x = s + 0.5 * t;
217 return(x*x/lambda);
218 }
219 }
220 }
221}
222
223std::ostream & RandGamma::put ( std::ostream & os ) const {
224 int pr=os.precision(20);
225 std::vector<unsigned long> t(2);
226 os << " " << name() << "\n";
227 os << "Uvec" << "\n";
228 t = DoubConv::dto2longs(defaultK);
229 os << defaultK << " " << t[0] << " " << t[1] << "\n";
230 t = DoubConv::dto2longs(defaultLambda);
231 os << defaultLambda << " " << t[0] << " " << t[1] << "\n";
232 os.precision(pr);
233 return os;
234#ifdef REMOVED
235 int pr=os.precision(20);
236 os << " " << name() << "\n";
237 os << defaultK << " " << defaultLambda << "\n";
238 os.precision(pr);
239 return os;
240#endif
241}
242
243std::istream & RandGamma::get ( std::istream & is ) {
244 std::string inName;
245 is >> inName;
246 if (inName != name()) {
247 is.clear(std::ios::badbit | is.rdstate());
248 std::cerr << "Mismatch when expecting to read state of a "
249 << name() << " distribution\n"
250 << "Name found was " << inName
251 << "\nistream is left in the badbit state\n";
252 return is;
253 }
254 if (possibleKeywordInput(is, "Uvec", defaultK)) {
255 std::vector<unsigned long> t(2);
256 is >> defaultK >> t[0] >> t[1]; defaultK = DoubConv::longs2double(t);
257 is >> defaultLambda>>t[0]>>t[1]; defaultLambda = DoubConv::longs2double(t);
258 return is;
259 }
260 // is >> defaultK encompassed by possibleKeywordInput
261 is >> defaultLambda;
262 return is;
263}
264
265} // namespace CLHEP
266
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
virtual double flat()=0
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
std::string name() const
Definition: RandGamma.cc:25
HepRandomEngine & engine()
Definition: RandGamma.cc:26
static void shootArray(const int size, double *vect, double k=1.0, double lambda=1.0)
Definition: RandGamma.cc:45
std::ostream & put(std::ostream &os) const
Definition: RandGamma.cc:223
virtual ~RandGamma()
Definition: RandGamma.cc:28
static double shoot()
void fireArray(const int size, double *vect)
Definition: RandGamma.cc:60
std::istream & get(std::istream &is)
Definition: RandGamma.cc:243
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
@ b
@ a