#  File src/library/stats/R/binom.test.R
#  Part of the R package, https://www.R-project.org
#
#  Copyright (C) 1995-2026 The R Core Team
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

binom.test <-
function(x, n, p = 0.5, alternative = c("two.sided", "less", "greater"),
         conf.level = 0.95, two.sided.method = c("minlike", "central"))
{
    DNAME <- deparse1(substitute(x))
    xr <- round(x)

    if(any(is.na(x) | (x < 0)) || max(abs(x-xr)) > 1e-7)
        stop("'x' must be nonnegative and integer")
    x <- xr
    if(length(x) == 2L) {
        ## x gives successes and failures
        n <- sum(x)
        x <- x[1L]
    }
    else if(length(x) == 1L) {
        ## x gives successes, n gives trials
        nr <- round(n)
        if((length(n) > 1L) || is.na(n) || (n < 1) || abs(n-nr) > 1e-7
           || (x > nr))
            stop("'n' must be a positive integer >= 'x'")
        DNAME <- paste(DNAME, "and", deparse1(substitute(n)))
        n <- nr
    }
    else
        stop("incorrect length of 'x'")

    if(!missing(p) && (length(p) > 1L || is.na(p) || p < 0 || p > 1))
        stop ("'p' must be a single number between 0 and 1")
    alternative <- match.arg(alternative)

    if(!missing(conf.level) &&
       (length(conf.level) != 1 || !is.finite(conf.level) ||
        conf.level < 0 || conf.level > 1))
        stop("'conf.level' must be a single number between 0 and 1")

    if(alternative == "two.sided")
        two.sided.method <- match.arg(two.sided.method)

    PVAL <- switch(alternative,
                   less = pbinom(x, n, p),
                   greater = pbinom(x - 1, n, p, lower.tail = FALSE),
                   two.sided =
                       .binom_test_two_sided_pval(x, n, p,
                                                  two.sided.method)
                   )
    ## Determine p s.t. Prob(B(n,p) >= x) = alpha.
    ## Use that for x > 0,
    ##   Prob(B(n,p) >= x) = pbeta(p, x, n - x + 1).
    p.L <- function(x, alpha) {
        if(x == 0)                      # No solution
            0
        else
            qbeta(alpha, x, n - x + 1)
    }
    ## Determine p s.t. Prob(B(n,p) <= x) = alpha.
    ## Use that for x < n,
    ##   Prob(B(n,p) <= x) = 1 - pbeta(p, x + 1, n - x).
    p.U <- function(x, alpha) {
        if(x == n)                      # No solution
            1
        else
            qbeta(1 - alpha, x + 1, n - x)
    }
    CINT <- switch(alternative,
                   less = c(0, p.U(x, 1 - conf.level)),
                   greater = c(p.L(x, 1 - conf.level), 1),
                   two.sided = {
                       alpha <- (1 - conf.level) / 2
                       c(p.L(x, alpha), p.U(x, alpha))
                   })
    attr(CINT, "conf.level") <- conf.level

    ESTIMATE <- x / n

    names(x) <- "number of successes"	# or simply "x" ??
    names(n) <- "number of trials"	# or simply "n" ??
    names(ESTIMATE) <-
    names(p) <- "probability of success"# or simply "p" ??

    structure(list(statistic = x,
                   parameter = n,
                   p.value = PVAL,
                   conf.int = CINT,
                   estimate = ESTIMATE,
                   null.value = p,
                   alternative = alternative,
                   method = "Exact binomial test",
                   data.name = DNAME),
              class = "htest")
}

.binom_test_two_sided_pval <-
function(x, n, p, method)
{
    switch(method,
           "central" = .binom_test_two_sided_pval_central(x, n, p),
           "minlike" = .binom_test_two_sided_pval_minlike(x, n, p))
}

.binom_test_two_sided_pval_central <-
function(x, n, p)
    min(2 * min(pbinom(x, n, p),
                pbinom(x - 1, n, p, lower.tail = FALSE)),
        1)

.binom_test_two_sided_pval_minlike <-
function(x, n, p, bisect = (n >= 1234))
{
    if(p == 0)
        return(x == 0)
    if(p == 1)
        return(x == n)
    log_dx <- dbinom(x, n, p, log = TRUE)
    if(!is.finite(log_dx))
        return(0)
    relErr <- 1 + 1e-07
    log_thresh <- log_dx + log(relErr)
    if(bisect) {
        ## Fast O(log n) path for large n.
        ## Contributed by Ravi Varadhan.
        mode <- floor((n + 1) * p)
        if(x == mode)
            1
        else {
            yL <- if(x < mode) x
                  else {
                      ## Left boundary: largest k <= mode
                      ## with dbinom(k) <= dbinom(x)*relErr 
                      lo <- 0L; hi <- mode
                      while(lo < hi) {
                          mid <- (lo + hi + 1L) %/% 2L
                          if(dbinom(mid, n, p, log = TRUE) > log_thresh)
                              hi <- mid - 1L
                          else
                              lo <- mid
                      }
                      lo
                  }
            yR <- if(x > mode) x
                  else {
                      ## Right boundary: smallest k >= mode
                      ## with dbinom(k) <= dbinom(x)*relErr
                      lo <- mode; hi <- n
                      while(lo < hi) {
                          mid <- (lo + hi) %/% 2L
                          if(dbinom(mid, n, p, log = TRUE) > log_thresh)
                              lo <- mid + 1L
                          else
                              hi <- mid
                      }
                      lo
                  }
            pbinom(yL, n, p) +
                pbinom(yR - 1, n, p, lower.tail = FALSE)
        }
    }
    else {
        ## Do
        ##   d <- dbinom(0 : n, n, p)
        ##   sum(d[d <= dbinom(x, n, p)])
        ## a bit more efficiently ...
        ## Note that we need a little fuzz.
        ## This is tricky: need to be sure
        ## only to sum values in opposite tail
        ## and not count x twice.
        ## For the binomial dist., the mode will
        ## equal the mean if it is an integer.
        m <- n * p
        if (x == m)
            1
        else if (x < m) {
            i <- seq.int(from = ceiling(m), to = n)
            y <- sum(dbinom(i, n, p, log = TRUE) <= log_thresh)
            pbinom(x, n, p) +
                pbinom(n - y, n, p, lower.tail = FALSE)
        } else {
            i <- seq.int(from = 0, to = floor(m))
            y <- sum(dbinom(i, n, p, log = TRUE) <= log_thresh)
            pbinom(y - 1, n, p) +
                pbinom(x - 1, n, p, lower.tail = FALSE)
        }
    }
}