logcf: Continued Fraction Approximation of Log-Related Power Series

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

Continued Fraction Approximation of Log-Related Power Series

Description

Compute a continued fraction approximation to the series (infinite sum)

\sum_{k=0}^\infty \frac{x^k}{i +k\cdot d} = \frac{1}{i} + \frac{x}{i+d} + \frac{x^2}{i+2*d} + \frac{x^3}{i+3*d} + \ldots

Needed as auxiliary function in log1pmx() and lgamma1p().

Usage

					
logcfR (x, i, d, eps, maxit = 10000L, trace = FALSE)
logcfR.(x, i, d, eps, maxit = 10000L, trace = FALSE)
logcf  (x, i, d, eps, trace = FALSE)

Arguments

x	numeric vector of values typically less than 1. "mpfr" (of potentially high precision, package Rmpfr) work in logcfR*(x,*).
i	positive numeric
d	non-negative numeric
eps	positive number, the convergence tolerance.
maxit	a positive integer, the maximal number of iterations or terms in the truncated series used.
trace	logical (or non-negative integer in the future) indicating if (and how much) diagnostic output should be printed to the console during the computations.

Details

logcfR.():

a pure R version where the iterations happen vectorized in x, only for those components x[i] they have not yet converged. This is particularly beneficial for not-very-short "mpfr" vectors x, and still conceptually equivalent to the logcfR() version.

logcfR():

a pure R version where each x[i] is treated separately, hence “properly” vectorized, but slowly so.

logcf():

only for numeric x, calls into (a clone of) R's own (non-API currently) logcf() C Rmathlib function.

Value

a numeric-alike vector with the same attributes as x. For the logcfR*() versions, an "mpfr" vector if x is one.

Note

Rescaling is done by (namespace hidden) “global” scalefactor which is 2^{256}, represented exactly (in double precision).

Author(s)

Martin Maechler, based on R's ‘nmath/pgamma.c’ implementation.

See Also

lgamma1p, log1pmx, and pbeta, whose prinicipal algorithm has evolved from TOMS 708.

Examples

x <- (-2:1)/2
logcf (x, 2,3, eps=1e-7, trace=TRUE) # shows iterations for each x[]
logcfR(x, 2,3, eps=1e-7, trace=TRUE) # 1 line per x[]
logcfR(x, 2,3, eps=1e-7, trace= 2  ) # shows iterations for each x[]

n <- 2049; x <- seq(-1,1, length.out = n)[-n] ; stopifnot(diff(x) == 1/1024)
plot(x, (lcf <- logcf(x, 2,3, eps=1e-12)), type="l", col=2)
lcR <- logcfR (x, 2,3, eps=1e-12); all.equal(lcf, lcR , tol=0)
lcR.<- logcfR.(x, 2,3, eps=1e-12); all.equal(lcf, lcR., tol=0)
stopifnot(exprs = {
  all.equal(lcf, lcR., tol=1e-14)# seen 0 (!)
  all.equal(lcf, lcR,  tol=1e-14)# seen 0 (!) -- failed for a while
})

l32 <- curve(logcf(x, 3,2, eps=1e-7), -3, 1)
abline(h=0,v=1, lty=3, col="gray50")
plot(y~x, l32, log="y", type = "o", main = "logcf(*, 3,2)  in log-scale")

DPQ documentation built on Nov. 3, 2024, 3 a.m.