ENH: special: add logiv and logive

.gitignore

@@ -244,6 +244,7 @@ scipy/special/_comb.c
 scipy/special/_ellip_harm_2.c
 scipy/special/_ellip_harm_2.h
 scipy/special/_logit.c
+scipy/special/_logiv.c
 scipy/special/_test_round.c
 scipy/special/_ufuncs.c
 scipy/special/_ufuncs.h

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,6 +657,11 @@
     spherical_kn
 )
 
+from ._log_iv import (
+    logiv,
+    logive
+)
+
 __all__ = _ufuncs.__all__ + _basic.__all__ + orthogonal.__all__ + [
     'SpecialFunctionWarning',
     'SpecialFunctionError',
@@ -671,6 +678,8 @@
     'spherical_yn',
     'spherical_in',
     'spherical_kn',
+    'logiv',
+    'logive',
 ]
 
 from scipy._lib._testutils import PytestTester

scipy/special/_log_iv.py

@@ -0,0 +1,70 @@
+"""
+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
+
+import numpy as np
+from scipy._lib._util import _lazywhere
+from scipy.special import iv, ive
+from ._ufuncs import _log_iv_asym
+
+__all__ = ['logiv', 'logive']
+
+
+def _log_value(val, a):
+    # Handle negative, 0, inf, nan
+    # inf => inf, 0 => -inf, negative, nan => nan
+    return _lazywhere(val > 0, (val,),
+        f=lambda v: np.log(v),
+        f2=lambda v: np.where(v == 0, -np.inf, np.nan))
+
+
+def logiv(v, z):
+    """
+    Log of modified Bessel function of the first kind of real order.
+
+    See `iv` for Parameters and return
+    """
+
+    # TODO: if logkv is added can support -v via DLMF 10.27.2
+
+    v = np.asarray(v)
+    z = np.asarray(z)
+
+    value = _lazywhere(
+        (v > 50) & (z > 0),
+        (v, z),
+        f=_log_iv_asym,
+        f2=lambda v, z: _log_value(val=iv(v, z), a=(v, z)))
+
+    # HACK: converts np.array(<scalar>) to <scalar>
+    return value + 0
+
+
+def logive(v, z):
+    """
+    Log of Exponentially scaled modified Bessel function of the first kind of real order.
+
+    See `ive` for Parameters and return
+    """
+
+    v = np.asarray(v)
+    z = np.asarray(z)
+
+    # scaled by exp(-abs(z.real))
+    value = _lazywhere(
+        (v > 50) & (z > 0),
+        (v, z),
+        f=lambda v, z: _log_iv_asym(v, z) - abs(z),
+        f2=lambda v, z: _log_value(val=ive(v, z), a=(v, z)))
+
+    # HACK: converts np.array(<scalar>) to <scalar>
+    return value + 0

scipy/special/_logiv.pxd

@@ -0,0 +1 @@
+cdef double log_iv_asym(double v, double z) nogil

scipy/special/_logiv.pyx

@@ -0,0 +1,59 @@
+import cython
+
+from libc.math cimport sqrt, log, log1p, M_PI
+
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+@cython.cdivision(True)
+cdef inline double log_iv_asym(double v, double z) nogil:
+    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))
+        double log_res = v*eta - log(2.0 * M_PI * v * inv_p) / 2
+
+        double p, p2, p4, p6, p8, p10
+        double pv, pv2, pv3, pv4, pv5
+
+    # 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/add_newdocs.py

@@ -19,6 +19,11 @@ def add_newdoc(name, doc):
     docdict[name] = doc
 
 
+add_newdoc("_log_iv_asym",
+    """
+    Internal function, use `ellip_harm` instead.
+    """)
+
 add_newdoc("_sf_error_test_function",
     """
     Private function; do not use.

scipy/special/functions.json

@@ -15,6 +15,11 @@
             "ellip_harmonic_unsafe": "ddddddd->d"
         }
     },
+    "_log_iv_asym": {
+        "_logiv.pxd": {
+            "log_iv_asym": "dd->d"
+        }
+    },
     "_igam_fac": {
         "cephes.h": {
             "igam_fac": "dd->d"

scipy/special/setup.py

@@ -109,6 +109,11 @@ def configuration(parent_package='',top_path=None):
                          sources=['_ellip_harm_2.c', 'sf_error.c',],
                          **cfg)
 
+    # logiv
+    config.add_extension('_logiv',
+                         sources=['_logiv.c'],
+                         define_macros=define_macros)
+
     # Cython API
     config.add_data_files('cython_special.pxd')
 

scipy/special/tests/test__log_iv.py

@@ -0,0 +1,54 @@
+import numpy as np
+from numpy.testing import assert_equal, assert_allclose
+from scipy.special import iv, logiv, logive
+from scipy.special._log_iv import _log_iv_asym
+
+
+def assert_iv(v, z, expected_iv):
+    # Underlying asymptotic implementation
+    assert_allclose(_log_iv_asym(v, z), expected_iv, rtol=1e-10)
+
+    # Wrapper which chooses between log(iv()) and log_iv_asym()
+    assert_allclose(logiv(v, z), expected_iv, rtol=3e-14)
+    assert_allclose(logive(v, z), expected_iv - abs(z), 3e-14)
+
+class TestLogIv(object):
+    def test_int(self):
+        assert_iv(55, 27, -22.01164919937816)
+
+    def test_float(self):
+        assert_iv(55.0, 27, -22.01164919937816)
+        assert_iv(55, 27.0, -22.01164919937816)
+        assert_iv(55.0, 27.0, -22.01164919937816)
+
+    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_large(self):
+        assert_iv(50, 1, np.log(iv(50, 1)))
+        assert_iv(100, 1, -433.0516183940659)
+        assert_iv(400, 1, -2277.758946766306)
+        assert_iv(800, 1, -5106.468163036196)
+
+        assert_iv(100, 10, np.log(iv(100, 10)))
+        assert_iv(100, 40, np.log(iv(100, 40)))
+        assert_iv(100, 80, np.log(iv(100, 80)))
+
+        sqrt2 = np.sqrt(2)
+        assert_iv(100, 10*sqrt2, np.log(iv(100, 10*sqrt2)))
+        assert_iv(100*sqrt2, 10, np.log(iv(100*sqrt2, 10)))
+
+        assert_iv(1000, 10, -4302.665291340331)
+        assert_iv(1000, 10*sqrt2, -3956.066726969071)
+        assert_iv(1000, 100, -1997.610772811001)
+        assert_iv(1000, 100*sqrt2, -1648.548946015068)
+        assert_iv(1000, 1000, 528.2938870936566)
+        assert_iv(1000, 1000*sqrt2, 1068.924421927351)
+
+        assert_iv(60, 100, 79.19174702982741)
+        assert_iv(60*sqrt2, 1000*sqrt2, 1407.121827193677)

scipy/special/tests/test_mpmath.py

@@ -1386,6 +1386,57 @@ def mp_igam_fac(a, x):
                             [Arg(0, 1e14, inclusive_a=False), Arg(0, 1e14)],
                             rtol=1e-10)
 
+    def test_logiv(self):
+        # small values both real and integer
+        small_values_real = np.arange(-2, 5, 0.74)
+        small_values_int = np.arange(-1, 5)
+
+        # values near the threshold for asym.
+        threshold = np.array([49, 49.87, 50.1, 50.78, 51])
+
+        # large values
+        large_values = np.array([123, 785.2, 967.2, 1300])
+
+        # all values
+        values = np.sort(np.concatenate((
+            small_values_real,
+            small_values_int,
+            threshold,
+            large_values)))
+
+        PREC = 3e-14
+        tests = 0
+        for v, z in itertools.product(values, repeat=2):
+            # These values are not supported by log_iv or iv.
+            if z <= 0 and abs(v) >= 50:
+                continue
+
+            # mpm freezes on these values
+            if v < -20 and 0 < z < 1:
+                continue
+
+            value = mpmath.besseli(v, z)
+            res = mpmath.log(value)
+            test = sc.logiv(v, z)
+
+            if abs(res.imag) >= PREC:
+                assert_(np.isnan(test))
+                continue
+
+            # ignore small imaginary parts
+            assert_(res.imag < PREC)
+            ans = float(res.real)
+
+            if v <= 50 and mpmath.isfinite(res):
+                # value might mismatch from overflow
+                ans = np.log(complex(value))
+
+            assert_allclose(ans, test, rtol=PREC)
+            tests += 1
+
+        # Verify many points were actually tested.
+        assert_(tests >= 530)
+
     def test_j0(self):
         # The Bessel function at large arguments is j0(x) ~ cos(x + phi)/sqrt(x)
         # and at large arguments the phase of the cosine loses precision.