[R] Lambert's W function

[Ben Bolker]

  This implementation, originally written by Nici Schraudolph, allows you
to choose which branch you want.  I've checked the answers for complex
arguments, non-systematically, against Mathematica's ProductLog function.

  Ben Bolker

lambertW = function(z,b=0,maxiter=10,eps=.Machine$double.eps,
  min.imag=1e-9) {
  if (any(round(Re(b)) != b))
    stop("branch number for W must be an integer")
  if (!is.complex(z) && any(z<0)) z=as.complex(z)
  ## series expansion about -1/e
  ##
  ## p = (1 - 2*abs(b)).*sqrt(2*e*z + 2);
  ## w = (11/72)*p;
  ## w = (w - 1/3).*p;
  ## w = (w + 1).*p - 1
  ##
  ## first-order version suffices:
  ##
  w = (1 - 2*abs(b))*sqrt(2*exp(1)*z + 2) - 1
  ## asymptotic expansion at 0 and Inf
  ##
  v = log(z + as.numeric(z==0 & b==0)) + 2*pi*b*1i;
  v = v - log(v + as.numeric(v==0))
  ## choose strategy for initial guess
  ##
  c = abs(z + exp(-1));
  c = (c > 1.45 - 1.1*abs(b));
  c = c | (b*Im(z) > 0) | (!Im(z) & (b == 1))
  w = (1 - c)*w + c*v
  ## Halley iteration
  ##
  for (n in 1:maxiter) {
    p = exp(w)
    t = w*p - z
    f = (w != -1)
    t = f*t/(p*(w + f) - 0.5*(w + 2.0)*t/(w + f))
    w = w - t
    if (abs(Re(t)) < (2.48*eps)*(1.0 + abs(Re(w)))
        && abs(Im(t)) < (2.48*eps)*(1.0 + abs(Im(w))))
      break
  }
  if (n==maxiter) warning(paste("iteration limit (",maxiter,
        ") reached, result of W may be inaccurate",sep=""))
  if (all(Im(w)<min.imag)) w = as.numeric(w)
  return(w)
}

---------

\name{lambertW}
\alias{lambertW}
\title{Lambert W function}
\description{
  Computes the Lambert W function, giving efficient solutions to the equation x*exp(x)==x
}
\usage{
lambertW(z, b = 0, maxiter = 10, eps = .Machine$double.eps, min.imag = 1e-09)
}
\arguments{
  \item{z}{(complex) vector of values for which to compute the function}
  \item{b}{(integer) vector of branches: b=0 specifies the principal
    branch, 0 and -1 are the ones that can take non-complex values}
  \item{maxiter}{maximum numbers of iterations for convergence}
  \item{eps}{convergence tolerance}
  \item{min.imag}{maximum magnitude of imaginary part to chop when
    returning solutions}
}
\details{
Compute the Lambert W function of z.  This function satisfies
W(z)*exp(W(z)) = z, and can thus be used to express solutions
of transcendental equations involving exponentials or logarithms.
The Lambert W function is also available in 
Mathematica (as the ProductLog function), and in Maple.
}
\value{
  Complex or real vector of solutions.
}
\references{Corless, Gonnet, Hare, Jeffrey, and Knuth (1996), "On the Lambert
W Function", Advances in Computational Mathematics 5(4):329-359}
\author{Nici Schraudolph <schraudo at inf.ethz.ch> (original
  version (c) 1998), Ben Bolker (R translation)
  }
\note{
This implementation should return values within 2.5*eps of its
counterpart in Maple V, release 3 or later.  Please report any
discrepancies to the author or translator.
}
\examples{
curve(lambertW(x),from=0,to=10)
pvec = seq(-1,1,length=40)
m = outer(pvec,pvec,function(x,y)Re(lambertW(x+y*1i)))
persp(pvec,pvec,m,
      theta=290,shade=0.5,zlab="lambertW")
num1 = uniroot(function(x) {x*exp(x)-1},lower=0,upper=1,tol=1e-9)
abs(lambertW(1)-num1$root)<1e-9
}
\keyword{math}

-- 
620B Bartram Hall                            bolker at zoo.ufl.edu
Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
Box 118525                                   (ph)  352-392-5697
Gainesville, FL 32611-8525                   (fax) 352-392-3704