docs: protected static double regularizedIncBetaQuadrature(double alpha, double beta, double x) {
      final double alphapbeta = alpha + beta;
      final double a1 = alpha - 1.0;
      final double b1 = beta - 1.0;
      final double mu = alpha / alphapbeta;
      final double lnmu = FastMath.log(mu);
      final double lnmuc = FastMath.log1p(-mu);
      double t = FastMath.sqrt(alpha * beta / (alphapbeta * alphapbeta * (alphapbeta + 1.0)));
      final double xu;
      if(x > alpha / alphapbeta) {
        if(x >= 1.0) {
          return 1.0;
        }
        xu = Math.min(1.0, Math.max(mu + 10.0 * t, x + 5.0 * t));
      }
      else {
        if(x <= 0.0) {
          return 0.0;
        }
        xu = Math.max(0.0, Math.min(mu - 10.0 * t, x - 5.0 * t));
      }
      double sum = 0.0;
      for(int i = 0; i < GAUSSLEGENDRE_Y.length; i++) {
        t = x + (xu - x) * GAUSSLEGENDRE_Y[i];
        sum += GAUSSLEGENDRE_W[i] * FastMath.exp(a1 * (FastMath.log(t) - lnmu) + b1 * (FastMath.log1p(-t) - lnmuc));
      }
      double ans = sum * (xu - x) * FastMath.exp(a1 * lnmu - GammaDistribution.logGamma(alpha) + b1 * lnmuc - GammaDistribution.logGamma(beta) + GammaDistribution.logGamma(alphapbeta));
      return ans > 0 ? 1.0 - ans : -ans;
    }