NaN value when using Lambert Function in Python - in Enthought Canopy

Question

I am trying to use the Lambert function in Python to solve a problem; however I am getting a NaN response, when using Canopy. My equation is as follows:

from scipy.special import lambertw

y=8.21016005323e+158

gama = -339.375260893

x = lambertw(y) + gama

print x

When I execute the same code in Matlab I get the value of x = 20.6524 which is the result I am looking for.

I am not sure what this NaN value is being caused by but I suspect it could be something to do with my enormous value for y. Is there any way that I can get Python to deal with this and give me the correct result?

Thanks

scipy.show_config()

   umfpack_info:
  NOT AVAILABLE
lapack_opt_info:
    libraries = ['mkl_lapack95_lp64', 'mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
    library_dirs = ['C:\\Users\\vagrant\\src\\master-env\\libs']
    define_macros = [('SCIPY_MKL_H', None)]
    include_dirs = ['C:\\Users\\vagrant\\src\\master-env\\include']
blas_opt_info:
    libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
    library_dirs = ['C:\\Users\\vagrant\\src\\master-env\\libs']
    define_macros = [('SCIPY_MKL_H', None)]
    include_dirs = ['C:\\Users\\vagrant\\src\\master-env\\include']
openblas_info:
  NOT AVAILABLE
lapack_mkl_info:
    libraries = ['mkl_lapack95_lp64', 'mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
    library_dirs = ['C:\\Users\\vagrant\\src\\master-env\\libs']
    define_macros = [('SCIPY_MKL_H', None)]
    include_dirs = ['C:\\Users\\vagrant\\src\\master-env\\include']
blas_mkl_info:
    libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
    library_dirs = ['C:\\Users\\vagrant\\src\\master-env\\libs']
    define_macros = [('SCIPY_MKL_H', None)]
    include_dirs = ['C:\\Users\\vagrant\\src\\master-env\\include']
mkl_info:
    libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
    library_dirs = ['C:\\Users\\vagrant\\src\\master-env\\libs']
    define_macros = [('SCIPY_MKL_H', None)]
    include_dirs = ['C:\\Users\\vagrant\\src\\master-env\\include']

Answer 1

There is an iteration in the lambertw code. Apparently it is not converging when given a large argument. (And as @unutbu's answer shows, whether or not it converges appears to depend on your configuration.)

Here's an alternative that works for (scalar) positive real arguments up to the maximum floating point value:

import numpy as np
from scipy.optimize import fsolve


def lw(x):
    """Lambert W function, for real x >= 0."""

    def func(w, x):
        return np.log(x) - np.log(w) - w

    if x == 0:
        return 0
    if x > 2.5:
        lnx = np.log(x)
        w0 = lnx - np.log(lnx)
    elif x > 0.25:
        w0 = 0.8 * np.log(x + 1)
    else:
        w0 = x * (1.0 - x)

    return fsolve(func, w0, args=(x,))[0]

For example:

In [79]: lw(8.21016005323e+158)
Out[79]: 360.02763631519991

In [80]: np.finfo(1.0).max
Out[80]: 1.7976931348623157e+308

In [81]: lw(np.finfo(1.0).max)
Out[81]: 703.22703310477016

---

Here's my configuration:

In [87]: scipy.show_config()
atlas_threads_info:
  NOT AVAILABLE
blas_opt_info:
    extra_link_args = ['-Wl,-framework', '-Wl,Accelerate']
    extra_compile_args = ['-msse3', '-I/System/Library/Frameworks/vecLib.framework/Headers']
    define_macros = [('NO_ATLAS_INFO', 3)]
atlas_blas_threads_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
lapack_opt_info:
    extra_link_args = ['-Wl,-framework', '-Wl,Accelerate']
    extra_compile_args = ['-msse3']
    define_macros = [('NO_ATLAS_INFO', 3)]
atlas_info:
  NOT AVAILABLE
lapack_mkl_info:
  NOT AVAILABLE
blas_mkl_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
mkl_info:
  NOT AVAILABLE