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