DOUBLE PRECISION FUNCTION LAMBERTW(X) - Lambert W-function

In this routine the Lambert W-function ${\rm LAMBERTW(X)}=W(x)$, i.e. the solution of

$$W(x)\exp(W(x))=x$$ (1)

is calculated.

The following simple, but accurate approximation is used:

$$W(x)\approx \left\{ \begin{array}{l@{\quad:\quad}l} 0.665\cd... ...{1}{\ln(x)})\cdot\ln(\ln(x)) & x > 500\,. \end{array} \right.$$ (2)