00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef BZ_RANDOM_GAMMA
00020 #define BZ_RANDOM_GAMMA
00021
00022 #ifndef BZ_RANDOM_UNIFORM
00023 #include <random/uniform.h>
00024 #endif
00025
00026 #ifndef BZ_RANDOM_NORMAL
00027 #include <random/normal.h>
00028 #endif
00029
00030 #ifndef BZ_RANDOM_EXPONENTIAL
00031 #include <random/exponential.h>
00032 #endif
00033
00034 #ifndef BZ_NUMINQUIRE_H
00035 #include <blitz/numinquire.h>
00036 #endif
00037
00038 BZ_NAMESPACE(ranlib)
00039
00040 template<typename T = double, typename IRNG = defaultIRNG,
00041 typename stateTag = defaultState>
00042 class Gamma : public UniformOpen<T,IRNG,stateTag>
00043 {
00044 public:
00045 typedef T T_numtype;
00046
00047 Gamma()
00048 {
00049 setMean(1.0);
00050 }
00051
00052 explicit Gamma(unsigned int i) :
00053 UniformOpen<T,IRNG,stateTag>(i)
00054 {
00055 setMean(1.0);
00056 };
00057
00058 Gamma(T mean)
00059 {
00060 setMean(mean);
00061 }
00062
00063 Gamma(T mean, unsigned int i) :
00064 UniformOpen<T,IRNG,stateTag>(i)
00065 {
00066 setMean(mean);
00067 };
00068
00069 T random();
00070
00071 void setMean(T mean)
00072 {
00073 BZPRECONDITION(mean >= 1.0);
00074 a = mean;
00075 }
00076
00077 protected:
00078 T ranf()
00079 {
00080 return UniformOpen<T,IRNG,stateTag>::random();
00081 }
00082
00083 T snorm()
00084 {
00085 return normal_.random();
00086 }
00087
00088 T sexpo()
00089 {
00090 return exponential_.random();
00091 }
00092
00093 T fsign(T num, T sign)
00094 {
00095
00096
00097 if ((sign>0.0L && num<0.0L)||(sign<0.0L && num>0.0L))
00098 return -num;
00099 else
00100 return num;
00101 }
00102
00103 NormalUnit<T,IRNG,sharedState> normal_;
00104 ExponentialUnit<T,IRNG,sharedState> exponential_;
00105
00106 T a;
00107 };
00108
00109 template<typename T, typename IRNG, typename stateTag>
00110 T Gamma<T,IRNG,stateTag>::random()
00111 {
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 static T q1 = 4.166669E-2;
00123 static T q2 = 2.083148E-2;
00124 static T q3 = 8.01191E-3;
00125 static T q4 = 1.44121E-3;
00126 static T q5 = -7.388E-5;
00127 static T q6 = 2.4511E-4;
00128 static T q7 = 2.424E-4;
00129 static T a1 = 0.3333333;
00130 static T a2 = -0.250003;
00131 static T a3 = 0.2000062;
00132 static T a4 = -0.1662921;
00133 static T a5 = 0.1423657;
00134 static T a6 = -0.1367177;
00135 static T a7 = 0.1233795;
00136 static T e1 = 1.0;
00137 static T e2 = 0.4999897;
00138 static T e3 = 0.166829;
00139 static T e4 = 4.07753E-2;
00140 static T e5 = 1.0293E-2;
00141 static T aa = 0.0;
00142 static T aaa = 0.0;
00143 static T sqrt32 = 5.656854249492380195206754896838792314280;
00144
00145
00146 static T sgamma,s2,s,d,t,x,u,r,q0,b,b0,si,c,v,q,e,w,p;
00147
00148 if(a == aa) goto S10;
00149 if(a < 1.0) goto S120;
00150
00151
00152
00153 aa = a;
00154 s2 = a-0.5;
00155 s = sqrt(s2);
00156 d = sqrt32-12.0*s;
00157 S10:
00158
00159
00160
00161
00162
00163 t = snorm();
00164 x = s+0.5*t;
00165 sgamma = x*x;
00166 if(t >= 0.0) return sgamma;
00167
00168
00169
00170 u = ranf();
00171 if(d*u <= t*t*t) return sgamma;
00172
00173
00174
00175 if(a == aaa) goto S40;
00176 aaa = a;
00177 r = 1.0/ a;
00178 q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
00179
00180
00181
00182
00183
00184 if(a <= 3.686) goto S30;
00185 if(a <= 13.022) goto S20;
00186
00187
00188
00189 b = 1.77;
00190 si = 0.75;
00191 c = 0.1515/s;
00192 goto S40;
00193 S20:
00194
00195
00196
00197 b = 1.654+7.6E-3*s2;
00198 si = 1.68/s+0.275;
00199 c = 6.2E-2/s+2.4E-2;
00200 goto S40;
00201 S30:
00202
00203
00204
00205 b = 0.463+s+0.178*s2;
00206 si = 1.235;
00207 c = 0.195/s-7.9E-2+1.6E-1*s;
00208 S40:
00209
00210
00211
00212 if(x <= 0.0) goto S70;
00213
00214
00215
00216 v = t/(s+s);
00217 if(fabs(v) <= 0.25) goto S50;
00218 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00219 goto S60;
00220 S50:
00221 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00222 S60:
00223
00224
00225
00226 if(log(1.0-u) <= q) return sgamma;
00227 S70:
00228
00229
00230
00231
00232
00233 e = sexpo();
00234 u = ranf();
00235 u += (u-1.0);
00236 t = b+fsign(si*e,u);
00237
00238
00239
00240 if(t < -0.7187449) goto S70;
00241
00242
00243
00244 v = t/(s+s);
00245 if(fabs(v) <= 0.25) goto S80;
00246 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00247 goto S90;
00248 S80:
00249 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00250 S90:
00251
00252
00253
00254 if(q <= 0.0) goto S70;
00255 if(q <= 0.5) goto S100;
00256
00257
00258
00259 if(q < 15.0) goto S95;
00260
00261
00262
00263
00264
00265
00266 if((q+e-0.5*t*t) > 87.49823) goto S115;
00267 if(c*fabs(u) > exp(q+e-0.5*t*t)) goto S70;
00268 goto S115;
00269 S95:
00270 w = exp(q)-1.0;
00271 goto S110;
00272 S100:
00273 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
00274 S110:
00275
00276
00277
00278 if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
00279 S115:
00280 x = s+0.5*t;
00281 sgamma = x*x;
00282 return sgamma;
00283 S120:
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 b0 = 1.0+0.3678794*a;
00298 S130:
00299 p = b0*ranf();
00300 if(p >= 1.0) goto S140;
00301 sgamma = exp(log(p)/ a);
00302 if(sexpo() < sgamma) goto S130;
00303 return sgamma;
00304 S140:
00305 sgamma = -log((b0-p)/ a);
00306 if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
00307 return sgamma;
00308
00309 }
00310
00311 BZ_NAMESPACE_END
00312
00313 #endif // BZ_RANDOM_GAMMA