Search code examples
pythonpython-3.xscipyscientific-computing

Use a higher number than the maximum float for the special function lambertw


I am using the special function lambertw (k=-1) in Python 3 and I need to use it with numbers higher/lower than the maximum/minimum float number (1.7976931348623157e+308).

What can I do?

Also I tried with "decimal", but it did not work, i. e.,

from decimal import Decimal
from scipy.special import lambertw

lambertw(Decimal('3.1E+600'))

obtained this,

Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/share/apps/sistema/Python-3.5.1/lib/python3.5/site-packages/scipy/special/lambertw.py", line 107, in lambertw
return _lambertw(z, k, tol)
TypeError: ufunc '_lambertw' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''    

Solution

  • The mpmath library in SymPy includes an implementation of lambertw. mpmath implements arbitrary precision floating point arithmetic.

    Here's an example. First, import mpmath from sympy, and set the digits of precision to 100 (chosen arbitrarily--change to fit your needs):

    In [96]: from sympy import mpmath
    
    In [97]: mpmath.mp.dps = 100
    

    Verify that the mpmath function gives the same results as scipy.special.lambertw:

    In [98]: from scipy.special import lambertw
    
    In [99]: lambertw(123.45)
    Out[99]: (3.5491328966138256+0j)
    
    In [100]: mpmath.lambertw(123.45)
    Out[100]: mpf('3.549132896613825444243187580460572741065183903716765715536934583554830913412258511917029758623080475405')
    

    Compute lambertw(3.1e600). The argument is entered as a string, because we can't represent 3.1e600 as a regular floating point value. mpmath will convert the string to a high-precision floating point value using the precision that we set earlier.

    In [101]: mpmath.lambertw('3.1e600')
    Out[101]: mpf('1375.455917376503282959382815269413629072666427317318260231463057587794635136887591876065911283365916388')
    

    We can also create a variable x to hold the input value, and then call mpmath.lambertw(x):

    In [102]: x = mpmath.mpf('3.1e600')
    
    In [103]: x
    Out[103]: mpf('3.099999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999e+600')
    
    In [104]: mpmath.lambertw(x)
    Out[104]: mpf('1375.455917376503282959382815269413629072666427317318260231463057587794635136887591876065911283365916388')
    

    The result can be represented as a regular floating point value, so we pass it to the builtin function float() to do the conversion:

    In [105]: float(mpmath.lambertw(x))
    Out[105]: 1375.455917376503