ENH: Add scipy.special.logiv and logive

scipy/special/__init__.py

@@ -76,6 +76,8 @@
    kve      -- Exponentially scaled modified Bessel function of the second kind.
    iv       -- Modified Bessel function of the first kind of real order.
    ive      -- Exponentially scaled modified Bessel function of the first kind.
+   logiv    -- Log Modified Bessel function of the first kind of real order.
+   logive   -- Log Exponentially scaled modified Bessel function of the first kind.
    hankel1  -- Hankel function of the first kind.
    hankel1e -- Exponentially scaled Hankel function of the first kind.
    hankel2  -- Hankel function of the second kind.
@@ -655,7 +657,10 @@
     spherical_kn
 )
 
-__all__ = _ufuncs.__all__ + _basic.__all__ + orthogonal.__all__ + [
+from . import log_iv
+from .log_iv import *
+
+__all__ = _ufuncs.__all__ + _basic.__all__ + orthogonal.__all__ + log_iv.__all__ + [
     'SpecialFunctionWarning',
     'SpecialFunctionError',
     'orthogonal',  # Not public, but kept in __all__ for back-compat

scipy/special/_logiv.pyx

@@ -0,0 +1,51 @@
+from libc.math cimport sqrt, log, log1p, M_PI
+
+def log_iv_asym(v, z):
+    r"""Asymptotic expansion for I_{v}(z)
+    for real $z > 0$ and $v\to +\infty$.
+    Based off DLMF 10.41
+
+    References
+    ----------
+    .. [dlmf] NIST Digital Library of Mathematical Functions
+           https://dlmf.nist.gov/10.41#ii
+    .. [a144617] Sequence A144617, The On-Line Encyclopedia of Integer Sequences,
+           https://oeis.org/A144617
+    """
+    cdef:
+        # DLMF 10.41.3 calculates I_v(v*z)
+        double x = z / v
+        # DLMF 10.41.8
+        double inv_p = sqrt(1.0 + x*x)
+        # DLMF 10.41.7
+        double eta = inv_p + log(x) - log1p(inv_p)
+
+    # e^(vn) / ((2piv)^(1/2)*p^(-1/2))
+    log_res = v*eta - log(2.0 * M_PI * v * inv_p)/2
+
+    # large-v asymptotic correction, DLMF 10.41.10
+    p = 1.0/ inv_p
+    p2 = p * p
+    p4 = p2 * p2
+    p6 = p4 * p2
+    p8 = p4 * p4
+    p10 = p6 * p4
+
+    # p over v
+    pv  = p / v
+    pv2 = pv * pv
+    pv3 = pv2 * pv
+    pv4 = pv2 * pv2
+    pv5 = pv3 * pv2
+
+    # terms from OEIS A144617 and A144618
+    u1 = (3.0 - 5.0*p2) / 24.0
+    u2 = (81.0 - 462.0*p2 + 385.0*p4) / 1152.0
+    u3 = (30375 - 369603*p2 + 765765*p4 - 425425*p6) / 414720.0
+    u4 = (4465125 - 94121676*p2 + 349922430*p4 - 446185740*p6 +
+          185910725*p8) / 39813120.0
+    u5 = (1519035525.0 - 49286948607.0*p2 + 284499769554.0*p4 - 614135872350.0*p6 +
+          566098157625.0*p8 - 188699385875.0*p10) / 6688604160.0
+    u_corr_i = 1.0 + u1 * pv + u2 * pv2 + u3 * pv3 + u4 * pv4 + u5 * pv5
+
+    return log_res + log(u_corr_i)

scipy/special/log_iv.py

@@ -0,0 +1,43 @@
+"""
+Log of Bessel functions
+
+References
+----------
+.. [dlmf] NIST Digital Library of Mathematical Functions
+       https://dlmf.nist.gov/10.41#ii
+.. [a144617] Sequence A144617, The On-Line Encyclopedia of Integer Sequences,
+       https://oeis.org/A144617
+"""
+#
+# Author:  Seth Troisi 2020
+
+from numpy import log
+from scipy.special import iv, ive
+from ._logiv import log_iv_asym
+
+__all__ = ['logiv', 'logive']
+
+def logiv(v, z):
+    """
+    Log of modified Bessel function of the first kind of real order.
+
+    See `iv` for Parameters and return
+    """
+
+    if abs(v) <= 50:
+        return log(iv(v, z))
+
+    return log_iv_asym(v, z)
+
+def logive(v, z):
+    """
+    Log of Exponentially scaled modified Bessel function of the first kind of real order.
+
+    See `ive` for Parameters and return
+    """
+
+    if abs(v) <= 50:
+        return log(ive(v, z))
+
+    # scaled by exp(-abs(z.real))
+    return log_iv_asym(v, z) - abs(z)

scipy/special/setup.py

@@ -139,6 +139,10 @@ def configuration(parent_package='',top_path=None):
     config.add_extension('_comb',
                          sources=['_comb.c'])
 
+    # logiv
+    config.add_extension('_logiv',
+                         sources=['_logiv.c'])
+
     # testing for _round.h
     config.add_extension('_test_round',
                          sources=['_test_round.c'],

scipy/special/tests/test_logiv.py

@@ -0,0 +1,32 @@
+import numpy as np
+from numpy.testing import assert_equal, assert_allclose
+from scipy.special import iv
+from scipy.special._logiv import log_iv_asym
+
+
+def assert_iv(v, z, expected_iv):
+    assert_allclose(log_iv_asym(v, z), expected_iv)
+
+
+class TestLogIvAsym(object):
+    def test_int(self):
+        assert_iv(55, 27, -22.011649199378154565)
+
+    def test_float(self):
+        assert_iv(55.0, 27, -22.011649199378154565)
+        assert_iv(55, 27.0, -22.011649199378154565)
+        assert_iv(55.0, 27.0, -22.011649199378154565)
+
+        assert_iv(np.sqrt(90), np.sqrt(80), 2.0740885600035209681)
+
+    def test_nan(self):
+        assert_iv(5, np.inf, np.nan)
+        assert_iv(5, np.nan, np.nan)
+        assert_iv(np.nan, 5, np.nan)
+
+    def test_inf(self):
+        assert_iv(np.inf, 5, -np.inf)
+
+    def test_special(self):
+        for v, z in [(9, 8), (10, 4), (15, 8), (12.3, 7.2), (100, 4.2)]:
+            assert_iv(v, z, np.log(iv(v, z)))