swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: b5e4bf28fcba9f5eaffbeecfb0bc307452d074ee authored by Hans W. Borchers on 01 November 2014, 00:00:00 UTC
version 1.7.7
version 1.7.7
Tip revision: b5e4bf2
expint.R
##
## e x p i n t . R Exponential Integral
##
expint <- function(x) {
stopifnot(is.numeric(x) || is.complex(x))
eps <- .Machine$double.eps
x <- c(x)
n <- length(x)
y <- numeric(n)
p <- c(-3.602693626336023e-09, -4.819538452140960e-07, -2.569498322115933e-05,
-6.973790859534190e-04, -1.019573529845792e-02, -7.811863559248197e-02,
-3.012432892762715e-01, -7.773807325735529e-01, 8.267661952366478e+00)
polyv <- polyval(p, Re(x))
# series expansion
k <- which(abs(Im(x)) <= polyv)
if (length(k) != 0) {
# initialization
egamma <- 0.57721566490153286061
xk <- x[k]
yk <- -egamma - log(xk +0i)
j <- 1
pterm <- xk
term <- xk
while (any(abs(term) > eps)) {
yk <- yk + term
j <- j + 1
pterm <- -xk * pterm / j
term <- pterm / j
}
y[k] <- yk
}
# continued fraction
k <- which( abs(Im(x)) > polyv )
if (length(k) != 0) {
m <- 1 # we're calculating E1(x)
# initialization
xk <- x[k]
nk <- length(xk)
am2 <- numeric(nk)
bm2 <- rep(1, nk)
am1 <- rep(1, nk)
bm1 <- xk;
f <- am1 / bm1
oldf <- rep(Inf, nk)
j <- 2
while (any(abs(f - oldf) > (100 * eps) * abs(f))) {
alpha <- m - 1 + (j/2)
# calculate the recursion coefficients
a <- am1 + alpha * am2
b <- bm1 + alpha * bm2
# save new normalized variables for next pass
am2 <- am1 / b
bm2 <- bm1 / b
am1 <- a / b
bm1 <- 1
f <- am1
j <- j + 1
# calculate the coefficients for j odd
alpha <- (j-1)/2
beta <- xk
a <- beta * am1 + alpha * am2
b <- beta * bm1 + alpha * bm2
am2 <- am1 / b
bm2 <- bm1 / b
am1 <- a / b
bm1 <- 1
oldf <- f
f <- am1
j <- j+1
}
y[k] <- exp(-xk) * f - 1i * pi * ((Re(xk) < 0) & (Im(xk) == 0))
}
if (all(Im(y) == 0)) y <- Re(y)
return(y)
}
expint_E1 <- expint # E1()
expint_Ei <- function(x) { # Ei()
stopifnot(is.numeric(x) || is.complex(x))
# y <- -expint(-x) + sign(Im(x)) * pi * 1i
y <- ifelse(sign(Im(x)) <= 0, -expint(-x) - pi*1i, -expint(-x) + pi*1i)
if (all(Im(y) == 0)) y <- Re(y)
return(y)
}
li <- function(x) {
stopifnot(is.numeric(x) || is.complex(x))
y <- expint_Ei(log(x + 0i))
if (all(Im(y) == 0)) y <- Re(y)
return(y)
}