docs: protected static double rawQuantile(double p, double alpha, double beta, final double logbeta) {
      // Initial estimate for x
      double x;
      {
        // Very fast approximation of y.
        double tmp = FastMath.sqrt(-2 * FastMath.log(p));
        double y = tmp - (2.30753 + 0.27061 * tmp) / (1. + (0.99229 + 0.04481 * tmp) * tmp);
  
        if(alpha > 1 && beta > 1) {
          double r = (y * y - 3.) / 6.;
          double s = 1. / (alpha + alpha - 1.);
          double t = 1. / (beta + beta - 1.);
          double h = 2. / (s + t);
          double w = y * FastMath.sqrt(h + r) / h - (t - s) * (r + 5. / 6. - 2. / (3. * h));
          x = alpha / (alpha + beta * FastMath.exp(w + w));
        }
        else {
          double r = beta + beta;
          double t = 1. / (9. * beta);
          final double a = 1. - t + y * FastMath.sqrt(t);
          t = r * a * a * a;
          if(t <= 0.) {
            x = 1. - FastMath.exp((FastMath.log1p(-p) + FastMath.log(beta) + logbeta) / beta);
          }
          else {
            t = (4. * alpha + r - 2.) / t;
            if(t <= 1.) {
              x = FastMath.exp((FastMath.log(p * alpha) + logbeta) / alpha);
            }
            else {
              x = 1. - 2. / (t + 1.);
            }
          }
        }
        // Degenerate initial approximations
        if(x < 3e-308 || x > 1 - 2.22e-16) {
          x = 0.5;
        }
      }
  
      // Newon-Raphson method using the CDF
      {
        final double ialpha = 1 - alpha;
        final double ibeta = 1 - beta;
  
        // Desired accuracy, as from GNU R adoption of AS 109
        final double acu = Math.max(1e-300, FastMath.pow(10., -13 - 2.5 / (alpha * alpha) - .5 / (p * p)));
        double prevstep = 0., y = 0., stepsize = 1;
  
        for(int outer = 0; outer < 1000; outer++) {
          // Current CDF value
          double ynew = cdf(x, alpha, beta);
          if(Double.isInfinite(ynew)) { // Degenerated.
            return Double.NaN;
          }
          // Error gradient
          ynew = (ynew - p) * FastMath.exp(logbeta + ialpha * FastMath.log(x) + ibeta * FastMath.log1p(-x));
          if(ynew * y <= 0.) {
            prevstep = Math.max(Math.abs(stepsize), 3e-308);
          }
          // Inner loop: try different step sizes: y * 3^-i
          double g = 1, xnew = 0.;
          for(int inner = 0; inner < 1000; inner++) {
            stepsize = g * ynew;
            if(Math.abs(stepsize) < prevstep) {
              xnew = x - stepsize; // Candidate x
              if(xnew >= 0. && xnew <= 1.) {
                // Close enough
                if(prevstep <= acu || Math.abs(ynew) <= acu) {
                  return x;
                }
                if(xnew != 0. && xnew != 1.) {
                  break;
                }
              }
            }
            g /= 3.;
          }
          // Convergence
          if(Math.abs(xnew - x) < 1e-15 * x) {
            return x;
          }
          // Iterate with new values
          x = xnew;
          y = ynew;
        }
      }
      // Not converged in Newton-Raphson
      throw new ArithmeticException("Beta quantile computation did not converge.");
    }