#' Generic infimum‑type quantile function
#'
#' @param F   A function that evaluates the distribution function
#'            (CDF).  It must be monotone non‑decreasing in its first
#'            argument and return values in $0,1$.  Additional
#'            arguments for the CDF can be supplied via `...`.
#' @param p   Numeric vector of probabilities.  Values must lie in
#'            $0,1$; values outside this range give `NaN`.
#' @param lower  Left end of the search interval.  By default
#'            `-Inf` is used; if the CDF is known to be zero below a
#'            finite value you can give that value to speed up the
#'            search.
#' @param upper  Right end of the search interval.  By default
#'            `Inf` is used; if the CDF is known to be one above a
#'            finite value you can give that value.
#' @param tol   Desired absolute tolerance for the quantile.  The
#'            algorithm stops when the interval width is ≤ `tol`.
#' @param max.iter Maximum number of bisection iterations (safety
#'            guard).  The default is `1000`.
#' @param ...   Additional arguments passed to `F`.
#'
#' @return A numeric vector of the same length as `p` containing the
#'         infimum‑type quantiles.
#' @examples
#' ## 1. Continuous normal distribution (compare with qnorm)
#' qinf(pnorm, p = c(0.025, 0.5, 0.975), lower = -10, upper = 10)
#' qnorm(c(0.025, 0.5, 0.975))
#'
#' ## 2. Discrete distribution: Binomial(10, 0.3)
#' Fbinom <- function(x, size, prob) pbinom(floor(x), size, prob)
#' qinf(Fbinom, p = seq(0, 1, 0.1), size = 10, prob = 0.3)
#' qbinom(seq(0, 1, 0.1), size = 10, prob = 0.3)   # built‑in quantile
#'
#' ## 3. Mixed distribution (continuous + point mass at 0)
#' Fmix <- function(x, sigma = 1) {
#'   0.3 * (x >= 0) + 0.7 * pnorm(x, sd = sigma)
#' }
#' qinf(Fmix, p = c(0.1, 0.3, 0.5, 0.9), lower = -5, upper = 5)
#' @export
qinf <- function(F,
                 p,
                 lower = -Inf,
                 upper =  Inf,
                 tol = .Machine$double.eps^0.5,
                 max.iter = 1000L,
                 ...) {
  
  # 1.  Input checks -----------------------------------------------------------
  if (!is.function(F))
    stop("'F' must be a function that evaluates a CDF")
  if (!is.numeric(p))
    stop("'p' must be numeric")
  if (any(is.na(p)))               # keep NA positions in the output
    warning("NA probabilities supplied – corresponding quantiles will be NA")
  if (any(p < 0 | p > 1, na.rm = TRUE))
    stop("probabilities must lie in [0,1]")
  
  # 2.  Helper function for a single probability ---------------------------------
  
  one_quantile <- function(p_i) {
    ## Edge cases: 0 → leftmost point where F(x) >= 0,
    ##               1 → rightmost point where F(x) >= 1
    if (is.na(p_i)) return(NA_real_)
    if (p_i == 0) {
      ## The infimum of the set {x : F(x) >= 0} is the lower bound of the
      ## support.  If the user gave a finite lower bound we return it,
      ## otherwise we try to locate it by moving left until F changes.
      if (is.finite(lower)) return(lower)
      ## Search leftwards until we see a change in the CDF (or hit -Inf)
      x <- 0
      step <- 1
      while (is.finite(x) && F(x, ...) == 0) {
        x <- x - step
        step <- step * 2
        if (x < -1e12) break   # give up – treat as -Inf
      }
      return(x)
    }
    if (p_i == 1) {
      if (is.finite(upper)) return(upper)
      ## Search rightwards until F reaches 1
      x <- 0
      step <- 1
      while (is.finite(x) && F(x, ...) < 1) {
        x <- x + step
        step <- step * 2
        if (x > 1e12) break
      }
      return(x)
    }
    
    ## 2.1  Initialise the bracketing interval =================================
    lo <- lower
    hi <- upper
    
    ## If the interval is infinite we first try to find a finite bracket.
    ## This is done by exponential expansion from 0 (or from the sign of p)
    ## until the CDF straddles p.
    if (!is.finite(lo) || !is.finite(hi)) {
      # start from 0 (or any convenient point)
      centre <- 0
      # expand left if needed
      if (!is.finite(lo)) {
        step <- 1
        while (F(centre - step, ...) >= p_i) step <- step * 2
        lo <- centre - step
      }
      # expand right if needed
      if (!is.finite(hi)) {
        step <- 1
        while (F(centre + step, ...) < p_i) step <- step * 2
        hi <- centre + step
      }
    }
    
    ## 2.2  Bisection loop =====================================================
    ## we keep the leftmost point where F(q) >= p
    iter <- 0L
    while ((hi - lo) > tol && iter < max.iter) {
      mid <- (lo + hi) / 2
      fmid <- F(mid, ...)
      if (is.na(fmid)) {
        ## If the CDF returns NA (e.g. because mid is outside the domain)
        ## we shrink the interval towards the side that is known to be
        ## finite.
        if (is.finite(lo)) hi <- mid else lo <- mid
      } else if (fmid >= p_i) {
        ## We have reached (or passed) the target probability – move the
        ## upper bound leftwards to keep the *infimum*.
        hi <- mid
      } else {
        ## Still below the target – move the lower bound rightwards.
        lo <- mid
      }
      iter <- iter + 1L
    }
    
    ## 2.3  Return the upper bound (the smallest x with F(x) >= p) =============
    hi
  }
  
  # 3.  Vectorised call – preserve names / attributes of p ---------------------
  out <- vapply(p, one_quantile, numeric(1), USE.NAMES = FALSE)
  if (!is.null(names(p))) names(out) <- names(p)
  out
}