22 #ifndef BZ_RANDOM_GAMMA
23 #define BZ_RANDOM_GAMMA
25 #ifndef BZ_RANDOM_UNIFORM
29 #ifndef BZ_RANDOM_NORMAL
33 #ifndef BZ_RANDOM_EXPONENTIAL
37 #ifndef BZ_NUMINQUIRE_H
38 #include <blitz/numinquire.h>
43 template<
typename T = double,
typename IRNG =
defaultIRNG,
76 BZPRECONDITION(mean >= 1.0);
100 if ((sign>0.0L && num<0.0L)||(sign<0.0L && num>0.0L))
112 template<
typename T,
typename IRNG,
typename stateTag>
125 static T q1 = 4.166669E-2;
126 static T q2 = 2.083148E-2;
127 static T q3 = 8.01191E-3;
128 static T q4 = 1.44121E-3;
129 static T q5 = -7.388E-5;
130 static T q6 = 2.4511E-4;
131 static T q7 = 2.424E-4;
132 static T a1 = 0.3333333;
133 static T a2 = -0.250003;
134 static T a3 = 0.2000062;
135 static T a4 = -0.1662921;
136 static T a5 = 0.1423657;
137 static T a6 = -0.1367177;
138 static T a7 = 0.1233795;
140 static T e2 = 0.4999897;
141 static T e3 = 0.166829;
142 static T e4 = 4.07753E-2;
143 static T e5 = 1.0293E-2;
146 static T sqrt32 = 5.656854249492380195206754896838792314280;
149 static T sgamma,s2,
s,d,
t,x,u,
r,q0,
b,b0,si,c,v,
q,e,w,
p;
151 if(
a == aa)
goto S10;
152 if(
a < 1.0)
goto S120;
169 if(t >= 0.0)
return sgamma;
174 if(d*u <= t*t*t)
return sgamma;
178 if(
a == aaa)
goto S40;
181 q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*
r;
187 if(
a <= 3.686)
goto S30;
188 if(
a <= 13.022)
goto S20;
208 b = 0.463+s+0.178*s2;
210 c = 0.195/s-7.9E-2+1.6E-1*
s;
215 if(x <= 0.0)
goto S70;
220 if(fabs(v) <= 0.25)
goto S50;
221 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
224 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
229 if(log(1.0-u) <=
q)
return sgamma;
243 if(t < -0.7187449)
goto S70;
248 if(fabs(v) <= 0.25)
goto S80;
249 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
252 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
257 if(q <= 0.0)
goto S70;
258 if(q <= 0.5)
goto S100;
262 if(q < 15.0)
goto S95;
269 if((q+e-0.5*t*t) > 87.49823)
goto S115;
270 if(c*fabs(u) > exp(q+e-0.5*t*t))
goto S70;
276 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*
q;
281 if(c*fabs(u) > w*exp(e-0.5*t*t))
goto S70;
300 b0 = 1.0+0.3678794*
a;
303 if(p >= 1.0)
goto S140;
304 sgamma = exp(log(p)/
a);
305 if(sexpo() < sgamma)
goto S130;
308 sgamma = -log((b0-p)/
a);
309 if(sexpo() < (1.0-
a)*log(sgamma))
goto S130;
316 #endif // BZ_RANDOM_GAMMA
Definition: exponential.h:33
T fsign(T num, T sign)
Definition: gamma.h:96
_bz_global blitz::IndexPlaceholder< 9 > r
Definition: indexexpr.h:265
_bz_global blitz::IndexPlaceholder< 11 > t
Definition: indexexpr.h:267
_bz_global blitz::IndexPlaceholder< 0 > i
Definition: indexexpr.h:256
T ranf()
Definition: gamma.h:81
_bz_global blitz::IndexPlaceholder< 8 > q
Definition: indexexpr.h:264
T a
Definition: gamma.h:109
ExponentialUnit< T, IRNG, sharedState > exponential_
Definition: gamma.h:107
Gamma(unsigned int i)
Definition: gamma.h:55
T sexpo()
Definition: gamma.h:91
T T_numtype
Definition: gamma.h:48
Gamma()
Definition: gamma.h:50
MersenneTwister defaultIRNG
Definition: default.h:120
Gamma(T mean, unsigned int i)
Definition: gamma.h:66
T random()
Definition: gamma.h:113
N_length & a
Definition: tvecglobs.h:47
const T2 & b
Definition: minmax.h:48
_bz_global blitz::IndexPlaceholder< 10 > s
Definition: indexexpr.h:266
void setMean(T mean)
Definition: gamma.h:74
T snorm()
Definition: gamma.h:86
Gamma(T mean)
Definition: gamma.h:61
sharedState defaultState
Definition: default.h:55
_bz_global blitz::IndexPlaceholder< 7 > p
Definition: indexexpr.h:263
NormalUnit< T, IRNG, sharedState > normal_
Definition: gamma.h:106