docs: @Reference(authors = "D. J. Best, D. E. Roberts", //
        title = "Algorithm AS 91: The percentage points of the χ² distribution", //
        booktitle = "Journal of the Royal Statistical Society. Series C (Applied Statistics)", //
        url = "https://doi.org/10.2307/2347113", //
        bibkey = "doi:10.2307/2347113")
    protected static double chisquaredProbitApproximation(final double p, double nu, double g) {
      final double EPS1 = 1e-14; // Approximation quality
      // Sanity checks
      if(Double.isNaN(p) || Double.isNaN(nu)) {
        return Double.NaN;
      }
      // Range check
      if(p <= 0) {
        return 0;
      }
      if(p >= 1) {
        return Double.POSITIVE_INFINITY;
      }
      // Invalid parameters
      if(nu <= 0) {
        return Double.NaN;
      }
      // Shape of gamma distribution, "XX" in AS 91
      final double k = 0.5 * nu;
  
      // For small chi squared values - AS 91
      final double logp = FastMath.log(p);
      if(nu < -1.24 * logp) {
        // FIXME: implement and use logGammap1 instead - more stable?
        //
        // final double lgam1pa = (alpha < 0.5) ? logGammap1(alpha) :
        // (FastMath.log(alpha) + g);
        // return FastMath.exp((lgam1pa + logp) / alpha + MathUtil.LOG2);
        // This is literal AS 91, above is the GNU R variant.
        return FastMath.pow(p * k * FastMath.exp(g + k * MathUtil.LOG2), 1. / k);
      }
      else if(nu > 0.32) {
        // Wilson and Hilferty estimate: - AS 91 at 3
        final double x = NormalDistribution.quantile(p, 0, 1);
        final double p1 = 2. / (9. * nu);
        final double a = x * FastMath.sqrt(p1) + 1 - p1;
        double ch = nu * a * a * a;
  
        // Better approximation for p tending to 1:
        if(ch > 2.2 * nu + 6) {
          ch = -2 * (FastMath.log1p(-p) - (k - 1) * FastMath.log(0.5 * ch) + g);
        }
        return ch;
      }
      else {
        // nu <= 0.32, AS 91 at 1
        final double C7 = 4.67, C8 = 6.66, C9 = 6.73, C10 = 13.32;
        final double ag = FastMath.log1p(-p) + g + (k - 1) * MathUtil.LOG2;
        double ch = 0.4;
        while(true) {
          final double p1 = 1 + ch * (C7 + ch);
          final double p2 = ch * (C9 + ch * (C8 + ch));
          final double t = -0.5 + (C7 + 2 * ch) / p1 - (C9 + ch * (C10 + 3 * ch)) / p2;
          final double delta = (1 - FastMath.exp(ag + 0.5 * ch) * p2 / p1) / t;
          ch -= delta;
          if(Math.abs(delta) <= EPS1 * Math.abs(ch)) {
            return ch;
          }
        }
      }
    }