double BetaIncompleteInvImpl()

in absl/random/internal/distribution_test_util.cc [256:363]

• 94 lines of code
• 20 McCabe index (conditional complexity)

double BetaIncompleteInvImpl(const double p, const double q, const double beta,
                             const double alpha) {
  if (alpha < 0.5) {
    // Inverse Incomplete beta function is symmetrical, return the complement.
    return 1. - BetaIncompleteInvImpl(q, p, beta, 1. - alpha);
  }
  const double kErr = 1e-14;
  double value = kErr;

  // Compute the initial estimate.
  {
    double r = std::sqrt(-std::log(alpha * alpha));
    double y =
        r - fma(r, 0.27061, 2.30753) / fma(r, fma(r, 0.04481, 0.99229), 1.0);
    if (p > 1. && q > 1.) {
      r = (y * y - 3.) / 6.;
      double s = 1. / (p + p - 1.);
      double t = 1. / (q + q - 1.);
      double h = 2. / s + t;
      double w =
          y * std::sqrt(h + r) / h - (t - s) * (r + 5. / 6. - t / (3. * h));
      value = p / (p + q * std::exp(w + w));
    } else {
      r = q + q;
      double t = 1.0 / (9. * q);
      double u = 1.0 - t + y * std::sqrt(t);
      t = r * (u * u * u);
      if (t <= 0) {
        value = 1.0 - std::exp((std::log((1.0 - alpha) * q) + beta) / q);
      } else {
        t = (4.0 * p + r - 2.0) / t;
        if (t <= 1) {
          value = std::exp((std::log(alpha * p) + beta) / p);
        } else {
          value = 1.0 - 2.0 / (t + 1.0);
        }
      }
    }
  }

  // Solve for x using a modified newton-raphson method using the function
  // BetaIncomplete.
  {
    value = std::max(value, kErr);
    value = std::min(value, 1.0 - kErr);

    const double r = 1.0 - p;
    const double t = 1.0 - q;
    double y;
    double yprev = 0;
    double sq = 1;
    double prev = 1;
    for (;;) {
      if (value < 0 || value > 1.0) {
        // Error case; value went infinite.
        return std::numeric_limits<double>::infinity();
      } else if (value == 0 || value == 1) {
        y = value;
      } else {
        y = BetaIncompleteImpl(value, p, q, beta);
        if (!std::isfinite(y)) {
          return y;
        }
      }
      y = (y - alpha) *
          std::exp(beta + r * std::log(value) + t * std::log(1.0 - value));
      if (y * yprev <= 0) {
        prev = std::max(sq, std::numeric_limits<double>::min());
      }
      double g = 1.0;
      for (;;) {
        const double adj = g * y;
        const double adj_sq = adj * adj;
        if (adj_sq >= prev) {
          g = g / 3.0;
          continue;
        }
        const double tx = value - adj;
        if (tx < 0 || tx > 1) {
          g = g / 3.0;
          continue;
        }
        if (prev < kErr) {
          return value;
        }
        if (y * y < kErr) {
          return value;
        }
        if (tx == value) {
          return value;
        }
        if (tx == 0 || tx == 1) {
          g = g / 3.0;
          continue;
        }
        value = tx;
        yprev = y;
        break;
      }
    }
  }

  // NOTES: See also: Asymptotic inversion of the incomplete beta function.
  // https://core.ac.uk/download/pdf/82140723.pdf
  //
  // NOTE: See the Boost library documentation as well:
  // https://www.boost.org/doc/libs/1_52_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_beta/ibeta_function.html
}