R/beta-fns.R

In DPQ: Density, Probability, Quantile ('DPQ') Computations

##' Qab := Q(a,b) := Gamma(a+b) / Gamma(b)  =  Gamma(a) / Beta(a,b)
##' <==>
##' log(Q(a,b)) = logQab(a,b)
##'             = lgamma(a+b) - lgamma(b) =
##'             = lgamma(a)   - lbeta(a,b)
logQab_asy <- function(a,b, k.max = 5, give.all = FALSE)
{
  ## Purpose:  log(Q(a,b)) asymptotically when max(a,b) -> Inf  &  a^2/b -> C
  ## -------------------------------------------------------------------------
  ## Arguments: a,b: as for  lbeta(.)
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date:  8 Jun 97, 17:28

  if(any(a<=0) || any(b<=0)) stop("a & b must be > 0")
  if(all(a >= b)) { cc <- b; b <- a; a <- cc } #- Now: a <= b
  else if(any(a>b)) stop('Require a <= b  or  a > b   for  ALL vector elements')
  if(length(b)>1) stop("this function only works for SCALAR b (> a)")
  na <- length(a)
  pa <- rbind(Qab_terms(a, k.max)) ##  --> pa = matrix   na  x  (k.max+1)
  ##          ~~~~~~~~~
  ki <- if(k.max >= 1) k.max:1 else numeric(0)
  if(give.all) {
    r <- matrix(0, nrow=na, ncol=k.max+1,
                dimnames = list(format(a), paste(k.max:0)))
    for(k in ki) r[,k] <- pa[,k] + r[,k+1]*(a-k-1)/b
  } else {
    r <- numeric(na)
    for(k in ki) r <- pa[,k] + r*(a-k-1)/b
  }
  a*log(b) + log(1 + a*(a-1)/b * r)
}

Qab_terms <- function(a, k)
{
  ## Purpose: Compute the terms used for  Qab, and beta function
  ## 		Qab := Q(a,b) := Gamma(a+b) / Gamma(b)
  ## -------------------------------------------------------------------------
  ## Arguments: a: the smaller of the arguments of  beta(a,b)
  ##            k: the number of terms in the series expansion
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 19 Jun 97, 11:08

  k <- as.integer(k)
  if(!is.numeric(a)) a <- as.numeric(a) #- leave integers alone..
  if(length(k) != 1 || k > 5 || k < 0)
    stop("'k' must be 1 number in { 0, 1,..,5 }")
  na <- length(a)
  Qab.trms <- expression(
      rep.int(1/2, na),	##  k=1
      a/8 - 1/24, 	##  k=2
      a*(a-1)/48,	##  k=3
      (2/5 + a*(1+3*a*(a-2)))/1152,  ## k=4
      a*(1 + a*(2+a*(a-5)))/5760     ## k=5
  )[seq_len(k)]

  vapply(Qab.trms, eval, numeric(na), envir=environment())
}

### NOTA BENE: All these  lbeta*() are -- I think -- not better than the
##  ---------  builtin  lbeta()  [which exists only since 1998]

lbeta_asy <- function(a,b, k.max = 5, give.all = FALSE)
{
  ## Purpose:  log(beta( a, b))  which also works for HUGE  a / b or b/a
  ## -------------------------------------------------------------------------
  ## Arguments: a,b: as for  lbeta(.)
  ##  log(B(a,b)) = log(G(a) G(b) / G(a+b)) =
  ##              = log(G(a)) - log( G(a+b)/G(b) ) = log(G(a)) - log(Q(a,b))
  lgamma(a) - logQab_asy(a,b, give.all = give.all)
}


### FIXME ==> ../tests/qbeta-dist.R
##            -------- ~~~~~~~~~~~~ at the end has already *several* such approximations !!! <<< FIXME

##' Taylor Approximation to  log(a*B(a,b))  for  a << 1 -- derived from a --> 0
##'
##' for small a,  B(a,b) ~ 1/a, hence  a*B(a,b) ~ a * 1/a = 1, hence log(a*B(a,b)) ~= 0
##'
##' log(a * beta(a,b)) == log(a) + log(beta(a,b)) == log(a) + lbeta(a,b), but the latter suffers cancellation !!
##' We use the Taylor series
##' log(a * B(a,b)) = a \sum_{n=0}^\infty (\psi^{(n)}(1) - \psi^{(n)}(b)) \frac{a^n}{(n+1)!}
##'              where \psi^{(n)} is the n-th derivative of the "psi gamma" function
##' derived from Abramowitz & Stegun, p.259, (6.3.14) series for \psi(1+z)
##' it's integral is \log \Gamma(z+1) = - \gamma z + \sum_{k=2}^\infty \zeta(k)/k (-z)^k   (G)
##' and as \Gamma(z+1) = z \Gamma(z), have \log \Gamma(y) = \log(\Gamma(y+1)/y) = - \log(y) + formula (G) above.
##'
##' Now, \log(a \times B(a,b)) = \log(a) + \log \frac{\Gamma(a) \Gamma(b)}{\Gamma(a+b)}
##'                            = \log(a) + \log \Gamma(a)         + \log (\Gamma(b) / \Gamma(a+b))
##'                            = \gamma a + \pi^2/12 a^2 +O(a^3)  + \log (\Gamma(b) / \Gamma(a+b))
##'  ...
##'                            = a \times \sum_{n=0}^\infty .. sum above

##' @title Accurate Approximation of log(a * beta(a,b)) for small a
##' @param a, b  numeric vectors: the non-negative arguments of Beta
##' @return
##' @author Martin Maechler
laBeta <- function(a,b, nT) {
    stopifnot(length(nT) == 1, nT == as.integer(nT), nT >= 1)

### ===> ../tests/qbeta-dist.R
##       =========~~~~~~~~~~~~ at the end has already *several* such approximations !!! <<< FIXME

}


## MM: this is unused (why ??) -- should be "uniformly better" than lbeta_asy() above
lbetaMM <- function(a,b, cutAsy = 1e-2, verbose = FALSE)
{
  ## Purpose:  log(beta( a, b))  which also works for HUGE  a / b or b/a
  ## -------------------------------------------------------------------------
  ## Arguments: a,b: as for  lbeta(.)
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date:  5 May 97, 17:17
  if(a <= 0 || b <= 0) stop("a & b must be > 0")
  if(length(a)>1 || length(b) > 1)
    stop("this function only works for SCALAR arguments")
  stopifnot(is.numeric(cutAsy), length(cutAsy) == 1, cutAsy >= 0)
  if(a>b) { cc <- b; b <- a; a <- cc } #- Now: a <= b
  if(a*a < b*cutAsy) {
    if(verbose)
      message("a=",formatC(a)," b=",formatC(b), " -- using asymptotic  lbeta(.)")
    lgamma(a) - logQab_asy(a,b)
  } else
    lgamma(a)+lgamma(b)-lgamma(a+b) ## was =: lbeta00(a, b)
}

lbetaM <-  function(a,b, k.max = 5, give.all = FALSE)
{
  ## Purpose:  log(beta( a, b))  which also works for HUGE  a / b or b/a
  ## -------------------------------------------------------------------------
  ## Arguments: a,b: as for  lbeta(.)
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date:  5 May 97, 17:17
  if(any(a<=0) || any(b<=0)) stop("a & b must be > 0")
  na <- length(a)
  nb <- length(b)
  if(nb == 1 || nb == na) {
    r <- numeric(na)
    for(i in 1:na) r[i] <- lbeta_asy(a[i], b[if(nb>1)i else 1],
                                     k.max= k.max, give.all = give.all)
  } else { ##--- return  outer(..)
    r <- matrix(0, nrow = nb, ncol = na)
    for(i in 1:na)r[,i] <- lbeta_asy(a[i], b, k.max= k.max, give.all = give.all)
  }
  r
}
lbetaI <- function(a, n)
{
  ## Purpose:  log(beta(a, n)) for (small) INTEGER n.
  ##    beta(a,n) = G(a) G(n) / G(a+n) = (n-1)! / (a*(a+1)*...*(a+n-1))
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 20 May 97, 10:08
  if(is.na(n <- as.integer(n)) || min(n) <= 0 || max(n) > 10000)
    stop("'n' must be positive (not too large) integer")
  r <- numeric(nn <- length(n))
  for (i in 1:nn)
    r[i] <- lgamma(n[i]) - sum(log(a+0:(n[i]-1)))
  r
}

betaI <- function(a, n)
{
  ## Purpose: beta(a, n) for (small) INTEGER n.
  ##    beta(a,n) = G(a) G(n) / G(a+n) = 1*2*...*(n-1) / (a*(a+1)*...*(a+n-1))
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 20 May 97, 10:08
  if(length(n) > 1) stop("'n' must have length 1.")
  if(is.na(n <- as.integer(n)) || n <= 0)
    stop("'n' must be positive (not too large) integer")
  ni <- seq_len(n-1L)
  prod(ni) / prod(a+c(0,ni))
}

G_half <- 1.7724538509055160272981674833411451827974 # sqrt(pi) = Gamma(1/2)

if(FALSE) ## not closed form smart formula I think (2019-08)
lbetaIhalf <- function(a, n)
{
  ## Purpose:  log(beta(a, n + 1/2)) for (small) INTEGER n.
  ##    beta(a,n + h) = G(a) G(n+h) / G(a+n+h) = (n-1)! / (a*(a+1)*...*(a+n-1))
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 20 May 97, 10:08
  if(length(n) > 1) stop("'n' must have length 1.")
  if(is.na(n <- as.integer(n)) || n <= 0)
    stop("'n' must be positive (not too large) integer")

  stop("unfinished -- FIXME")
}


### pbeta()  --- is now used from TOMS 708 exclusively
### -------  ~/R/D/r-devel/R/src/nmath/toms708.c
###                                    ~~~~~~~~~

## Provide TOMS 708 version of  expm1(), called REXP() in Fortran, then rexpm1() in our C code
## Keeping name '*exmp1()' because, indeed, it should be close to expm1()
##
## EVALUATION OF THE FUNCTION EXP(X) - 1
rexpm1 <- function(x)
{
    ## Vectorize in 'x' -- and work with NaN/NA
    r <- x
    notNA <- !is.na(x)
    sml <- abs(x) <= 0.15

    if(any(s <- sml & notNA))
      r[s] <- # |x| <= 0.15
        local({ x <- x[s]
            ## minimax rational approximation, according to Didonato & Morris (1986)
            ## "Computation of the Incomplete Gamma Function Ratios and their Inverse"
            ## ACM Trans. on Math. Softw. 12, 377--393  // https://dl.acm.org/doi/10.1145/22721.23109
            p1 = 9.14041914819518e-10;
            p2 = .0238082361044469;
            q1 = -.499999999085958;
            q2 = .107141568980644;
            q3 = -.0119041179760821;
            q4 = 5.95130811860248e-4;
            ##
            x * (((p2 * x + p1) * x + 1.) /
                 ((((q4 * x + q3) * x + q2) * x + q1) * x + 1.))
        })
    if(any(B <- !sml & notNA))
      r[B] <- # |x| > 0.15
        local({ x <- x[B]
            w = exp(x)
            ## ifelse(x > 0,
            ##        w * (0.5 - 1./w + 0.5),
            ##        w - 0.5 - 0.5)
            if(any(p <- x > 0))
                w[p] <- w[p] * (0.5 - 1/w[p] + 0.5)
            if(any(n <- !p)) # n = (x <= 0)
                w[n] <- (w[n] - 0.5) - 0.5
            w
        })
    r
} # rexpm1

## Provide TOMS 708 version of  rlog1(), called RLOG1() in Fortran
##
rlog1 <- function(x)
{
## /* -----------------------------------------------------------------------
##  *             Evaluation of the function  x - ln(1 + x)
##  *                                        ~=~ - log1pmx(x)  in ./pgamma.c
##  * ----------------------------------------------------------------------- */
    ## Vectorize in 'x' -- and work with NaN/NA
    r <- x
    notNA <- !is.na(x)
    sml <- -0.39 <= x & x < 0.57 # x inside [-0.39, 0.57)
    if(any(s <- sml & notNA))
      r[s] <- # x in [-0.39, 0.57)
        local({ x <- x[s]
            ## only for "mid": |x| <= 0.18 ; Argument Reduction
            h <- x
            w1 <- rep(0, length(x))
            if (any(L <- x < -0.18)) { ## "Left": x in [-0.39, -0.18)  ("L10")
                h[L] <- h. <- (x[L] + .3) / .7
                a = .0566749439387324
                w1[L] <- a - h. * .3
            }
            if (any(R <- x > 0.18)) { ## "Right": x in (0.18, 0.57) ("L20")
                h[R] <- h. <- x[R] * .75 - .25
                b = .0456512608815524
                w1[R] <- b + h. / 3.
            }
            ##  L30: "Series Expansion" i.e., rational approximation
            r = h / (h + 2.)
            t = r * r

            p0 = .333333333333333;
            p1 = -.224696413112536;
            p2 = .00620886815375787;
            q1 = -1.27408923933623;
            q2 = .354508718369557;
            w <- ((p2 * t + p1) * t + p0) / ((q2 * t + q1) * t + 1.)

            t * 2 * (1/(1 - r) - r * w) + w1
        })

    if(any(B <- !sml & notNA))
      r[B] <- { #  direct evaluation: x - log(1+x)
        x <- x[B]
        x - log((x + 0.5) + 0.5)
      }
    r
} ## rlog1



## From ~/R/D/r-devel/R/src/nmath/pgamma.c :
##                                --------
##/* Compute  log(gamma(a+1))  accurately also for small a (0 < a < 0.5). */

## ---> lgamma1p() below

scalefactor <- 2^256 # == 2^2^2^3 == double prec. exact{checking via Rmpfr and gmp::as.bigz(2)^256}

## smallest exponent E for which exp(E) is "above underflow" {exp(E) == .Machine$double.xmin}
M_minExp <- log(2) * .Machine$double.min.exp # ~= -708.396

##/* If |x| > |k| * M_cutoff,  then  log[ exp(-x) * k^x ] =~=  -x */
M_cutoff <- log(2) * .Machine$double.max.exp / .Machine$double.eps ## = 3.196577e18


##' a useful 1-string  format(summary(<logical>)): __TODO: add to {sfsmisc} __
formatSummL <- function(x) {
    stopifnot(is.logical(x))
    if(!length(x)) return("logi(0)")
    ## else non-empty -- potentially has  {TRUE, FALSE, NA}
    tab <- table(x, exclude=NULL)
    if(any(I <- is.na(names(tab))))# "NA", not <NA>
        names(tab)[I] <- "NA"
    nt <- sort(names(tab), decreasing = TRUE)# typical starting with "TRUE"
    k <- length(
        snt <- if(any(x & !is.na(x))) sub("FALSE", "F.", nt, fixed=TRUE) else nt)
    m <- rbind(unname(tab[nt]), " ", snt
             , if(k >= 2) c(rep(", ", k-2), ", and ", "") # else NULL
             , deparse.level=0L)
    paste(m, collapse="")
}

## This is the "fast pure R" version that does iterations
## in parallel *where applicable*
## The next one, logcfR(), will entirely vectorize, ==> iterations separately
logcfR. <- function (x, i, d, eps, # ~ relative tolerance
                     maxit = 10000L, trace = FALSE)
{
    ## Continued fraction for calculation of  sum_{k=0}^Inf x^k/(i+k*d) =
    ##		1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ...
    ##
    ## auxiliary in log1pmx() and lgamma1p()
    stopifnot (i > 0, d >= 0, eps > 0,
               length(i) == 1, length(d) == 1, length(eps) == 1
               ## x < 1 for (i,d = 3,2) -- there x = -2 is fine, but in general?
               ## |x| < 1 # otherwise it does not (always?) converge,
               ##            ??
               )
    if(any_mpfr(x, i, d)) {
        stopifnot(requireNamespace("Rmpfr"))
        mpfr <- Rmpfr::mpfr ; getPrec <- Rmpfr::getPrec
        prec <- max(getPrec(x), getPrec(i), getPrec(d))
        mean <- Rmpfr::mean
        ## this *is* necessary / improving results!
        if(!inherits(x, "mpfr")) x <- mpfr(x, prec)
        if(!inherits(i, "mpfr")) i <- mpfr(i, prec)
        if(!inherits(d, "mpfr")) d <- mpfr(d, prec)
        asN <- Rmpfr::asNumeric
    } else asN <- as.numeric

    c1 <- 2 * d  # scalar
    c2 <- i + d  #  "
    c4 <- c2 + d #  "
    ## these will all depend on x :
    a1 <- c2
    b1 <- i * (c2 - i * x)
    r <- d * d * x
    A2 <- c4 * c2 - r
    B2 <- c4 * b1 - i * r

    it <- 0L
    if(trace) cat(sprintf("logcf(x[], i=%g, d=%g, eps=%g)", i, d, eps),
                  "  iterations:\n")
    n <- length(r <- x) # r[1:n] : the result (also good for 'mpfr' numbers)
    do.x <- rep(TRUE, n)
    while (any(needIt <- abs(A2 * b1 - a1 * B2) > abs(eps * b1 * B2))) {
        ## needIt, x, (a1,b1, A2, B2) are all of "shrinking length"
        ## (r, do.x) must stay at length == n
        it <- it+1L
        if(trace)
            cat(sprintf("it=%2d: needIt: %19s; ", it, formatSummL(needIt)))
        if(length(iFini <- which(!needIt))) {
            ## compute iFi, the indices in the *original* length-n vectors
            iFi <- which(do.x)[iFini]
            r   [iFi] <- A2[iFini] / B2[iFini]
            do.x[iFi] <- FALSE
            ## and shorten the rest:
            iK <- which(needIt) # the indices to keep
             x <-  x[iK]
            if(length(a1) == length(b1)) # (not when it = 1)
            a1 <- a1[iK]
            b1 <- b1[iK]
            A2 <- A2[iK]
            B2 <- B2[iK]
        }
        ## else : all(needIt) --> no finished parts, no shortening

        if(trace) cat(sprintf("length(x[<todo>])=%2d, ", sum(do.x)))

        c3 <- c2*c2 * x
	c2 <- c2 + d
	c4 <- c4 + d
	a1 <- c4 * A2 - c3 * a1
	b1 <- c4 * B2 - c3 * b1

	c3 <- c1*c1 * x
	c1 <- c1 + d
	c4 <- c4 + d
	A2 <- c4 * a1 - c3 * A2
	B2 <- c4 * b1 - c3 * B2

        m.B2 <- exp(mean(log(abs(B2))))
        if(trace) cat(sprintf(" m.B2=%12g", asN(m.B2)))
	if (m.B2 > scalefactor) {
            if(trace) cat(sprintf("  Lrg m.B2"))
	    a1 <- a1 / scalefactor
	    b1 <- b1 / scalefactor
	    A2 <- A2 / scalefactor
	    B2 <- B2 / scalefactor
	}
        else if (m.B2 < 1 / scalefactor) {
            if(trace) cat(sprintf("  Sml m.B2"))
	    a1 <- a1 * scalefactor
	    b1 <- b1 * scalefactor
	    A2 <- A2 * scalefactor
	    B2 <- B2 * scalefactor
	}
        if(trace) cat("\n")
        if(it > maxit) {
            warning("non-convergence in logcf(), iter > ", maxit)
            break
        }
    }
    if(trace && it <= maxit) cat("  logcf(*) end: after", it, "iterations.\n")

    r[do.x] <- A2 / B2
    ## and return
    r
} ## logcfR.()

## This is the entirely vectorizing "pure R" version
logcfR <- function (x, i, d, eps, # ~ relative tolerance
                    maxit = 10000L, trace = FALSE)
{
    ## Continued fraction for calculation of  sum_{k=0}^Inf x^k/(i+k*d) =
    ##		1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ...
    ##
    ## auxiliary in log1pmx() and lgamma1p()
    stopifnot (i > 0, d >= 0, eps > 0,
               length(i) == 1, length(d) == 1, length(eps) == 1
               ## x < 1 for (i,d = 3,2) -- there x = -2 is fine, but in general?
               ## |x| < 1 # otherwise it does not (always?) converge,
               ##            ??
               )
    if(any_mpfr(x, i, d)) {
        stopifnot(requireNamespace("Rmpfr"))
        mpfr <- Rmpfr::mpfr ; getPrec <- Rmpfr::getPrec
        prec <- max(getPrec(x), getPrec(i), getPrec(d))
        ## this *is* necessary / improving results!
        if(!inherits(x, "mpfr")) x <- mpfr(x, prec)
        if(!inherits(i, "mpfr")) i <- mpfr(i, prec)
        if(!inherits(d, "mpfr")) d <- mpfr(d, prec)
        asN <- Rmpfr::asNumeric
    } else asN <- as.numeric

    r <- x
    for(ii in seq_along(x)) {
        xi <- x[ii]
        c1 <- 2 * d
        c2 <- i + d
        c4 <- c2 + d
        a1 <- c2
        b1 <- i * (c2 - i * xi)
        ddx <- d * d * xi
        A2 <- c4 * c2 - ddx
        B2 <- c4 * b1 - i * ddx

        it <- 0L
        if(trace) cat(sprintf("logcf(x[%2d]=%8g, i=%g, d=%g, eps=%g)", ii, xi, i, d, eps),
                      if(trace >= 2)"  iterations:\n")
        while (abs(A2 * b1 - a1 * B2) > eps * abs(b1 * B2)) {
            c3 <- c2*c2*xi
            c2 <- c2 + d
            c4 <- c4 + d
            a1 <- c4 * A2 - c3 * a1
            b1 <- c4 * B2 - c3 * b1

            c3 <- c1 * c1 * xi
            c1 <- c1 + d
            c4 <- c4 + d
            A2 <- c4 * a1 - c3 * A2
            B2 <- c4 * b1 - c3 * B2

            it <- it+1L
            if(trace >= 2) cat(sprintf("it=%2d: ==> B2=%13g", it, asN(B2)))
            if (abs(B2) > scalefactor) {
                if(trace >= 2) cat(sprintf("  Lrg m.B2\n%27s", ""))
                a1 <- a1 / scalefactor
                b1 <- b1 / scalefactor
                A2 <- A2 / scalefactor
                B2 <- B2 / scalefactor
            }
            else if (abs(B2) < 1 / scalefactor) {
                if(trace >= 2) cat(sprintf("  Sml m.B2\n%27s", ""))
                a1 <- a1 * scalefactor
                b1 <- b1 * scalefactor
                A2 <- A2 * scalefactor
                B2 <- B2 * scalefactor
            }
            if(trace >= 2) cat(sprintf(" --> crit. |A2*b1 - a1*B2|/|b1*B2| = %g\n",
                                  asN(abs(A2 * b1 - a1 * B2)/abs(b1 * B2))))
            if(it > maxit) {
                warning("non-convergence in logcf(), iter > ", maxit)
                break
            }
        }
        if(trace && it <= maxit) cat(sprintf("  logcf(*) end: after %3d iterations.\n", it))
        r[ii] <- A2 / B2
        ##       -------
    } ## end for(ii ..)
    r
} ## logcfR() {"properly" vectorized}


logcf <- function (x, i, d, eps, trace = FALSE) {
    .Call(C_R_logcf, x, i, d, eps, trace)
}

## Accurate calculation of log(1+x)-x, particularly for small x.
## See also R-interface to R's C API [src/nmath/pgamma.c] log1pmx()
##     --> log1pmxC() in >> ./utils.R >> ../src/DPQ-misc.c
log1pmx <- function(x, tol_logcf = 1e-14,
                    eps2 = 0.01,
                    minL1 = -0.79149064, ## << was hard-wired 'minLog1Value' in R's source of log1pmx()
                    trace.lcf = FALSE,
                    logCF = if(is.numeric(x)) logcf else logcfR.)
{
    stopifnot(is.numeric(eps2), eps2 >= 0, is.numeric(minL1), -1 <= minL1, minL1 < 0)# < -1/4 ?
    r <- x
    if(any(c1 <- (x > 1 | x < minL1)))
        r[c1] <- log1p(x[c1]) - x[c1]
    ## else { ## ##/* expand in [x/(2+x)]^2 */
    if(any(c2 <- !c1)) {
        x <- x[c2]
	term <- x / (2 + x)
	y <- term * term
        r[c2] <- term * { ## not using ifelse(), rather what works with "mpfr"
            ## ifelse(abs(x) < eps2,
            ##        (((2 / 9 * y + 2 / 7) * y + 2 / 5) * y + 2 / 3) * y - x,
            ##        2 * y * logcf(y, 3, 2, tol_logcf) - x)
            A <- x
            if(any(isSml <- abs(x) <= eps2)) {
                i <- which(isSml)
                y. <- y[i]
                z <- if(is.numeric(x)) 2 else 2 + 0*y. # becomes y-like (e.g. mpfr)
                A[i] <- (((z / 9 * y. + z / 7) * y. + z / 5) * y. + z / 3) * y. - x[i]
            }
            if(length(iLrg <- which(!isSml))) {
                y. <- y[iLrg]
                A[iLrg] <- 2 * y. *
                    logCF(y., 3, 2, tol_logcf, trace = trace.lcf) - x[iLrg]
            }
            A
        }
    }
    r
}

lgamma1p <- function(a, tol_logcf = 1e-14, f.tol = 1., ...)
{
    ## Compute  log(gamma(a+1))  accurately also for small a (0 < a < 0.5).

    eulers_const <- 0.5772156649015328606065120900824024

    ##/* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */
    N <- 40
    coeffs <-
        c(0.3224670334241132182362075833230126e-0, ##/* = (zeta(2)-1)/2 */
          0.6735230105319809513324605383715000e-1, ##/* = (zeta(3)-1)/3 */
          0.2058080842778454787900092413529198e-1,
          0.7385551028673985266273097291406834e-2,
          0.2890510330741523285752988298486755e-2,
          0.1192753911703260977113935692828109e-2,
          0.5096695247430424223356548135815582e-3,
          0.2231547584535793797614188036013401e-3,
          0.9945751278180853371459589003190170e-4,
          0.4492623673813314170020750240635786e-4,
          0.2050721277567069155316650397830591e-4,
          0.9439488275268395903987425104415055e-5,
          0.4374866789907487804181793223952411e-5,
          0.2039215753801366236781900709670839e-5,
          0.9551412130407419832857179772951265e-6,
          0.4492469198764566043294290331193655e-6,
          0.2120718480555466586923135901077628e-6,
          0.1004322482396809960872083050053344e-6,
          0.4769810169363980565760193417246730e-7,
          0.2271109460894316491031998116062124e-7,
          0.1083865921489695409107491757968159e-7,
          0.5183475041970046655121248647057669e-8,
          0.2483674543802478317185008663991718e-8,
          0.1192140140586091207442548202774640e-8,
          0.5731367241678862013330194857961011e-9,
          0.2759522885124233145178149692816341e-9,
          0.1330476437424448948149715720858008e-9,
          0.6422964563838100022082448087644648e-10,
          0.3104424774732227276239215783404066e-10,
          0.1502138408075414217093301048780668e-10,
          0.7275974480239079662504549924814047e-11,
          0.3527742476575915083615072228655483e-11,
          0.1711991790559617908601084114443031e-11,
          0.8315385841420284819798357793954418e-12,
          0.4042200525289440065536008957032895e-12,
          0.1966475631096616490411045679010286e-12,
          0.9573630387838555763782200936508615e-13,
          0.4664076026428374224576492565974577e-13,
          0.2273736960065972320633279596737272e-13,
          0.1109139947083452201658320007192334e-13 ##/* = (zeta(40+1)-1)/(40+1) */
          )

    c <- 0.2273736845824652515226821577978691e-12##/* zeta(N+2)-1 */

    r <- a
    ## if(abs(a) >= 0.5)
    a.lrg <- (abs(a) >= 0.5)
    r[a.lrg] <- lgamma(a[a.lrg] + 1)
    ## else {
    if(any(a.sml <- !a.lrg)) {
        a. <- a[a.sml]
        ## Abramowitz & Stegun 6.1.33 : for |x| < 2,
        ## <==> log(gamma(1+x)) = -(log(1+x) - x) - gamma*x + x^2 * sum_{n=0..Inf} c_n (-x)^n
        ## where c_n := (Zeta(n+2) - 1)/(n+2)  = coeffs[n]
        ##
        ## Continued fraction convergence acceleration is used to compute
        ## lgam(x) :=  sum_{n=0..Inf} c_n (-x)^n
        lgam <- c * logcf(-a. / 2, N + 2, 1, tol_logcf)
        for (i in N:1)
            lgam <- coeffs[i] - a. * lgam

        ## return
        r[a.sml] <- (a. * lgam - eulers_const) * a. -
            log1pmx(a., tol_logcf = f.tol * tol_logcf, ...)
    }
    r
}## lgamma1p


lgamma1p_series <- function(x, k) {
    stopifnot(k == as.integer(k), 1 <= k, k <= 11)
    ## From Maple : ~/maple/gamma-asympt.txt

    ##  lG := series(ln(GAMMA(x+1)), x, 11);
    ##                   1    2  2   1          3    1    4  4
    ##  lG := -gamma x + -- Pi  x  - - Zeta(3) x  + --- Pi  x
    ##                   12          3              360

    ##       1          5    1     6  6   1          7     1     8  8
    ##     - - Zeta(5) x  + ---- Pi  x  - - Zeta(7) x  + ----- Pi  x
    ##       5              5670          7              75600

    ##       1          9     1      10  10    / 11\
    ##     - - Zeta(9) x  + ------ Pi   x   + O\x  /
    ##       9              935550

    ## Maple> hG := convert(convert(lG, polynom), horner);
    ## gives
    ## (-gamma +
    ##   (Pi*Pi/12 +
    ##    (-Zeta(3)/3+
    ##     (Pi*Pi*Pi*Pi/360+
    ##      (-Zeta(5)/5+
    ##       (Pi*Pi*Pi*Pi*Pi*Pi/5670+
    ##        (-Zeta(7)/7+
    ##         (Pi*Pi*Pi*Pi*Pi*Pi*Pi*Pi/75600+
    ##          (-Zeta(9)/9 +
    ##           Pi^10/935550 * x) * x) * x) * x) * x) * x) * x) * x) * x) * x

    useM <- (inherits(x, "mpfr"))
    if(useM) {
        stopifnot(requireNamespace("Rmpfr"))
        prec <- max(Rmpfr::getPrec(x))
        mpfr <- Rmpfr::mpfr
        zeta <- Rmpfr::zeta
    }
    ## "Euler's constant" gamma
    gamma <- if(useM) Rmpfr::Const("gamma", prec) else 0.57721566490153286
    if(useM) pi <- Rmpfr::Const("pi", prec) # else use R's pi # 3.1415926535897932

    if(k >= 2) {
        px <- (pi*x)^2 # = pi^2 x^2
        px4 <- px/4    # = pi^2 x^2 / 4
        if(k >= 3) {
            z3 <- if(useM) zeta(mpfr(3, prec)) else 1.202056903159594285
            x2 <- x*x
            if(k >= 5) {
                z5 <- if(useM) zeta(mpfr(5, prec)) else 1.036927755143369927
                if(k >= 6) {
                    I <- if(useM) mpfr(1, prec) else 1 # use I/3 = 1/3 in high precision
                    if(k >= 7) {
                        z7 <- if(useM) zeta(mpfr(7, prec)) else 1.008349277381922827
                        ## if(k >= 9) {
                        ##     z9 <- if(useM) zeta(mpfr(9, prec)) else 1.002008392826082214
                        ##     if(k >= 11) {
                        ##         z11 <- if(useM) zeta(mpfr(11, prec)) else 1.000494188604119464
                        ##     }
                        ## }

                    }
                }
            }
        }
    }
   ## return
    switch(k,
           -gamma * x,
           -gamma * x + px4/3,    # k = 2
           -gamma * x + px4/3 - x*x2/3*z3,    # k = 3
           -gamma * x - x*x2*(z3/3) + px4/3*(1 + px4/7.5),    # k = 4
           ##   pi^4/360 = (pi^2/4)^2 / 3 / 7.5  (3 * 16 * 7.5 == 360)
           ## k = 5 :
           -gamma * x - x*x2*(z3/3 + x2*(z5/5)) + px4/3*(1 + px4/7.5),
           ## k = 6 : (fractions: multiple of 2^-k are exact also in double prec)
           ##         360/16 = 22.5   5670/16 = 354.375
           -gamma * x - x*x2*(z3/3 + x2*(z5/5))            + px4*(1/3 + px4*(I/22.5 + px/354.375)),
           -gamma * x - x*x2*(z3/3 + x2*(z5/5+ x2*(z7/7))) + px4*(1/3 + px4*(I/22.5 + px/354.375)), # k = 7
           stop("Currently, only k <= 7  is implemented, but k=",k))
}

algdiv <- function(a,b) .Call(C_R_algdiv, as.double(a), as.double(b))

bpser <- function(a,b, x, log.p=FALSE, eps = 1e-15, verbose=FALSE, warn=TRUE) { # ../src/bpser.c
    .Call(C_R_bpser, a, b, x, eps, log.p, verbose, warn)
    ## return( list(r = r_, err = ier_) )
}



qnormAppr <- function(p) {
##' original (getting inaccurate when p --> 1)
  ## qnorm: normal quantile approximation for qnorm(p),  for p > 1/2
  ## -- to be used in  qbeta(.)
  ## The relative error of this approximation is quite ASYMMETRIC: mainly < 0
  ## NB:  This is Abramowitz & Stegun 26.2.22
    .Deprecated("qnormUappr")
  a0 <- 2.30753
  a1 <- 0.27061
  b1 <- 0.99229
  b2 <- 0.04481
  t <- sqrt(-2*log(1-p))
  t - (a0 + a1 * t) / (1 + (b1 + b2 * t) * t)
}

qnormUappr <- function(p,
                       lp = .DT_Clog(p, lower.tail=lower.tail, log.p=log.p),
                                        # ~= log(1-p) -- independent of lower.tail, log.p
                       lower.tail=FALSE, log.p = missing(p), tLarge = 1e10)
{
    ## qnorm: normal quantile approximation;
    ## -- to be used in  qbeta(.)
    ## The relative error of this approximation is quite ASYMMETRIC: mainly < 0
    ## NB:  This is Abramowitz & Stegun 26.2.22; improved with swapping etc
    a0 <- 2.30753
    a1 <- 0.27061
    b1 <- 0.99229
    b2 <- 0.04481

    if(use.lp <- missing(p)) { ## log.p unused;  lp = log(1-p)  <==>  e^lp = 1-p  <==>  p = 1 - e^lp
        if(lower.tail) stop("lower.tail=TRUE not allowed when p is not specified, but lp is")
        p. <- -expm1(lp)
        ## swap <- TRUE -- will swap *all* below
    } else {
        p. <- .D_qIv(p, log.p)
        ## swap p <--> 1-p -- so we are where approximation is better
        swap <- if(lower.tail) p. < 1/2 else p. > 1/2 # logical vector
        p[swap] <- if(log.p) log1mexp(-p[swap]) else 1 - p[swap]
    }
    R <- t <- sqrt(-2 * lp)
    t.ok <- t < tLarge ## e.g., for p == 1, t = Inf would give R = NaN
    ## for t >= tLarge:  R = t numerically in formula below
    t <- t[t.ok]
    R[t.ok] <- t - (a0 + a1 * t) / (1 + (b1 + b2 * t) * t)
    if(use.lp) R <- -R else R[swap] <- -R[swap]
    R[p. == 1/2] <- 0
    R
}

qnormUappr6 <- function(p,
                        lp = .DT_Clog(p, lower.tail=lower.tail, log.p=log.p),
                                        # ~= log(1-p) -- independent of lower.tail, log.p
                        lower.tail=FALSE, log.p = missing(p), tLarge = 1e10)
{
    ## qnorm: normal quantile approximation;
    ## ....
    ## NB:  This is Abramowitz & Stegun 26.2.23 (6 coefficients)
    ##                                  ~~~~~~~  =
    c0 <- 2.515517;  d1 <- 1.432788
    c1 <-  .802853;  d2 <-  .189269
    c2 <-  .010328;  d3 <-  .001308

    if(use.lp <- missing(p)) { ## log.p unused;  lp = log(1-p)  <==>  e^lp = 1-p  <==>  p = 1 - e^lp
        if(lower.tail) stop("lower.tail=TRUE not allowed when p is not specified, but lp is")
        p. <- -expm1(lp)
        ## swap <- TRUE -- will swap *all* below
    } else {
        p. <- .D_qIv(p, log.p)
        ## swap p <--> 1-p -- so we are where approximation is better
        swap <- if(lower.tail) p. < 1/2 else p. > 1/2 # logical vector
        p[swap] <- if(log.p) log1mexp(-p[swap]) else 1 - p[swap]
    }
    R <- t <- sqrt(-2 * lp)
    ##          vvvvv <=== check which one is optimal !!  (use Rmpfr qnormR()}
    t.ok <- t < tLarge ## e.g., for p == 1, t = Inf would give R = NaN
    ## for t >= tLarge:  R = t numerically in formula below
    t <- t[t.ok]
    R[t.ok] <- t - (c0 + (c1 + c2*t) * t) / (1 + (d1 + (d2 + d3*t) * t) * t)
    if(use.lp) R <- -R else R[swap] <- -R[swap]
    R[p. == 1/2] <- 0
    R
}


##              --------
##===> more:  ./norm_f.R <<<<
##              --------


qbetaAppr.1 <- function(a, p, q, lower.tail=TRUE, log.p=FALSE,
                        y = qnormUappr(a, lower.tail=lower.tail, log.p=log.p))
{
  ## Purpose: Approximate  qbeta(a, p,q) -- Abramowitz & Stegun (26.5.22)
  ##          qbeta(.) takes this only when  p>1 & q>1 : "Carter(1947), see AS 109, remark '5.'"
  ## -------------------------------------------------------------------------
  ## Arguments: a: percentage point, only used via qnorm*(a,..);  (p,q): beta parameters

  r <- (y * y - 3) / 6 ## lambda
  s <- 1 / (p + p - 1)
  t <- 1 / (q + q - 1)
  h <- 2 / (s + t)
  w <- y * sqrt(h + r) / h - (t - s) * (r + 5 / 6 - 2 / (3 * h))
  p / (p + q * exp(w + w))
}

qbetaAppr.3 <- function(a, p, q, lower.tail=TRUE, log.p=FALSE, logbeta = lbeta(p,q))
{
   ##  a=alpha
    ## Purpose: Approximate  qbeta(a, p,q) -- for small a [log.p=TRUE: a far negative]
    ##  Inversion of   I_x(a,b) ~= x^a / (a B(a,b))

    ## CARE:  log(p) + logbeta == log(p * beta(p,q)) suffers from cancellation
    ## ----   for small p
    ## ---> look at Qab() above or *rather* also see experiments c12pBeta() , p.err.pBeta()
    ## in  ../tests/qbeta-dist.R
    ##              ============

    ## log.a = log("a") = log(pr_{lower_tail})
    log.a <- if(lower.tail) .D_log(a, log.p=log.p) else .D_LExp(a, log.p=log.p)

    pmin(1, exp(log.a +  log(p) + logbeta) / p)
}

qbetaAppr.2 <- function(a, p, q, lower.tail=TRUE, log.p=FALSE, logbeta = lbeta(p,q))
{
    ## Purpose: Approximate  qbeta(a, p,q)
    ## log(1-"a") = log(pr_{upper_tail}):
    l1ma <- if(lower.tail) .D_LExp(a, log.p=log.p) else .D_log(a, log.p=log.p)
    ## pmax(0, ....) -- needed when  (l1ma + log(q) + logbeta) / q > 0, e.g., for
    ## e.g. qbetaAppr.2(1/4, 1/2, 1/2) gives -0.388
    -expm1((l1ma + log(q) + logbeta) / q)
}

qbetaAppr.4 <- function(a, p, q, lower.tail=TRUE, log.p=FALSE,
                        y = qnormUappr(a, lower.tail=lower.tail, log.p=log.p),
                        verbose=getOption("verbose"))
{
    ## Purpose: Approximate  qbeta(a, p,q); 'a' is only used via qnorm(a,..)
    r <- q + q
    t <- 1 / (9 * q)
    t <- r * (1 - t + y * sqrt(t))^3 # = \chi^2_alpha .. which "must be > 0", but I find: no!
    if(verbose) cat("chi^2[alpha] =: t = ",format(t), if(t <= 0)"t < 0 !!","\n")
    t <- (4 * p + r - 2) / t
    if(verbose) cat("t = (1 + x0)/(1 - x0) = ",format(t), if(t <= 1) "t <= 1 !!","\n")
    1 - 2 / (t + 1)
}

qbetaAppr.5 <- function(a, p, q, lower.tail=TRUE, log.p=FALSE,
                        y = qnormUappr(a, lower.tail=lower.tail, log.p=log.p),
                        logbeta = lbeta(p,q),
                        verbose=getOption("verbose"))
{
    ## Purpose: Approximate  qbeta(a, p,q) for p << q , really from  p --> 0 expansion
    ##
    ## based on future good  laBeta() ==robust== log(p * B(p,q))
    ##
    ## see MM's notes (on paper!)
}


qbetaAppr <- function(a, p, q, lower.tail=TRUE, log.p=FALSE,
                      y = qnormUappr(a, lower.tail=lower.tail, log.p=log.p),
                      logbeta = lbeta(p,q),
                      verbose = getOption("verbose") && length(a) == 1)
{
    ## Purpose: Approximate  qbeta(a, p,q) --- for  a <= 1/2
    ## -------------------------------------------------------------------------
    ## Arguments:
    ## -------------------------------------------------------------------------
    ## Author: After qbeta.c in  R  (--> AS 109...)
    if(length(p)!=1 || length(q)!=1)
        stop("'p' & 'q' must have length 1 !!")
    if(any(a > 1/2)) warning("a[.] > 1/2 (not thought for!)")

    if (p > 1 && q > 1) { ## Abramowitz & Stegun(26.5.22)
        if(verbose) cat("p,q > 1: Using qbetaAppr.1(): Abramowitz-Stegun (26.5.22)\n")
        ## 'a' is not needed here, just 'y' is
        qbetaAppr.1(,p,q, lower.tail=lower.tail, log.p=log.p, y=y)
    } else {
        if(verbose) cat("p or q <= 1 : ")
        r <- q + q
        ## more "robust" version of   t <- 1 / (9 * q) ; t <- r * (1 - t + y * sqrt(t))^3
        st <- 1/(3 * sqrt(q))
        t <- r * (1 - st*(st + y))^3

        ## NOTE: length(t) == length(y) == length(a)  is fulfilled

        ## vectorized in t { ==> logical vectors instead of simple if(..) .. else if(..) .. else ..
        ans <- t
        if(any(neg <- t <= 0)) { ## forget y(a) and t(q, y) from above; qbetaAppr.2():
            if(verbose)   cat("t <= 0:  using qbetaAppr.2()\n")
            ## log(1-"a") = log(pr_{upper_tail}):
            l1ma <- if(lower.tail) .D_LExp(a[neg], log.p=log.p) else .D_log(a[neg], log.p=log.p)
            ## FIXME -- still ? -- : "catastrophic"  in -expm1(xx) when xx > 0 ==> -expm1(xx) = 1-e^xx < 0 "nonsense"
            ans[neg] <- - expm1((l1ma + log(q) + logbeta) / q)
        }
        if(any(!neg)) { ##else
            t <- (4 * p + r - 2) / t
            if (any(L1 <- !neg & t <= 1)) {
                if(verbose) cat("t > 0; t' <= 1: qbetaAppr.3()\n")
                ## ans[L1] <- exp((log(a[L1] * p) + logbeta) / p)
                ans[L1] <- qbetaAppr.3(a[L1], p, q, lower.tail=lower.tail, log.p=log.p, logbeta=logbeta)
                ## ==> linear relationship on log-log scale:
                ##   log(ans) = 1/p * log(a) + (log p + log beta(p,q))/p
            }
            if (any(L2 <- !neg & t > 1)) { # the other case
                if(verbose) cat("t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1)\n")
                ans[L2] <- 1 - 2 / (t[L2] + 1)
            }
        }
        ans
    }
}

## Hmm, this is (sometimes!) not quite equivalent to C's :
## ex.  qbeta.R(0.5078, .01, 5) -> 2.77558e-15 but qbeta() now correctly gives 4.651e-31
qbeta.R	 <-  function(alpha, p, q,
                      lower.tail = TRUE, log.p = FALSE,
		      logbeta = lbeta(p,q),
		      low.bnd = 3e-308, up.bnd = 1-2.22e-16,
                      method = c("AS109", "Newton-log"),
                      tol.outer = 1e-15,
                      ## FIXME: (a,p,q) : and then uses (a, pp) .. hmm
		      f.acu = function(a,p,q) max(1e-300, 10^(-13- 2.5/pp^2 - .5/a^2)),
		      fpu = .Machine$ double.xmin,
		      qnormU.fun = function(u, lu) qnormUappr(p=u, lp=lu)
		    , R.pre.2014 = FALSE
		    , verbose = getOption("verbose") ## FALSE, TRUE, or 0, 1, 2, ..
		    , non.finite.report = verbose
		      )
{
### ----- following the C code in  ~/R/D/r-devel/R/src/nmath/qbeta.c
###
### ----------> make use of 'log.p' in	qnormUappr(* , log.p)	as well!

    ## Purpose:	 qbeta(.) in R -- edited;  verbose output
    ## -------------------------------------------------------------------------
    ## Arguments: alpha: percentage;  (p,q) Beta parameter
    ##	  low.bnd, up.bnd: algorithm parameter
    ##	acu, fpu:	accuracy;     .Machine$ double.xmin = 2.22e-308 for IEEE
    ##		SAE below is the most negative decimal exponent which does not
    ##		cause an underflow; a value of -308 or thereabouts will often be
    ##		OK in double precision.
    ##		sae <- -37 ;  fpu <- 10 ^ sae
    ##	  qnormU.fun(p) ~= qnorm(1 - p) = qnorm(p, ..., lower.tail=FALSE)
    ## -------------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 30 Apr 1997, 12:55
    ## ------ edited C-code  qbeta.c  into  R-code

    if(length(alpha)>1) stop("only works for length(1) 'alpha'")

    if (low.bnd >= up.bnd || low.bnd < 0 || up.bnd > 1)
	stop("MUST  0 <= low.bnd < up.bnd <= 1")

    ## Test for admissibility of parameters
    if (p < 0 || q < 0 || alpha < 0. || alpha > 1.) stop("DOMAIN_ERROR")

    ##  R_Q_P01_boundaries(alpha, 0, 1) :=
    if(log.p) {
        if(alpha == -Inf) return(if(lower.tail) 1 else 0)
        if(alpha ==   0 ) return(if(lower.tail) 0 else 1)
    } else {
        if (alpha == 0 || alpha == 1) return(if(lower.tail) alpha else 1-alpha)
    }

    p. <- .DT_qIv(alpha, lower.tail, log.p=log.p) # = lower_tail prob (in any case)

    if(R.pre.2014) {
        if(log.p && (p. == 0. || p. == 1.))
	return(p.) ## better than NaN or infinite loop;

    } else { ## can and *must* do better:
        if(log.p && (p. == 0. || p. == 1.))
            warning("p. is 0 or 1")
    }

    ## inflection point of pbeta() (outside [0,1] when p < 1 < q or q < 1 < p)
    ## x.ip <- (p-1) / (p+q -2)
    ## y.ip <- pbeta(x.ip, p, q)

    ## change tail if necessary, such that  beta(a,p,q)	 with  0 = a <= 1/2
    ## la := log(a), but without numerical cancellation:
    if (alpha <= 1/2) {
	a <- p.
	la <- if(lower.tail) .D_log(alpha, log.p) else .D_LExp(alpha, log.p)
        pp <- p; qq <- q; swap.tail <- FALSE
    } else {
	a <- .DT_CIv(alpha, lower.tail, log.p)
	la <- if(lower.tail) .D_LExp(alpha, log.p) else .D_log(alpha, log.p)
        pp <- q; qq <- p; swap.tail <- TRUE
    }
    if(verbose) cat(if(swap.tail)"SWAP tail" else "no swap", "\n")
    stopifnot( all.equal(log(a), la) )

    acu <- f.acu(a,p,q)
    ## set the exponent of acu to -2r-2 for r digits of accuracy
    if(acu <= 0 || acu > .1) acu   <- 1.0e-32 # 32: 15 digits accu.

    if(verbose) cat("logbeta=", formatC(logbeta,digits=19),
		    "; acu=",formatC(acu, digits=6) ,"\n")

    ## calculate the initial approximation 'xinbta' [FIXME: never designed for log.p=TRUE]
    xinbta <- qbetaAppr(a, pp,qq, y = qnormU.fun(u=a, lu=la), logbeta=logbeta,
                         verbose = verbose)
    xinbta0 <- xinbta # in case we want get back to it
    if(verbose) cat("initial 'xinbta' =", formatC(xinbta0, digits=18))
    if(!is.finite(xinbta)) {
	warning("** qbeta.R(): qbetaAppr() was not finite!   xinbta := 0.5")
        xinbta <- 0.5
    }

    ## Solve for x by a modified Newton-Raphson method,
    ## using the function  pbeta(.)

    r <- 1 - pp
    t <- 1 - qq
    yprev <- 0
    adj <- 1
    ##not used: prev <- 1
    if (xinbta < low.bnd) { xinbta <- 0.5; if(verbose)cat(" -- low.bnd  -> xi= 0.5\n") } else
    if (xinbta >  up.bnd) { xinbta <- 0.5; if(verbose)cat(" -- up.bnd   -> xi= 0.5\n") } else
    if(verbose) cat("\n")

    ##-- for single precision -- from CORRECTED	 AS 109:
    ##R	 iex <- as.integer( max(-5/pp^2 - 1/a^2 - 13, sae))
    ##NN acu <- max(fpu, 10 ^ (-5/pp^2 - 1/a^2 - 13))

    finish <- FALSE
    o.it <- 0
    method <- match.arg(method)

    switch(method,
           "AS109" =
    while(!finish) { ##-- Outer Loop
	y <- pbeta(xinbta, pp, qq)
	if(verbose) cat("y=pbeta(x..)=", formatC(y, digits=15, width=18, flag='-'))
	if(!is.finite(y)) {
            if(non.finite.report)
	    cat("  y=pbeta(): not finite: 'DOMAIN_ERROR'\n",
		"  xinbta =", format01prec(xinbta, digits=18),"\n")
	    return(xinbta) ##browser()
	}
	y <- (y - a) * exp(logbeta + r * log(xinbta) + t * log1p(-xinbta))
	if(verbose) cat(" y.n=(y-a)*e^..=", formatC(y))
	if(!is.finite(y)) {
            if(non.finite.report)
	    cat("  y = (y-a)*exp(...) not finite:\n",
		"  xinbta =", format01prec(xinbta, digits=18),"\n")
 	    return(xinbta) ##browser()
	}

	if (y * yprev <= 0) prev <- max(fpu, abs(adj))
	g <- 1
	in.it <- 0
	if(verbose>=2) cat("\n	Inner loop: ")
	for(in.it in 1:1000) {
	    adj <- g * y
	    if (abs(adj) < prev) { #-- current adjustment < last one which gave sign change
		if(verbose>=2)cat("<")
		tx <- xinbta - adj
		if (0 <= tx && tx <= 1) {
		    if (prev <= acu || abs(y) <= acu) { ## goto L_converged
			finish <- TRUE; break
		    }
		    if (tx != 0 && tx != 1) break
		}
	    } else if(verbose>=2) cat(".")
	    if(!finish) g <- g / 3
	} ##-- end inner loop
	if(verbose>=2) cat(if(in.it>10) "\n" else " ")
	if(verbose) cat(" ", in.it,"it --> adj=g*y=", formatC(adj),"\n")
	if (abs(tx - xinbta) <= tol.outer * xinbta) break # goto L_converged;
	if(verbose >= 2)
            cat("--Outer Loop: |tx - xinbta| > eps; tx=",formatC(tx, digits=18),
                " tx-xinbta =", formatC(tx-xinbta),"\n")
	xinbta <- tx
	yprev <- y
	o.it <- o.it + 1
    }, ##- end outer loop {method "AS109"}

    "Newton-log" = #--------------------------------- new --------
        ## but really, I want  Newton on   log-*log* scale --
    while(!finish) { ##-- Outer Loop
	y <- pbeta(xinbta, pp, qq, log.p = TRUE)
	if(verbose) cat("y=pbeta(xi,*, log.p=TRUE)=",
                        formatC(y, digits=15, width=18, flag='-'))
	if(!is.finite(y)) {
	    cat("  y=pbeta(): not finite: 'DOMAIN_ERROR'\n",
		"  xinbta =",formatC(xinbta,digits=18),"\n")
	    return(xinbta) ##browser()
	}
        ## y := g(xinbta) / g'(xinbta)  where g(x) =  log(alpha) - log pbeta(x, ..)
        ##                              ==>  g'(x) = - dbeta(x,.) / pbeta(x,.)
	y <- (y - la) * exp(y + (logbeta + r * log(xinbta) + t * log1p(-xinbta)))
	if(verbose) cat(" y.n=(y-a)*e^..=", formatC(y))
	if(!is.finite(y)) {
	    cat("  y = (y-a)*exp(...) not finite:\n",
		"  xinbta =",formatC(xinbta,digits=18),"\n")
 	    return(xinbta) ##browser()
	}

## 	if (y * yprev <= 0) prev <- max(fpu, abs(adj))
## 	g <- 1
## 	in.it <- 0
## 	if(verbose>=2) cat("\n	Inner loop: ")
## 	for(in.it in 1:1000) {
## 	    adj <- g * y
## 	    if (abs(adj) < prev) { #-- current adjustment < last one which gave sign change
## 		if(verbose>=2)cat("<")
## 		tx <- xinbta - adj
## 		if (0 <= tx && tx <= 1) {
## 		    if (prev <= acu || abs(y) <= acu) { ## goto L_converged
## 			finish <- TRUE; break
## 		    }
## 		    if (tx != 0 && tx != 1) break
## 		}
## 	    } else if(verbose>=2) cat(".")
## 	    if(!finish) g <- g / 3
## 	} ##-- end inner loop
## 	if(verbose>=2) cat(if(in.it>10)"\n" else " ")
## 	if(verbose) cat(" ",in.it,"it --> adj=g*y=",formatC(adj),"\n")
 	if(verbose) cat("\n")
        tx <- xinbta - y
	if (abs(tx - xinbta) < tol.outer * xinbta) break # goto L_converged;
	if(verbose>=2) cat("--Outer Loop: |tx - xinbta| > eps; tx=",formatC(tx, digits=18),
                           " tx-xinbta =", formatC(tx-xinbta),"\n")
	xinbta <- tx
## 	yprev <- y
	o.it <- o.it + 1
    }, ##- end outer loop, method = "Newton-log"
    ## otherwise:
    stop("unknown method ", dQuote(method)))

    ## L_converged:
    if(verbose) cat(" ", o.it,"outer iterations\n")
    if (swap.tail) 1 - xinbta else xinbta
}