CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandGamma.cc,v 1.6 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandGamma --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // John Marraffino - Created: 12th May 1998 00012 // M Fischler - put and get to/from streams 12/13/04 00013 // M Fischler - put/get to/from streams uses pairs of ulongs when 00014 // + storing doubles avoid problems with precision 00015 // 4/14/05 00016 // ======================================================================= 00017 00018 #include "CLHEP/Random/defs.h" 00019 #include "CLHEP/Random/RandGamma.h" 00020 #include "CLHEP/Random/DoubConv.hh" 00021 #include <cmath> // for std::log() 00022 00023 namespace CLHEP { 00024 00025 std::string RandGamma::name() const {return "RandGamma";} 00026 HepRandomEngine & RandGamma::engine() {return *localEngine;} 00027 00028 RandGamma::~RandGamma() { 00029 } 00030 00031 double RandGamma::shoot( HepRandomEngine *anEngine, double k, 00032 double lambda ) { 00033 return genGamma( anEngine, k, lambda ); 00034 } 00035 00036 double RandGamma::shoot( double k, double lambda ) { 00037 HepRandomEngine *anEngine = HepRandom::getTheEngine(); 00038 return genGamma( anEngine, k, lambda ); 00039 } 00040 00041 double RandGamma::fire( double k, double lambda ) { 00042 return genGamma( localEngine.get(), k, lambda ); 00043 } 00044 00045 void RandGamma::shootArray( const int size, double* vect, 00046 double k, double lambda ) 00047 { 00048 for( double* v = vect; v != vect + size; ++v ) 00049 *v = shoot(k,lambda); 00050 } 00051 00052 void RandGamma::shootArray( HepRandomEngine* anEngine, 00053 const int size, double* vect, 00054 double k, double lambda ) 00055 { 00056 for( double* v = vect; v != vect + size; ++v ) 00057 *v = shoot(anEngine,k,lambda); 00058 } 00059 00060 void RandGamma::fireArray( const int size, double* vect) 00061 { 00062 for( double* v = vect; v != vect + size; ++v ) 00063 *v = fire(defaultK,defaultLambda); 00064 } 00065 00066 void RandGamma::fireArray( const int size, double* vect, 00067 double k, double lambda ) 00068 { 00069 for( double* v = vect; v != vect + size; ++v ) 00070 *v = fire(k,lambda); 00071 } 00072 00073 double RandGamma::genGamma( HepRandomEngine *anEngine, 00074 double a, double lambda ) { 00075 /************************************************************************* 00076 * Gamma Distribution - Rejection algorithm gs combined with * 00077 * Acceptance complement method gd * 00078 *************************************************************************/ 00079 00080 static double aa = -1.0, aaa = -1.0, b, c, d, e, r, s, si, ss, q0, 00081 q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875, 00082 q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332, 00083 q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320, 00084 a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867, 00085 a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581, 00086 a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866, 00087 e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848, 00088 e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826, 00089 e7 = 0.000247453; 00090 00091 double gds,p,q,t,sign_u,u,v,w,x; 00092 double v1,v2,v12; 00093 00094 // Check for invalid input values 00095 00096 if( a <= 0.0 ) return (-1.0); 00097 if( lambda <= 0.0 ) return (-1.0); 00098 00099 if (a < 1.0) 00100 { // CASE A: Acceptance rejection algorithm gs 00101 b = 1.0 + 0.36788794412 * a; // Step 1 00102 for(;;) 00103 { 00104 p = b * anEngine->flat(); 00105 if (p <= 1.0) 00106 { // Step 2. Case gds <= 1 00107 gds = std::exp(std::log(p) / a); 00108 if (std::log(anEngine->flat()) <= -gds) return(gds/lambda); 00109 } 00110 else 00111 { // Step 3. Case gds > 1 00112 gds = - std::log ((b - p) / a); 00113 if (std::log(anEngine->flat()) <= ((a - 1.0) * std::log(gds))) return(gds/lambda); 00114 } 00115 } 00116 } 00117 else 00118 { // CASE B: Acceptance complement algorithm gd 00119 if (a != aa) 00120 { // Step 1. Preparations 00121 aa = a; 00122 ss = a - 0.5; 00123 s = std::sqrt(ss); 00124 d = 5.656854249 - 12.0 * s; 00125 } 00126 // Step 2. Normal deviate 00127 do { 00128 v1 = 2.0 * anEngine->flat() - 1.0; 00129 v2 = 2.0 * anEngine->flat() - 1.0; 00130 v12 = v1*v1 + v2*v2; 00131 } while ( v12 > 1.0 ); 00132 t = v1*std::sqrt(-2.0*std::log(v12)/v12); 00133 x = s + 0.5 * t; 00134 gds = x * x; 00135 if (t >= 0.0) return(gds/lambda); // Immediate acceptance 00136 00137 u = anEngine->flat(); // Step 3. Uniform random number 00138 if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance 00139 00140 if (a != aaa) 00141 { // Step 4. Set-up for hat case 00142 aaa = a; 00143 r = 1.0 / a; 00144 q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) * 00145 r + q3) * r + q2) * r + q1) * r; 00146 if (a > 3.686) 00147 { 00148 if (a > 13.022) 00149 { 00150 b = 1.77; 00151 si = 0.75; 00152 c = 0.1515 / s; 00153 } 00154 else 00155 { 00156 b = 1.654 + 0.0076 * ss; 00157 si = 1.68 / s + 0.275; 00158 c = 0.062 / s + 0.024; 00159 } 00160 } 00161 else 00162 { 00163 b = 0.463 + s - 0.178 * ss; 00164 si = 1.235; 00165 c = 0.195 / s - 0.079 + 0.016 * s; 00166 } 00167 } 00168 if (x > 0.0) // Step 5. Calculation of q 00169 { 00170 v = t / (s + s); // Step 6. 00171 if (std::fabs(v) > 0.25) 00172 { 00173 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v); 00174 } 00175 else 00176 { 00177 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * 00178 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v; 00179 } // Step 7. Quotient acceptance 00180 if (std::log(1.0 - u) <= q) return(gds/lambda); 00181 } 00182 00183 for(;;) 00184 { // Step 8. Double exponential deviate t 00185 do 00186 { 00187 e = -std::log(anEngine->flat()); 00188 u = anEngine->flat(); 00189 u = u + u - 1.0; 00190 sign_u = (u > 0)? 1.0 : -1.0; 00191 t = b + (e * si) * sign_u; 00192 } 00193 while (t <= -0.71874483771719); // Step 9. Rejection of t 00194 v = t / (s + s); // Step 10. New q(t) 00195 if (std::fabs(v) > 0.25) 00196 { 00197 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v); 00198 } 00199 else 00200 { 00201 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * 00202 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v; 00203 } 00204 if (q <= 0.0) continue; // Step 11. 00205 if (q > 0.5) 00206 { 00207 w = std::exp(q) - 1.0; 00208 } 00209 else 00210 { 00211 w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * 00212 q + e1) * q; 00213 } // Step 12. Hat acceptance 00214 if ( c * u * sign_u <= w * std::exp(e - 0.5 * t * t)) 00215 { 00216 x = s + 0.5 * t; 00217 return(x*x/lambda); 00218 } 00219 } 00220 } 00221 } 00222 00223 std::ostream & RandGamma::put ( std::ostream & os ) const { 00224 int pr=os.precision(20); 00225 std::vector<unsigned long> t(2); 00226 os << " " << name() << "\n"; 00227 os << "Uvec" << "\n"; 00228 t = DoubConv::dto2longs(defaultK); 00229 os << defaultK << " " << t[0] << " " << t[1] << "\n"; 00230 t = DoubConv::dto2longs(defaultLambda); 00231 os << defaultLambda << " " << t[0] << " " << t[1] << "\n"; 00232 os.precision(pr); 00233 return os; 00234 #ifdef REMOVED 00235 int pr=os.precision(20); 00236 os << " " << name() << "\n"; 00237 os << defaultK << " " << defaultLambda << "\n"; 00238 os.precision(pr); 00239 return os; 00240 #endif 00241 } 00242 00243 std::istream & RandGamma::get ( std::istream & is ) { 00244 std::string inName; 00245 is >> inName; 00246 if (inName != name()) { 00247 is.clear(std::ios::badbit | is.rdstate()); 00248 std::cerr << "Mismatch when expecting to read state of a " 00249 << name() << " distribution\n" 00250 << "Name found was " << inName 00251 << "\nistream is left in the badbit state\n"; 00252 return is; 00253 } 00254 if (possibleKeywordInput(is, "Uvec", defaultK)) { 00255 std::vector<unsigned long> t(2); 00256 is >> defaultK >> t[0] >> t[1]; defaultK = DoubConv::longs2double(t); 00257 is >> defaultLambda>>t[0]>>t[1]; defaultLambda = DoubConv::longs2double(t); 00258 return is; 00259 } 00260 // is >> defaultK encompassed by possibleKeywordInput 00261 is >> defaultLambda; 00262 return is; 00263 } 00264 00265 } // namespace CLHEP 00266