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")
    public static double quantile(double p, double k, double theta) {
      final double EPS2 = 5e-7; // final precision of AS 91
      final int MAXIT = 1000;
  
      // Avoid degenerates
      if(!(p >= 0) || p > 1 || Double.isNaN(k) || Double.isNaN(theta)) {
        return Double.NaN;
      }
      // Range check
      if(p == 0) {
        return 0;
      }
      if(p == 1) {
        return Double.POSITIVE_INFINITY;
      }
      // Shape parameter check
      if(k < 0 || theta <= 0) {
        return Double.NaN;
      }
      // Corner case - all at 0
      if(k == 0) {
        return 0.;
      }
  
      int max_newton_iterations = 1;
      // For small values, ensure some refinement iterations
      if(k < 1e-10) {
        max_newton_iterations = 7;
      }
  
      final double g = logGamma(k); // == logGamma(v/2)
  
      // Phase I, an initial rough approximation
      // First half of AS 91
      double ch = chisquaredProbitApproximation(p, 2 * k, g);
      // Second half of AS 91 follows:
      // Refine ChiSquared approximation
      chisq: {
        if(Double.isInfinite(ch)) {
          // Cannot refine infinity
          max_newton_iterations = 0;
          break chisq;
        }
        if(ch < EPS2) {
          // Do not iterate, but refine with Newton method
          max_newton_iterations = 20;
          break chisq;
        }
        if(p > 1 - 1e-14 || p < 1e-100) {
          // Not in appropriate value range for AS 91
          max_newton_iterations = 20;
          break chisq;
        }
  
        // Phase II: Iteration
        final double c = k - 1;
        final double ch0 = ch; // backup initial approximation
        for(int i = 1; i <= MAXIT; i++) {
          final double q = ch; // previous approximation
          final double p1 = 0.5 * ch;
          final double p2 = p - regularizedGammaP(k, p1);
          if(Double.isInfinite(p2) || ch <= 0) {
            ch = ch0;
            max_newton_iterations = 27;
            break chisq;
          }
          { // Taylor series of AS 91: iteration via "goto 4"
            final double t = p2 * FastMath.exp(k * MathUtil.LOG2 + g + p1 - c * FastMath.log(ch));
            final double b = t / ch;
            final double a = 0.5 * t - b * c;
            final double s1 = (210. + a * (140. + a * (105. + a * (84. + a * (70. + 60. * a))))) / 420.;
            final double s2 = (420. + a * (735. + a * (966. + a * (1141. + 1278 * a)))) / 2520.;
            final double s3 = (210. + a * (462. + a * (707. + 932. * a))) / 2520.;
            final double s4 = (252. + a * (672. + 1182. * a) + c * (294. + a * (889. + 1740. * a))) / 5040.;
            final double s5 = (84. + 2264. * a + c * (1175. + 606. * a)) / 2520.;
            final double s6 = (120. + c * (346. + 127. * c)) / 5040.;
            ch += t * (1 + 0.5 * t * s1 - b * c * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
          }
          if(Math.abs(q - ch) < EPS2 * ch) {
            break chisq;
          }
          // Divergence treatment, from GNU R
          if(Math.abs(q - ch) > 0.1 * Math.abs(ch)) {
            ch = ((ch < q) ? 0.9 : 1.1) * q;
          }
        }
        LoggingUtil.warning("No convergence in AS 91 Gamma probit.");
        // no convergence in MAXIT iterations -- but we add Newton now...
      }
      double x = 0.5 * ch / theta;
      if(max_newton_iterations > 0) {
        // Refine result using final Newton steps.
        // TODO: add unit tests that show an improvement! Maybe in logscale only?
        x = gammaQuantileNewtonRefinement(FastMath.log(p), k, theta, max_newton_iterations, x);
      }
      return x;
    }