42#ifndef BZ_RANDOM_UNIFORM
46#ifndef BZ_NUMINQUIRE_H
47 #include <blitz/numinquire.h>
52template<
typename T = double,
typename IRNG =
defaultIRNG,
75 infnty = 0.3 * blitz::huge(T());
76 minlog = 0.085 * blitz::tiny(T());
90template<
typename T,
typename IRNG,
typename stateTag>
101 static T olda = minlog;
102 static T oldb = minlog;
103 static T genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
109 qsame = (olda == aa) && (oldb == bb);
113 BZPRECHECK((aa > minlog) && (bb > minlog),
114 "Beta RNG: parameters must be >= " << minlog << endl
115 <<
"Parameters are aa=" << aa <<
" and bb=" << bb << endl);
120 if (!(min(aa,bb) > 1.0))
130 beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
139 v = beta*log(u1/(1.0-u1));
141 if(v > expmax)
goto S55;
147 if(w > infnty/a)
goto S55;
154 r = gamma*v-1.3862944;
159 if(s+2.609438 >= 5.0*z)
goto S70;
174 if(alpha/(b+w) < minlog)
goto S40;
175 if(r+alpha*log(alpha/(b+w)) < t)
goto S40;
180 if(!(aa == a))
goto S80;
198 k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
199 k2 = 0.25+(0.5+0.25/delta)*b;
207 if(u1 >= 0.5)
goto S130;
213 if(0.25*u2+z-y >= k1)
goto S120;
220 if(!(z <= 0.25))
goto S160;
221 v = beta*log(u1/(1.0-u1));
226 if(a > 1.0)
goto S135;
228 if(v > expmax)
goto S132;
233 if(w > expmax)
goto S140;
238 if(v > expmax)
goto S140;
240 if(w > infnty/a)
goto S140;
257 if(z >= k2)
goto S120;
263 v = beta*log(u1/(1.0-u1));
265 if(a > 1.0)
goto S175;
267 if(v > expmax)
goto S172;
272 if(w > expmax)
goto S180;
277 if(v > expmax)
goto S180;
279 if(w > infnty/a)
goto S180;
297 if(alpha/(b+w) < minlog)
goto S120;
298 if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z))
goto S120;
303 if(!(a == aa))
goto S210;
T random()
Definition: beta.h:91
T minlog
Definition: beta.h:87
void setParameters(T a, T b)
Definition: beta.h:71
Beta(T a, T b)
Definition: beta.h:59
T aa
Definition: beta.h:86
Beta(T a, T b, unsigned int i)
Definition: beta.h:64
T T_numtype
Definition: beta.h:57
T infnty
Definition: beta.h:87
T ranf()
Definition: beta.h:81
T expmax
Definition: beta.h:87
T bb
Definition: beta.h:86
MersenneTwister defaultIRNG
Definition: default.h:120