36 #ifndef BZ_RANDOM_BETA
37 #define BZ_RANDOM_BETA
39 #ifndef BZ_RANDOM_UNIFORM
43 #ifndef BZ_NUMINQUIRE_H
44 #include <blitz/numinquire.h>
49 template<typename T =
double, typename IRNG =
defaultIRNG,
60 infnty = 0.3 * huge(T());
61 minlog = 0.085 * tiny(T());
67 void setParameters(T
a, T b)
83 template<
typename T,
typename IRNG,
typename stateTag>
94 static T olda = minlog;
95 static T oldb = minlog;
96 static T genbet,
a,alpha,b,beta,delta,gamma,k1,k2,
r,
s,
t,u1,u2,v,w,y,z;
102 qsame = (olda == aa) && (oldb == bb);
106 BZPRECHECK((aa > minlog) && (bb > minlog),
107 "Beta RNG: parameters must be >= " << minlog << endl
108 <<
"Parameters are aa=" << aa <<
" and bb=" << bb << endl);
113 if (!(min(aa,bb) > 1.0))
123 beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
132 v = beta*log(u1/(1.0-u1));
134 if(v > expmax)
goto S55;
140 if(w > infnty/a)
goto S55;
147 r = gamma*v-1.3862944;
152 if(s+2.609438 >= 5.0*z)
goto S70;
167 if(alpha/(b+w) < minlog)
goto S40;
168 if(r+alpha*log(alpha/(b+w)) < t)
goto S40;
173 if(!(aa == a))
goto S80;
191 k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
192 k2 = 0.25+(0.5+0.25/delta)*b;
200 if(u1 >= 0.5)
goto S130;
206 if(0.25*u2+z-y >= k1)
goto S120;
213 if(!(z <= 0.25))
goto S160;
214 v = beta*log(u1/(1.0-u1));
219 if(a > 1.0)
goto S135;
221 if(v > expmax)
goto S132;
226 if(w > expmax)
goto S140;
231 if(v > expmax)
goto S140;
233 if(w > infnty/a)
goto S140;
250 if(z >= k2)
goto S120;
256 v = beta*log(u1/(1.0-u1));
258 if(a > 1.0)
goto S175;
260 if(v > expmax)
goto S172;
265 if(w > expmax)
goto S180;
270 if(v > expmax)
goto S180;
272 if(w > infnty/a)
goto S180;
290 if(alpha/(b+w) < minlog)
goto S120;
291 if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z))
goto S120;
296 if(!(a == aa))
goto S210;
308 #endif // BZ_RANDOM_BETA