A numerical power series evaluation needs a stopping rule. For and integer , the absolute remainder after terms satisfies the polylogarithm series remainder boundIndeed for , and the remaining powers sum as a geometric series. This also covers negative . The code sums terms in vectors of at most 4096 entries, checking the bound after each block; there is no loop over individual terms.
polylog <- function(a, z, tol = 1e-10, max_terms = 1000000) {
scalar <- function(v) {
is.numeric(v) && !is.complex(v) && length(v) == 1L &&
is.finite(v)
}
if (!scalar(a) || a < 0 || a != floor(a))
stop("a must be a nonnegative integer")
if (!scalar(z) || abs(z) >= 1)
stop("z must be real with abs(z) < 1")
if (!scalar(tol) || tol <= 0 || tol >= 1)
stop("tol must lie strictly between 0 and 1")
if (!scalar(max_terms) || max_terms < 1 ||
max_terms != floor(max_terms))
stop("max_terms must be a positive integer")
if (z == 0) return(0)
r <- abs(z)
N <- 0
ans <- 0
repeat {
if (N >= max_terms)
stop("requested truncation tolerance not reached")
B <- min(4096, max_terms - N)
k <- seq.int(N + 1, N + B)
ans <- ans + sum(z^k * exp(-a * log(k)))
N <- N + B
log_bound <- (N + 1) * log(r) - a * log(N + 1) - log1p(-r)
if (log_bound <= log(tol)) return(ans)
}
}This accepts every positive integer order in the question, and also order zero if is the convention adopted. The bound controls absolute truncation error; floating-point rounding is a separate source of error. Close to , convergence can be slow. If the resource limit is reached, the function reports failure instead of silently claiming the requested accuracy. Useful exact checks are and .