# [R] Lambert's W function

Stefan Böhringer commercial at s-boehringer.de
Wed Nov 26 11:25:03 CET 2003

```On Tue, 2003-11-25 at 12:09, Robin Hankin wrote:
> Hello List
>
> does anyone have an R function for the Lambert W function?  I need
> complex arguments.
>
> [the Lamert W function W(z) satisfies
>
> W(z)*exp(W(z)) = z
>
> but I could'nt even figure out how to use uniroot() for complex z]
>
>

There are several 'branches' of W(z) for complex arguments, since the
inverse of f(z) = z * exp(z) does not exist. Therefore you have to
constrain root-finding methods appropriately.

Here is an implementation for real valued arguments, if that is of any
help to you:

WE1 <- function(x) {
if (1) {
o <- optim(c(1),
# optimize squared difference of z and log(t) + t
function(t) { (x - (log(t) + t))**2 },
# the gradient of the above function
#function(t) { 2(log(t) + t - x)(1/t - 1) },
method="L-BFGS-B", lower=exp(-70),
control = list(lmm=40, maxit=10000, factr=1e9));
} else {
o <- optim(c(1),
# optimize squared difference of z and log(t) + t
function(t) { (x - (log(t) + t))**2 },
method="BFGS");
}
o\$par
}

W2 <- function(x) {
eps <- 10^-5;	# precision

if (x == 0) {
wnew <- 0
} else {
# initial guess
wold <- if (x < 2.5) {
(x + 4/3*x^2)/(1 + 7/3*x + 5/6*x^2)
} else {
log(x)
};
if (x < 0) print("error");
lx <- log(x);
wnew <- 2 * wold;
if (x != Inf) {
while (abs( (wnew - wold)/wold ) > eps) {
wold <- wnew
zn <- lx - wold - log(wold);
y <- 2*(1 + wold)*((1 + wold) + 2/3*zn) - zn;
wnew <- wold*(1 + (zn * y)/((1 + wold)*(y - zn)));
}
} else {
wnew <- Inf;
}
}
wnew
}
WE <- function(z) {
if ( z < -1 ) { W2(exp(z)) }
else { WE1(z) };
}
W <- function(z) { WE(log(z)) }

Stefan

```