19 #include "CLHEP/Random/RandBinomial.h"
20 #include "CLHEP/Random/DoubConv.h"
21 #include "CLHEP/Utility/thread_local.h"
30 RandBinomial::~RandBinomial() {
35 return genBinomial( anEngine, n, p );
39 HepRandomEngine *anEngine = HepRandom::getTheEngine();
40 return genBinomial( anEngine, n, p );
43 double RandBinomial::fire(
long n,
double p ) {
44 return genBinomial( localEngine.get(),
n, p );
47 void RandBinomial::shootArray(
const int size,
double* vect,
50 for(
double* v = vect; v != vect+size; ++v )
54 void RandBinomial::shootArray( HepRandomEngine* anEngine,
55 const int size,
double* vect,
58 for(
double* v = vect; v != vect+size; ++v )
59 *v =
shoot(anEngine,n,p);
62 void RandBinomial::fireArray(
const int size,
double* vect)
64 for(
double* v = vect; v != vect+size; ++v )
65 *v = fire(defaultN,defaultP);
68 void RandBinomial::fireArray(
const int size,
double* vect,
71 for(
double* v = vect; v != vect+size; ++v )
95 #define C1 8.33333333333333333e-02 // +1/12
96 #define C3 -2.77777777777777778e-03 // -1/360
97 #define C5 7.93650793650793651e-04 // +1/1260
98 #define C7 -5.95238095238095238e-04 // -1/1680
100 static const double c[31] = { 0.0,
101 8.106146679532726e-02, 4.134069595540929e-02,
102 2.767792568499834e-02, 2.079067210376509e-02,
103 1.664469118982119e-02, 1.387612882307075e-02,
104 1.189670994589177e-02, 1.041126526197209e-02,
105 9.255462182712733e-03, 8.330563433362871e-03,
106 7.573675487951841e-03, 6.942840107209530e-03,
107 6.408994188004207e-03, 5.951370112758848e-03,
108 5.554733551962801e-03, 5.207655919609640e-03,
109 4.901395948434738e-03, 4.629153749334029e-03,
110 4.385560249232324e-03, 4.166319691996922e-03,
111 3.967954218640860e-03, 3.787618068444430e-03,
112 3.622960224683090e-03, 3.472021382978770e-03,
113 3.333155636728090e-03, 3.204970228055040e-03,
114 3.086278682608780e-03, 2.976063983550410e-03,
115 2.873449362352470e-03, 2.777674929752690e-03,
120 r = 1.0 / (double) k; rr = r * r;
121 return(r*(
C1 + rr*(
C3 + rr*(
C5 + rr*
C7))));
126 double RandBinomial::genBinomial( HepRandomEngine *anEngine,
long n,
double p )
173 #define C1_3 0.33333333333333333
174 #define C5_8 0.62500000000000000
175 #define C1_6 0.16666666666666667
178 static CLHEP_THREAD_LOCAL
long int n_last = -1
L, n_prev = -1
L;
179 static CLHEP_THREAD_LOCAL
double par,np,p0,q,p_last = -1.0, p_prev = -1.0;
180 static CLHEP_THREAD_LOCAL
long b,
m,
nm;
181 static CLHEP_THREAD_LOCAL
double pq, rc, ss, xm, xl, xr, ll, lr, c,
184 long bh,i, K, Km, nK;
185 double f, rm, U, V, X, T, E;
187 if (n != n_last || p != p_last)
197 if( np <= 0.0 )
return (-1.0);
203 p0=std::exp(n*std::log(q));
204 bh=(
long int)(np+10.0*std::sqrt(np*q));
209 rc = (n + 1.0) * (pq = par / q);
211 i = (
long int) (2.195*std::sqrt(ss) - 4.6*q);
213 xl = (double) (m - i);
214 xr = (double) (m + i + 1
L);
215 f = (rm - xl) / (rm - xl*par); ll = f * (1.0 + 0.5*f);
216 f = (xr - rm) / (xr * q); lr = f * (1.0 + 0.5*f);
217 c = 0.134 + 20.5/(15.3 + (double) m);
220 p2 = p1 * (1.0 + c + c);
225 if( np <= 0.0 )
return (-1.0);
245 pk=(double)(((n-K+1)*par*pk)/(K*q));
248 return ((p>0.5) ? (double)(n-K):(double)K);
253 V = anEngine->flat();
254 if ((U = anEngine->flat() * p4) <= p1)
256 K=(
long int) (xm - U + p1*V);
257 return ((p>0.5) ? (double)(n-K):(double)K);
262 if ((V = V*c + 1.0 - std::fabs(xm - X)/p1) >= 1.0)
continue;
267 if ((X = xl + std::log(V)/ll) < 0.0)
continue;
273 if ((K = (
long int) (xr - std::log(V)/lr)) >
n)
continue;
278 if ((Km = std::labs(K - m)) <=
DMAX_KM || Km + Km + 2
L >= ss)
287 if ((f *= (rc / ++i - pq)) < V)
break;
294 if ((V *= (rc / ++i - pq)) > f)
break;
304 T = - Km * Km / (ss + ss);
305 E = (Km / ss) * ((Km * (Km *
C1_3 +
C5_8) +
C1_6) / ss + 0.5);
306 if (V <= T - E)
break;
309 if (n != n_prev || par != p_prev)
315 ch = xm * std::log((m + 1.0)/(pq * nm)) +
322 if (V <= ch + (n + 1.0)*std::log((
double) nm / (
double) nK) +
323 (K + 0.5)*std::log(nK * pq / (K + 1.0)) -
328 return ((p>0.5) ? (
double)(n-K):(
double)K);
331 std::ostream & RandBinomial::put ( std::ostream & os )
const {
332 int pr=os.precision(20);
333 std::vector<unsigned long> t(2);
334 os <<
" " <<
name() <<
"\n";
335 os <<
"Uvec" <<
"\n";
336 t = DoubConv::dto2longs(defaultP);
337 os << defaultN <<
" " << defaultP <<
" " << t[0] <<
" " << t[1] <<
"\n";
342 std::istream & RandBinomial::get ( std::istream & is ) {
345 if (inName !=
name()) {
346 is.clear(std::ios::badbit | is.rdstate());
347 std::cerr <<
"Mismatch when expecting to read state of a "
348 <<
name() <<
" distribution\n"
349 <<
"Name found was " << inName
350 <<
"\nistream is left in the badbit state\n";
353 if (possibleKeywordInput(is,
"Uvec", defaultN)) {
354 std::vector<unsigned long> t(2);
355 is >> defaultN >> defaultP;
356 is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t);
ThreeVector shoot(const G4int Ap, const G4int Af)
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
static double StirlingCorrection(long int k)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments