19 #ifndef BZ_RANDOM_GAMMA
20 #define BZ_RANDOM_GAMMA
22 #ifndef BZ_RANDOM_UNIFORM
26 #ifndef BZ_RANDOM_NORMAL
30 #ifndef BZ_RANDOM_EXPONENTIAL
34 #ifndef BZ_NUMINQUIRE_H
35 #include <blitz/numinquire.h>
40 template<typename T =
double, typename IRNG =
defaultIRNG,
61 BZPRECONDITION(mean >= 1.0);
73 return normal_.random();
78 return exponential_.random();
81 T fsign(T num, T sign)
85 if ((sign>0.0L && num<0.0L)||(sign<0.0L && num>0.0L))
97 template<
typename T,
typename IRNG,
typename stateTag>
110 static T q1 = 4.166669E-2;
111 static T q2 = 2.083148E-2;
112 static T q3 = 8.01191E-3;
113 static T q4 = 1.44121E-3;
114 static T q5 = -7.388E-5;
115 static T q6 = 2.4511E-4;
116 static T q7 = 2.424E-4;
117 static T a1 = 0.3333333;
118 static T a2 = -0.250003;
119 static T a3 = 0.2000062;
120 static T a4 = -0.1662921;
121 static T a5 = 0.1423657;
122 static T a6 = -0.1367177;
123 static T a7 = 0.1233795;
125 static T e2 = 0.4999897;
126 static T e3 = 0.166829;
127 static T e4 = 4.07753E-2;
128 static T e5 = 1.0293E-2;
131 static T sqrt32 = 5.656854249492380195206754896838792314280;
134 static T sgamma,s2,
s,d,
t,x,u,
r,q0,b,b0,si,c,v,
q,e,w,
p;
136 if(
a == aa)
goto S10;
137 if(
a < 1.0)
goto S120;
154 if(t >= 0.0)
return sgamma;
159 if(d*u <= t*t*t)
return sgamma;
163 if(
a == aaa)
goto S40;
166 q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*
r;
172 if(
a <= 3.686)
goto S30;
173 if(
a <= 13.022)
goto S20;
193 b = 0.463+s+0.178*s2;
195 c = 0.195/s-7.9E-2+1.6E-1*
s;
200 if(x <= 0.0)
goto S70;
205 if(fabs(v) <= 0.25)
goto S50;
206 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
209 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
214 if(log(1.0-u) <=
q)
return sgamma;
228 if(t < -0.7187449)
goto S70;
233 if(fabs(v) <= 0.25)
goto S80;
234 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
237 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
242 if(q <= 0.0)
goto S70;
243 if(q <= 0.5)
goto S100;
247 if(q < 15.0)
goto S95;
254 if((q+e-0.5*t*t) > 87.49823)
goto S115;
255 if(c*fabs(u) > exp(q+e-0.5*t*t))
goto S70;
261 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*
q;
266 if(c*fabs(u) > w*exp(e-0.5*t*t))
goto S70;
285 b0 = 1.0+0.3678794*
a;
288 if(p >= 1.0)
goto S140;
289 sgamma = exp(log(p)/
a);
290 if(sexpo() < sgamma)
goto S130;
293 sgamma = -log((b0-p)/
a);
294 if(sexpo() < (1.0-
a)*log(sgamma))
goto S130;
301 #endif // BZ_RANDOM_GAMMA