BUG: fix special case z = -1/e in lambertw

[ilanschnell]

This PR fixes: scipy/scipy#24770

After inspecting the output of scipy's lambertw function near -1/e+0j, I found that only the branch point -1/e itself is the problem, as it will cause division by zero in Halley's method. When changing even the least significant (of the real or imaginary part of z), the method converges to the correct result.

[fbourgey]

Thanks @ilanschnell. A test in tests/xsf_tests would be nice to have, notably for values around -1/e.

[WarrenWeckesser]

Thanks @ilanschnell. However, returning -1 isn't quite right. The double precision floating point value -np.exp(-1) (which is -0.36787944117144233) is less than the mathematical value of $$-e^{-1}$$, and the difference is enough to make the correct value of the double precision calculation of lambertw(-np.exp(-1)) a complex value with a nonzero imaginary part. The following code uses mpmath to compute the correct double precision values of lambertw(-np.exp(-1), k) for the branches k=-1 and k=0:

In [314]: from mpmath import mp

In [315]: mp.dps = 200

In [316]: bp = -np.exp(-1)

In [317]: complex(mp.lambertw(bp, k=-1))
Out[317]: (-1-8.220079714836618e-09j)

In [318]: complex(mp.lambertw(bp, k=0))
Out[318]: (-1+8.220079714836618e-09j)

This is the same kind of phenomenon observed when computing, for example, np.sin(np.pi). The correct value of the latter is 1.2246467991473532e-16, not 0, because np.pi does not equal $$\pi$$:

In [319]: np.sin(np.pi).item()
Out[319]: 1.2246467991473532e-16

In [320]: float(mp.sin(np.pi))
Out[320]: 1.2246467991473532e-16

[ilanschnell]

Thanks Warren, this is a good point. I've add the imaginary value.
EXPN1 is defined as 0.36787944117144232159553, which by far exceeds the double floating point precision.
This makes we wonder, it I should change the condition from z == -EXPN1 to z == -0.36787944117144233, or change the definition to: EXPN1 = 0.36787944117144233.

I'm working on adding tests/xsf_tests/test_lambertw.cpp now.

[WarrenWeckesser]

EXPN1 is defined as 0.36787944117144232159553, which by far exceeds the double floating point precision.

It is a double precision floating point constant. The compiler won't store more than what fits in a double precision float, so writing it that way is no more precise than 0.36787944117144233. You can see the same thing in Python:

In [327]: 0.36787944117144232159553
Out[327]: 0.36787944117144233

So there shouldn't be any need to change the literal constant or the expression that checks z == -EXPN1. I guess the definition of EXPN1 could be changed to 0.36787944117144233, so future readers of the code see a value that is representative of the actual value that the compiler will use, but that is not necessary in this PR.

If the literal value ended with the l suffix, the compiler would parse it as a long double, but there is no point in doing that, since the precision of long double is platform-dependent (and is the same as double on some platforms), and rest of the variables are just double precision.

[ilanschnell]

While working on the new tests, I realized that when lowering the tolerance lambertw still gives nan in some cases when z is close to -1/e:

In [27]: from scipy.special import lambertw

In [28]: lambertw(-0.36787944, -1)
Out[28]: np.complex128(-1.0000000095529311+0j)

In [29]: lambertw(-0.36787944, -1, 1e-12)
Out[29]: np.complex128(nan+nanj)

Therefore the comment

 * TODO: use a series expansion when extremely close to the branch point
 * at `-1/e` and make sure that the proper branch is chosen there.

remains valid, However, as no nan-values are returned when using scipy's default tolerance (1e-8), I consider this less critical now. I'm happy to work on the series expansion near -1/e in a separate PR once this PR gets merged.