In all papers I've seen on the actual computation of the Lambert W function, there is a Halley's method iteration. I expanded and factored terms to create new variables in such a way as to make the computation a lot easier on the machine, but the following is exactly the same as the expression for $\omega_{n+1}$ found on Wikipedia, say. As far as I've been able to experimentally tell, it precisely and extremely quickly converges for all complex values, except those that lie on the line $(-\infty,-\frac{1}{e})$ and have zero imaginary part. As soon as I enter in such a value I get terribe results - and I understand that the branch cut for the principal branch lies on this line, which maybe explains the difficulties, but I do not understand how to edit the Halley's method to all of a sudden converge for these values. I also note that if $\Re(z)<-\frac{1}{e}$, (but $\Im(z)\neq0$) then it requires many iterations, up to $100$, to converge, instead of the normal $5$ - it converges very slowly for the left half of the complex plane - and converges to the wrong branch! - and not at all for the negative real line. How might we correct this? Corless et al. knew what they were doing when they stated this Halley's method, so they must have had a solution!
Letting $\omega\exp(\omega)=z$, then $\lim_{n\to\infty}\omega_n=\omega$ - unless $\Re(z)\in(-\infty,-\frac{1}{e}),\Im(z)=0$. I use $\omega_0=\ln|\Re(z)|+0i$ as that generally follows the trend of the Lambert W for positive $z$ - and I assumed by taking the absolute value that that would achieve a similar similarity, but perhaps it doesn't, and this explains why it converges to the wrong branch for the negative half-plane. If so, how might I better set an initial guess to guarantee the iteration always iterates within the range of the principal branch? I understand the branches are separated by the lines $-\Re(\omega)\cot \Re(\omega)+\Re(\omega)i$, but I clearly don't understand how to use this fact to keep my values in the right branch.
$$\begin{align*} \\ &\psi_n=\exp(\omega_n),\space\kappa_n=\omega_n+1,\space\varphi_n=z\cdot\kappa_n,\space\mu_n=\psi_n\cdot\kappa_n \\ \\ &\delta_n=2\cdot\frac{\varphi_n-\mu\cdot\omega}{\mu\cdot\kappa + \psi_n+\varphi_n+z} \\ \\ &\psi_{n+1}=\psi_n\cdot\exp(\delta_n) \\ \\ &\kappa_{n+1}=\kappa_n+\delta_n \\ \\ &\varphi_{n+1}=\varphi_n+z\cdot\delta_n \\ \\ &\mu_{n+1}=\psi_{n+1}\cdot\kappa_{n+1} \\ \\ &\omega_{n+1}=\omega_n+\delta_n \end{align*} $$
The only reason I show the method I am using in my code instead of the original one is because, although I have stepped through the expansion and shown the two iterations to be identical, perhaps there is some quirk of my characterisation of it that leads to problems.
Many thanks.