Selection of weps in the Lambert W function

[Yosef]

Hi,

I'm in the process of implementing the Lambert W function for SciPy,
using the mpmath implementation as a reference. There's one thing I
don't understand about it, though, which is the selection of the
iteration precision:

weps = ctx.ldexp(ctx.eps, 15)

I understand that it is meant to be 15 binary orders higher than the
machine epsilon (or just the arbitrary working epsilon?) but I don't
understand why.

If I'm working with machine precision, say float 64-bit (complex 128-
bit) will the same scheme be required?

Thanks for your help,
Yosef.

[Fredrik Johansson]

Hi Yosef,

> I'm in the process of implementing the Lambert W function for SciPy,
> using the mpmath implementation as a reference.

Please note that I just committed a bug fix to the mpmath svn trunk
which you should incorporate.

> There's one thing I
> don't understand about it, though, which is the selection of the
> iteration precision:
>
> weps = ctx.ldexp(ctx.eps, 15)
>
> I understand that it is meant to be 15 binary orders higher than the
> machine epsilon (or just the arbitrary working epsilon?) but I don't
> understand why.
>
> If I'm working with machine precision, say float 64-bit (complex 128-
> bit) will the same scheme be required?

This epsilon is actually 5 orders smaller than the target epsilon
because the working precision is increased by 20 earlier on. For
machine precision, this approach probably needs to be modified. For
the most part, convergence is fast (cubic convergence) and stable, so
it would probably be sufficient to iterate until the error is less
than, say, 1e-8, and then follow up with a single extra iteration.

Unfortunately, the convergence can be slow and near the singular point
-1/e and possibly (for all branches except the 0th) 0. Close to -1/e
in particular, you may want to use a series expansion. Mpmath should
really do this too -- it's just a bit of work to implement and test
for arbitrary precision, but should be easy to do for your purposes. A
series for the 0th branch is given here:
http://functions.wolfram.com/ElementaryFunctions/ProductLog/06/01/02/.
It converges slowly, but you only want to use it very close to -1/e,
where only a few terms are needed.

Fredrik

[Yosef]

On 6 נובמבר, 10:56, Fredrik Johansson <fredrik.johans...@gmail.com>
wrote:

Thanks for the guidance. A first submission of the SciPy code, if
anyone is interested, can be found here:
http://projects.scipy.org/scipy/ticket/1049