| 610 | } |
| 611 | |
| 612 | RAND_INT_TYPE random_binomial_btpe(bitgen_t *bitgen_state, RAND_INT_TYPE n, |
| 613 | double p, binomial_t *binomial) { |
| 614 | double r, q, fm, p1, xm, xl, xr, c, laml, lamr, p2, p3, p4; |
| 615 | double a, u, v, s, F, rho, t, A, nrq, x1, x2, f1, f2, z, z2, w, w2, x; |
| 616 | RAND_INT_TYPE m, y, k, i; |
| 617 | |
| 618 | if (!(binomial->has_binomial) || (binomial->nsave != n) || |
| 619 | (binomial->psave != p)) { |
| 620 | /* initialize */ |
| 621 | binomial->nsave = n; |
| 622 | binomial->psave = p; |
| 623 | binomial->has_binomial = 1; |
| 624 | binomial->r = r = MIN(p, 1.0 - p); |
| 625 | binomial->q = q = 1.0 - r; |
| 626 | binomial->fm = fm = n * r + r; |
| 627 | binomial->m = m = (RAND_INT_TYPE)floor(binomial->fm); |
| 628 | binomial->p1 = p1 = floor(2.195 * sqrt(n * r * q) - 4.6 * q) + 0.5; |
| 629 | binomial->xm = xm = m + 0.5; |
| 630 | binomial->xl = xl = xm - p1; |
| 631 | binomial->xr = xr = xm + p1; |
| 632 | binomial->c = c = 0.134 + 20.5 / (15.3 + m); |
| 633 | a = (fm - xl) / (fm - xl * r); |
| 634 | binomial->laml = laml = a * (1.0 + a / 2.0); |
| 635 | a = (xr - fm) / (xr * q); |
| 636 | binomial->lamr = lamr = a * (1.0 + a / 2.0); |
| 637 | binomial->p2 = p2 = p1 * (1.0 + 2.0 * c); |
| 638 | binomial->p3 = p3 = p2 + c / laml; |
| 639 | binomial->p4 = p4 = p3 + c / lamr; |
| 640 | } else { |
| 641 | r = binomial->r; |
| 642 | q = binomial->q; |
| 643 | fm = binomial->fm; |
| 644 | m = binomial->m; |
| 645 | p1 = binomial->p1; |
| 646 | xm = binomial->xm; |
| 647 | xl = binomial->xl; |
| 648 | xr = binomial->xr; |
| 649 | c = binomial->c; |
| 650 | laml = binomial->laml; |
| 651 | lamr = binomial->lamr; |
| 652 | p2 = binomial->p2; |
| 653 | p3 = binomial->p3; |
| 654 | p4 = binomial->p4; |
| 655 | } |
| 656 | |
| 657 | /* sigh ... */ |
| 658 | Step10: |
| 659 | nrq = n * r * q; |
| 660 | u = next_double(bitgen_state) * p4; |
| 661 | v = next_double(bitgen_state); |
| 662 | if (u > p1) |
| 663 | goto Step20; |
| 664 | y = (RAND_INT_TYPE)floor(xm - p1 * v + u); |
| 665 | goto Step60; |
| 666 | |
| 667 | Step20: |
| 668 | if (u > p2) |
| 669 | goto Step30; |
no test coverage detected