(The analysis that follows is a light introduction. For more austere results, consult Paper 5).

Lambert's W function is the function that satisfies the functional
equation W(z)*e^{W(z)}=z (1)

Various properties of this function have been presented in the articles: Infinite Exponentials, A Maclaurin Expansion for Lambert's W function and Numerically Approximating the Principal Branch of the Complex Lambert's W function.

This function will henceforth be denoted as W, unless explicitly named. W is multi-valued so it is usually denoted as W(k,z), with k in Z, denoting the appropriate branch and W(0,z)=W(z) denoting the principal branch.

**Preliminaries**

Let us return to equation (1) for a moment. As W satisfies it, it
is obvious that if we have an equation that can be brought into the form
x*A^{x}=y (2), W can solve it as follows:

x*A^{x}=y, =>

x*e^{log(A)*x}=y, =>

log(A)*x*e^{log(A)*x}=log(A)*y, =>

log(A)*x=W(log(A)*y), =>

x=W(log(A)*y)/log(A).

Now let us verify that x as given above, satisfies equation (2):

W(log(A)*y)/log(A)*A^{W(log(A)*y)/log(A)}=

W(log(A)*y)/log(A)*e^{W(log(A)*y)}=

log(A)*y/log(A), (using the definition (1) with z=log(A)*y)

=y.

A quite large family of exponential equations that are otherwise
unsolvable algebraically can be brought into form (2) or a similar form
and as such can be solved by the W function. For example,
the equation c^{x}=x,
can be brought to a form closely resembling (2), and thus we can solve
it analytically by W:

x=c^{x}, =>

x*c^{-x}=1, =>

-x*e^{-x*log(c)}=-1, =>

-x*log(c)*e^{-x*log(c)}=-log(c), =>

-x*log(c)=W(-log(c)), =>

x=-W(-log(c))/log(c).

A similar treatment^{[1]} solves equations of the form: A*c^{x}=B*x,
A,B=/=1, A=/=0. (3)

The principal motivation behind such manipulations is to solve
exponential equations that are otherwise impossible to solve
algebraically. Knowing the behavior of W one can determine immediately
for example whether equation 2^{x}=x has any real solutions
without resorting to numerical methods. In the latter case, solving the equation using W, one gets W(-ln(2))/ln(2).

As W is real valued for x in [-1/e,+∞) and as -ln(2)<-1/e, which
is the branch point of the W, it follows that there are no real values
that satisfy this equation (this could have been determined
by analyzing the function f(x)=2^{x}-x).

After working with several exponential equations, one starts noticing that there exist equations which cannot be brought into form (3). Examples would be the equations:

x*e^{ex}=y.

x*e^{x*ex}=y.

x*2^{x*7x}=y.

Trying to bring any of the above into form (3), would be an exercise in futile circularity. A little tinkering reveals the most general form of equation which is "solvable" by W for z:

a*z^{m}*e^{b*zn}=y.

The question that raises its head now, is, what's keeping us from defining a "function", similar to W, to solve these equations? Well, one has to be very careful with what one means by "function" here. How do we know that if such objects exist, they are well defined?

__Notation__

Let us start with some consistent notation. Assume we are in possession of the following:

S={f_{i}(x):R->R}, i in {1,2,3,...n}, n in N, finite
natural.

We define a function F:R->R, recursively:

F(n,x)={e^{f1(x)}, if n=1, (e^{fn(x)})^{F(n-1,x)}, if n>1}.

The function that will ultimately interest us (these definitions can later be extended from C to C) will be the "inverse" of:

x*F(n,x), n in N, x in R. (4)

The above function in its nonrecursive form, given n in N, it looks like this:

x*(e^{fn(x)})^{(efn-1(x))(efn-2(x))...(ef1(x))}.

The reason the author chose the indexes going from top to bottom, is to conform to the rules of iterated exponentiation and because as such it can be easily defined in Maple.

Now, given S and momentarily disregarding issues of existence and whether such an object exists altogether, we introduce a new "object", called (the) (possibly multi valued) "inverse" of x*F(n,x), by:

HW(f_{1},f_{2},...f_{n};y), (5)

In other words, we hope that this "object" should satisfy:

HW(f_{1},f_{2},...f_{n};y)*F(n,HW(f_{1},f_{2},...f_{n};y))=y. (6)

__Examples__

If n=1, then for f_{1}(x)=x, x*F(1,x)=x*e^{x}, which is the left side of (2), the inverse of
which is the real W.

If n=2, then for f_{1}(x)=log(2)*x and f_{2}(x)=1, x*F(2,x)=x*e^{2x}, a function whose inverse cannot be found via W algebraically.

A more exotic example: If n=3, then for f_{1}(x)=log(3)*x, f_{2}(x)=log(5)*x
and f_{3}(x)=log(7)*x^{2}, x*F(3,x)=x*7^{x2*5x*3x}, a
function whose inverse might be even more exotic.

Using the above notation, the "inverses" of the above examples, would be denoted as:

HW(x;y)=W(y)

HW(log(2)*x,1;y) and

HW(log(3)*x,log(5)*x,log(7)*x^{2};y)

What should the f_{i} look like so that we are guaranteed the
existence of HW? For the real case the answer us easy.

__Lemma #1:__

If f_{i}, i in {1,2,3,...n}, n in N, take values from D <= R
and are such that x*F(n,x) is defined in all of D and is 1-1 there,then HW exists there.

__Proof:__

Standard Calculus. Care has to be taken to ensure that x*F(n,x) is 1-1
in D. For example, if n=2, f_{1}=x and f_{2}=1/x, D
cannot contain any full neighborhood of x~0.6590460684, as this x
is one of the positive roots of the equation e^{x}*x-e^{x}+x=0
and that's where the (local) minimum of x*F(n,x) occurs.

__Numerical Examples__

Let n=1, 2, 3, ..., N and look at x*F(n,x), with f_{i}:R^{+}->R,
f_{i}(x)=x, for all i in {1,2,...,n}. We are then effectively
looking at the inverse of x*{^{n}(e^{x})}, using the tetration operator. The main idea here, is that if we have an analytic f(x) and the
equation f(x)=y, perhaps if we solve the equation T_{m}(x)=y
(for a certain finite m), where T_{m}(x) is the Taylor
polynomial of order m, we may be able to obtain approximations to a
root a of f(x)=y. Once we have such an approximation, we can perhaps
write a=f^{-1}(y), under suitable circumstances.

Our strategy is therefore: Given a certain analytic f and an equation f(x)=y:

- Solve the equation T
_{m}(x)=y. - Test all the solutions of the above equation and see if any of them are close to a root.

Let's try some Maple code and leave the burden of why it works for later:

> N:=6;

> for i from 1 to N do

> f[i]:=z;

> od:

> F:=proc(n)

> if n=1 then exp(f[1]);

> else exp(f[n])^F(n-1);

> fi:

> end:

And the "inverse":

> HW:=proc(y,m,n) #y:value of inverse at.n:Height of tower. m:accuracy.

> local s,p,sol,i,aprx,dy,dist,solution;

> dist:=infinity;#assume infinite distance

> s:=series(z*F(n),z,m);#find Taylor polynomial

> p:=convert(s,polynom);

> sol:={fsolve(p=y,z,complex)};#solve Taylor polynomial

> for i from 1 to nops(sol)do#check all roots one by one

> aprx:=evalf(subs(z=op(i,sol),z*F(n)));#get an estimate

> dy:=evalf(abs(y-aprx));#get epsilon

> if dy<=dist then#if closer

> solution:=op(i,sol);#this root is a better one

> dist:=dy;#update distance

> fi;

> od;

> if Im(solution)<>0 then#imaginary solution

> if Im(y)=0 then #pick Arg(z)>0

> solution:=Re(solution)+abs(Im(solution))*I;

> else #Im(y)<>0#original argument complex, so apply
counterclockwise continuity

> if Im(y)>0 then

> solution:=Re(solution)+abs(Im(solution))*I;

> else

> solution:=Re(solution)-abs(Im(solution))*I;

> fi;

> fi;

> fi;

> solution;

> end:

> plot({seq('HW'(n,y,35),n=1..6)},y=0..3);

From top to bottom, they are the functions:

HW(x;y),

HW(x,x;y),

...

HW(x,x,x,x,x,x;y).

The top function (HW(x;y)) is an old friend: The real W. The
functions appear to converge. We don't know this yet, but note that
informally, when ^{n}(e^{x}) converges, it converges
to e^{-W(-x)}^{[2]}. Accordingly, we are
looking at the inverse of x*e^{-W(-x)}. And here we get a
surprise with Maple:

> solve(x*exp(-W(-x))=y,x);

y/e^{y}, which makes perfect sense, since this is log(z)
whenever the infinite tower z^{zz...}
converges^{[3]}.

Informally thus, HW(x,x,x,...;y)=y/e^{y}.

Note that according to the same reference, we have convergence for positive real y, whenever y<=1. And here are the graphs along with the limit:

__The Complex Case__

Before proceeding we need to explain the code a bit:

First, the choice of "solution" is ambiguous in the do loop, whenever both "solution" and "conjugate(solution)" are in the solution set, since both minimize epsilon. For this, we pick among the two using counter-clockwise continuity, by defining our branch cut to be the negative real axis. If y is in R and "solution" is in C, we pick the one with Arg(solution)>0. If y is in C, we pick between "solution" and "conjugate(solution)" the one whose Arg has the same sign as Arg(y). This guarantees that if we approach the negative real axis from above, the solution's Arg will be positive, while if we approach the negative Real axis from below, the solution's Arg will be negative.

Second, how do we know that increasing m will increase the accuracy of the found root?

__Proof:__

Fix n and let g(z)=z*F(n,z), so we have the equation g(z)=y. Let f(z)=g(z)-y, g_{k}(z)=T_{k}(z), where T_{k}(z)
is the Taylor polynomial of degree k of g(z) and let f_{k}(z)=g_{k}(z)-y. f_{k}(z)->f(z) uniformly and all functions are entire. In particular g(z) is entire and non-constant, so according to Picard's
Little Theorem,=> (Ay)(with at most one exception)(Ea):g(a)=y,=> (Ea)(except in at most one case): f(a)=0.

Once we have a root a, Hurwitz's
Root Theorem provides for a neighborhood |x-a|<ä, and a
number m=N, such that within the neighborhood |x-a|<ä, => (Ak>m=N):f_{k}(z)=0 has a root in that neighborhood. In particular, (Ea_{T}): T_{k}(a_{T})=y, and a_{T}
satisfies: |a_{T}-a|<ä.

Using the continuity of g(z),=> |y-g(a_{T})|<epsilon, and the code above picks the a_{T} that minimizes epsilon. Note that if a_{Tk} is in {z: f_{k}(z)=0}, the uniform convergence of f_{k}(z) guarantees that lim_{k->+∞}inf{|a-a_{Tk}|}=0, so increasing the degree of the Taylor polynomials gives us better and
better approximations of the root a.

Let's see if we can solve some complicated exponential equations
this way. In each case, we have to run the do loop and then substitute
for the functions f_{i} that interest us. Then we run the
iterated exponential code.

If our equation is: x*e^{ex}=3, we are interested
in HW(x,1;3):

> for i from 1 to N do

> pf[i]:=z;

> od:

> f[2]:=1;

> F:=proc(n)

> if n=1 then

> exp(f[1]);

> else

> exp(f[n])^F(n-1);

> fi:

> end:

Our y is 3:

> a:=HW(3,10,2);

a := .5397051810

> evalf(subs(z=a,z*F(2)));

3.000469496

Suppose our equation had been instead: e^{ex}=x. This equation can be immediately brought into the form: x*e^{-ex}=1. In this case we run the code similarly as we are interested in
HW(x,-1;1):

> f[2]:=-1;

Our y now is 1:

> a:=HW(1,10,2);

a := .3144595775 + 1.338932374*I

> evalf(subs(z=a,z*F(2)));

1.003995243 + .007277260574*I

Let's increase the accuracy a bit:

> a:=HW(1,30,2);

a := .3181315052 + 1.337235701*I

> evalf(subs(z=a,z*F(2)));

.9999999991 + .5033177302*10^{-9}*I

A more exotic example: 2^{73x}=x. This equation can be immediately brought into the form: x*2^{-73x}=1, so we are looking for HW(x*log(3),log(7),-log(2);1):

> f[1]:=log(3);

> f[2]:=log(7);

> f[3]:=-log(2);

Our y is again 1:

> a:=HW(1,10,3);

a := -.1907010314 + .7599046786*I

> evalf(subs(z=a,1/F(3)));

-.5821011572 + 2.103917908*I

Let's increase the accuracy a bit:

> a:=HW(1,40,3);

a := .3517427618+.5499294499*I

> evalf(subs(z=a,1/F(3)));

.3453506993+.5666889887*I

Convergence is a little slow in the above case.

A Lambert W-like example: x*e^{x*ex}=24, so we are
looking for HW(x,x;24):

> f[1]:=z;

f[2]:=z;

Our y is 24:

> a:=HW(24,35,2);

a := 1.068699460

> evalf(subs(z=a,z*F(2)));

24.00008895

The above Maple code is an attempt to "invert" z*F(n,z) locally. For a global complex inverse, the procedure is pretty much straightforward. One finds the points where df/dz=0, and one cuts the graph, making sure the cut passes through the images of those points. For example, for W:

d(z*e^{z})/dz=0,=>

e^{z}+z*e^{z}=0,=>

e^{z}*(1+z)=0,=>

z=-1.

So, in the case of the W, the inverse is valid provided we have a
cut that starts at (-1)*e^{-1}=-1/e. Indeed that's where the branch cut of the principal branch of W
starts, extending towards the negative Real axis.

In the code above there is no such proviso, however the code behaves correctly and gives the desired branch cut. Let's verify this. First, the graph of the W as it is coded internally by Maple:

> p[0]:=complexplot3d(W(w),w=-1-I..1+I):

> display(p[0]);

Then, our code:

> p[1]:=complexplot3d('HW'(x+y*I,10,1),-1-I..1+I):

> display(p[1]);

Does our Maple code provide a branch cut starting at -1/e? Let's see:

First, let's check the value of our code @ -1/e:

> HW(-exp(-1),20,1);

-1.000020402

> evalf(W(-exp(-1)));

-1.

Let's move a bit towards 0:

> HW(-exp(-1)+0.01,20,1);

-.7832291994

> evalf(W(-exp(-1)+0.01));

-.7832291993

Now let's move å above and below the Real axis:

> HW(-exp(-1)+0.01+0.0001*I,20,1);

-.7832263387 + .001009590183*I

> evalf(W(-exp(-1)+0.01+0.0001*I));

-.7832263386 + .001009590182*I

> HW(-exp(-1)+0.01-0.0001*I,20,1);

-.7832263387 - .001009590183*I

> evalf(W(-exp(-1)+0.01-0.0001*I));

-.7832263386 - .001009590182*I

It sure looks as if our "function" is continuous at -1/e+å for z approaching above and below. No branch cut there.

Now let's move a little left of -1/e and approach from above and below:

> HW(-exp(-1)-0.01+0.0001*I,20,1);

-.9832468650 + .2314374776*I

> evalf(W(-exp(-1)-0.01+0.0001*I));

-.9809718562 + .2310842101*I

> HW(-exp(-1)-0.01-0.0001*I,20,1);

-.9832468650 - .2314374776*I

> evalf(W(-exp(-1)-0.01-0.0001*I));

-.9809718562 - .2310842101*I

Success! We have a definite positive and negative imaginary quantity there which is not close to 0, which seems to validate that we are approaching a cut!

With these thoughts in mind, let's see the principal branches of some of these "inverses" extended to the complex plane:

> HW(z,z;w):

> p[2]:=complexplot3d('HW'(x+I*y,10,2),-1-I..1+I):

> display(p[2]);

> HW(z,z,z;w):

> p[3]:=complexplot3d('HW'(x+I*y,10,3),-1-I..1+I):

> display(p[3]);

All three together:

> display({p[1],p[2],p[3]});

Let's look at them from a different angle:

Finally, all three, along with the limit function y/e^{y}:

> p[-1]:=complexplot3d(z*exp(-z),z=-1-I..1+I):

> display({p[-1],p[1],p[2],p[3]});

- It is instructive for the reader to try to solve the equation e
^{x}+a*x+b=0, for example, using only definition (1) (answer: x=-W(e^{(-b/a)}/a)-b/a). - See A Series Expansion for
^{n}(e^{x}). - See A Slightly Deeper Analysis of Infinite Exponentials.