docs: protected static double gammaQuantileNewtonRefinement(final double logpt, final double k, final double theta, final int maxit, double x) {
      final double EPS_N = 1e-15; // Precision threshold
      // 0 is not possible, try MIN_NORMAL instead
      if(x <= 0) {
        x = Double.MIN_NORMAL;
      }
      // Current estimation
      double logpc = logcdf(x, k, theta);
      if(x == Double.MIN_NORMAL && logpc > logpt * (1. + 1e-7)) {
        return 0.;
      }
      if(logpc == Double.NEGATIVE_INFINITY) {
        return 0.;
      }
      // Refine by newton iterations
      for(int i = 0; i < maxit; i++) {
        // Error of current approximation
        final double logpe = logpc - logpt;
        if(Math.abs(logpe) < Math.abs(EPS_N * logpt)) {
          break;
        }
        // Step size is controlled by PDF:
        final double g = logpdf(x, k, theta);
        if(g == Double.NEGATIVE_INFINITY) {
          break;
        }
        final double newx = x - logpe * FastMath.exp(logpc - g);
        // New estimate:
        logpc = logcdf(newx, k, theta);
        if(Math.abs(logpc - logpt) > Math.abs(logpe) || (i > 0 && Math.abs(logpc - logpt) == Math.abs(logpe))) {
          // no further improvement
          break;
        }
        x = newx;
      }
      return x;
    }