00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "CLHEP/Random/RandBinomial.h"
00020 #include "CLHEP/Random/DoubConv.h"
00021 #include <algorithm>
00022 #include <cmath>
00023
00024 namespace CLHEP {
00025
00026 std::string RandBinomial::name() const {return "RandBinomial";}
00027 HepRandomEngine & RandBinomial::engine() {return *localEngine;}
00028
00029 RandBinomial::~RandBinomial() {
00030 }
00031
00032 double RandBinomial::shoot( HepRandomEngine *anEngine, long n,
00033 double p ) {
00034 return genBinomial( anEngine, n, p );
00035 }
00036
00037 double RandBinomial::shoot( long n, double p ) {
00038 HepRandomEngine *anEngine = HepRandom::getTheEngine();
00039 return genBinomial( anEngine, n, p );
00040 }
00041
00042 double RandBinomial::fire( long n, double p ) {
00043 return genBinomial( localEngine.get(), n, p );
00044 }
00045
00046 void RandBinomial::shootArray( const int size, double* vect,
00047 long n, double p )
00048 {
00049 for( double* v = vect; v != vect+size; ++v )
00050 *v = shoot(n,p);
00051 }
00052
00053 void RandBinomial::shootArray( HepRandomEngine* anEngine,
00054 const int size, double* vect,
00055 long n, double p )
00056 {
00057 for( double* v = vect; v != vect+size; ++v )
00058 *v = shoot(anEngine,n,p);
00059 }
00060
00061 void RandBinomial::fireArray( const int size, double* vect)
00062 {
00063 for( double* v = vect; v != vect+size; ++v )
00064 *v = fire(defaultN,defaultP);
00065 }
00066
00067 void RandBinomial::fireArray( const int size, double* vect,
00068 long n, double p )
00069 {
00070 for( double* v = vect; v != vect+size; ++v )
00071 *v = fire(n,p);
00072 }
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092 static double StirlingCorrection(long int k)
00093 {
00094 #define C1 8.33333333333333333e-02 // +1/12
00095 #define C3 -2.77777777777777778e-03 // -1/360
00096 #define C5 7.93650793650793651e-04 // +1/1260
00097 #define C7 -5.95238095238095238e-04 // -1/1680
00098
00099 static double c[31] = { 0.0,
00100 8.106146679532726e-02, 4.134069595540929e-02,
00101 2.767792568499834e-02, 2.079067210376509e-02,
00102 1.664469118982119e-02, 1.387612882307075e-02,
00103 1.189670994589177e-02, 1.041126526197209e-02,
00104 9.255462182712733e-03, 8.330563433362871e-03,
00105 7.573675487951841e-03, 6.942840107209530e-03,
00106 6.408994188004207e-03, 5.951370112758848e-03,
00107 5.554733551962801e-03, 5.207655919609640e-03,
00108 4.901395948434738e-03, 4.629153749334029e-03,
00109 4.385560249232324e-03, 4.166319691996922e-03,
00110 3.967954218640860e-03, 3.787618068444430e-03,
00111 3.622960224683090e-03, 3.472021382978770e-03,
00112 3.333155636728090e-03, 3.204970228055040e-03,
00113 3.086278682608780e-03, 2.976063983550410e-03,
00114 2.873449362352470e-03, 2.777674929752690e-03,
00115 };
00116 double r, rr;
00117
00118 if (k > 30L) {
00119 r = 1.0 / (double) k; rr = r * r;
00120 return(r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
00121 }
00122 else return(c[k]);
00123 }
00124
00125 double RandBinomial::genBinomial( HepRandomEngine *anEngine, long n, double p )
00126 {
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 #define C1_3 0.33333333333333333
00173 #define C5_8 0.62500000000000000
00174 #define C1_6 0.16666666666666667
00175 #define DMAX_KM 20L
00176
00177 static long int n_last = -1L, n_prev = -1L;
00178 static double par,np,p0,q,p_last = -1.0, p_prev = -1.0;
00179 static long b,m,nm;
00180 static double pq, rc, ss, xm, xl, xr, ll, lr, c,
00181 p1, p2, p3, p4, ch;
00182
00183 long bh,i, K, Km, nK;
00184 double f, rm, U, V, X, T, E;
00185
00186 if (n != n_last || p != p_last)
00187 {
00188 n_last = n;
00189 p_last = p;
00190 par=std::min(p,1.0-p);
00191 q=1.0-par;
00192 np = n*par;
00193
00194
00195
00196 if( np <= 0.0 ) return (-1.0);
00197
00198 rm = np + par;
00199 m = (long int) rm;
00200 if (np<10)
00201 {
00202 p0=std::exp(n*std::log(q));
00203 bh=(long int)(np+10.0*std::sqrt(np*q));
00204 b=std::min(n,bh);
00205 }
00206 else
00207 {
00208 rc = (n + 1.0) * (pq = par / q);
00209 ss = np * q;
00210 i = (long int) (2.195*std::sqrt(ss) - 4.6*q);
00211 xm = m + 0.5;
00212 xl = (double) (m - i);
00213 xr = (double) (m + i + 1L);
00214 f = (rm - xl) / (rm - xl*par); ll = f * (1.0 + 0.5*f);
00215 f = (xr - rm) / (xr * q); lr = f * (1.0 + 0.5*f);
00216 c = 0.134 + 20.5/(15.3 + (double) m);
00217
00218 p1 = i + 0.5;
00219 p2 = p1 * (1.0 + c + c);
00220 p3 = p2 + c/ll;
00221 p4 = p3 + c/lr;
00222 }
00223 }
00224 if (np<10)
00225 {
00226 double pk;
00227
00228 K=0;
00229 pk=p0;
00230 U=anEngine->flat();
00231 while (U>pk)
00232 {
00233 ++K;
00234 if (K>b)
00235 {
00236 U=anEngine->flat();
00237 K=0;
00238 pk=p0;
00239 }
00240 else
00241 {
00242 U-=pk;
00243 pk=(double)(((n-K+1)*par*pk)/(K*q));
00244 }
00245 }
00246 return ((p>0.5) ? (double)(n-K):(double)K);
00247 }
00248
00249 for (;;)
00250 {
00251 V = anEngine->flat();
00252 if ((U = anEngine->flat() * p4) <= p1)
00253 {
00254 K=(long int) (xm - U + p1*V);
00255 return ((p>0.5) ? (double)(n-K):(double)K);
00256 }
00257 if (U <= p2)
00258 {
00259 X = xl + (U - p1)/c;
00260 if ((V = V*c + 1.0 - std::fabs(xm - X)/p1) >= 1.0) continue;
00261 K = (long int) X;
00262 }
00263 else if (U <= p3)
00264 {
00265 if ((X = xl + std::log(V)/ll) < 0.0) continue;
00266 K = (long int) X;
00267 V *= (U - p2) * ll;
00268 }
00269 else
00270 {
00271 if ((K = (long int) (xr - std::log(V)/lr)) > n) continue;
00272 V *= (U - p3) * lr;
00273 }
00274
00275
00276 if ((Km = std::labs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss)
00277 {
00278
00279
00280 f = 1.0;
00281 if (m < K)
00282 {
00283 for (i = m; i < K; )
00284 {
00285 if ((f *= (rc / ++i - pq)) < V) break;
00286 }
00287 }
00288 else
00289 {
00290 for (i = K; i < m; )
00291 {
00292 if ((V *= (rc / ++i - pq)) > f) break;
00293 }
00294 }
00295 if (V <= f) break;
00296 }
00297 else
00298 {
00299
00300
00301 V = std::log(V);
00302 T = - Km * Km / (ss + ss);
00303 E = (Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5);
00304 if (V <= T - E) break;
00305 if (V <= T + E)
00306 {
00307 if (n != n_prev || par != p_prev)
00308 {
00309 n_prev = n;
00310 p_prev = par;
00311
00312 nm = n - m + 1L;
00313 ch = xm * std::log((m + 1.0)/(pq * nm)) +
00314 StirlingCorrection(m + 1L) + StirlingCorrection(nm);
00315 }
00316 nK = n - K + 1L;
00317
00318
00319
00320 if (V <= ch + (n + 1.0)*std::log((double) nm / (double) nK) +
00321 (K + 0.5)*std::log(nK * pq / (K + 1.0)) -
00322 StirlingCorrection(K + 1L) - StirlingCorrection(nK)) break;
00323 }
00324 }
00325 }
00326 return ((p>0.5) ? (double)(n-K):(double)K);
00327 }
00328
00329 std::ostream & RandBinomial::put ( std::ostream & os ) const {
00330 int pr=os.precision(20);
00331 std::vector<unsigned long> t(2);
00332 os << " " << name() << "\n";
00333 os << "Uvec" << "\n";
00334 t = DoubConv::dto2longs(defaultP);
00335 os << defaultN << " " << defaultP << " " << t[0] << " " << t[1] << "\n";
00336 os.precision(pr);
00337 return os;
00338 }
00339
00340 std::istream & RandBinomial::get ( std::istream & is ) {
00341 std::string inName;
00342 is >> inName;
00343 if (inName != name()) {
00344 is.clear(std::ios::badbit | is.rdstate());
00345 std::cerr << "Mismatch when expecting to read state of a "
00346 << name() << " distribution\n"
00347 << "Name found was " << inName
00348 << "\nistream is left in the badbit state\n";
00349 return is;
00350 }
00351 if (possibleKeywordInput(is, "Uvec", defaultN)) {
00352 std::vector<unsigned long> t(2);
00353 is >> defaultN >> defaultP;
00354 is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t);
00355 return is;
00356 }
00357
00358 is >> defaultP;
00359 return is;
00360 }
00361
00362
00363 }