ENH: special: extend sici/shichi to complex arguments

scipy/special/_sici.pxd

@@ -0,0 +1,149 @@
+# Implementation of sin/cos/sinh/cosh integrals for complex arguments
+#
+# Sources
+# [1] Fredrik Johansson and others. mpmath: a Python library for
+#     arbitrary-precision floating-point arithmetic (version 0.19),
+#     December 2013. http://mpmath.org/.
+# [2] NIST, "Digital Library of Mathematical Functions",
+#     http://dlmf.nist.gov/
+import cython
+cimport numpy as np
+
+cimport sf_error
+
+cdef extern from "numpy/npy_math.h":
+    double NPY_PI
+    double NPY_PI_2
+    double NPY_EULER
+
+from _complexstuff cimport (
+    npy_cdouble_from_double_complex, double_complex_from_npy_cdouble,
+    inf, nan, zabs, zlog, zpack)
+
+cdef extern from "specfun_wrappers.h":
+    np.npy_cdouble cexpi_wrap(np.npy_cdouble) nogil
+
+DEF MAXITER = 100
+DEF TOL = 2.220446092504131e-16
+
+
+cdef inline double complex zexpi(double complex z) nogil:
+    cdef np.npy_cdouble r
+    r = cexpi_wrap(npy_cdouble_from_double_complex(z))
+    return double_complex_from_npy_cdouble(r)
+
+
+@cython.cdivision(True)
+cdef inline void power_series(int sgn, double complex z,
+                             double complex *s, double complex *c) nogil:
+    """DLMF 6.6.5 and 6.6.6. If sgn = -1 computes si/ci, and if sgn = 1
+    computes shi/chi.
+
+    """
+    cdef:
+        int n
+        double complex fac, term1, term2
+
+    fac = z
+    s[0] = fac
+    c[0] = 0        
+    for n in range(1, MAXITER):
+        fac *= sgn*z/(2*n)
+        term2 = fac/(2*n)
+        c[0] += term2
+        fac *= z/(2*n + 1)
+        term1 = fac/(2*n + 1)
+        s[0] += term1
+        if zabs(term1) < TOL*zabs(s[0]) and zabs(term2) < TOL*zabs(c[0]):
+            break
+
+
+cdef inline int csici(double complex z,
+                      double complex *si, double complex *ci) nogil:
+    """Compute sin/cos integrals at complex arguments. The algorithm
+    largely follows that of [1].
+
+    """
+    cdef double complex jz, term1, term2
+
+    if z == inf:
+        si[0] = NPY_PI_2
+        ci[0] = 0
+        return 0
+    elif z == -inf:
+        si[0] = -NPY_PI_2
+        ci[0] = 1j*NPY_PI
+        return 0
+    elif zabs(z) < 0.8:
+        # Use the series to avoid cancellation in si
+        power_series(-1, z, si, ci)
+        if z == 0:
+            sf_error.error("sici", sf_error.DOMAIN, NULL)
+            ci[0] = zpack(-inf, nan)
+        else:
+            ci[0] += NPY_EULER + zlog(z)
+        return 0
+
+    # DLMF 6.5.5/6.5.6 plus DLMF 6.4.4/6.4.6/6.4.7
+    jz = 1j*z
+    term1 = zexpi(jz)
+    term2 = zexpi(-jz)
+    si[0] = -0.5j*(term1 - term2)
+    ci[0] = 0.5*(term1 + term2)
+    if z.real == 0:
+        if z.imag > 0:
+            ci[0] += 1j*NPY_PI_2
+        elif z.imag < 0:
+            ci[0] -= 1j*NPY_PI_2
+    elif z.real > 0:
+        si[0] -= NPY_PI_2
+    else:
+        si[0] += NPY_PI_2
+        if z.imag >= 0:
+            ci[0] += 1j*NPY_PI
+        else:
+            ci[0] -= 1j*NPY_PI
+
+    return 0
+
+
+cdef inline int cshichi(double complex z,
+                        double complex *shi, double complex *chi) nogil:
+    """Compute sinh/cosh integrals at complex arguments. The algorithm
+    largely follows that of [1].
+
+    """
+    cdef double complex term1, term2
+
+    if z == inf:
+        shi[0] = inf
+        chi[0] = inf
+        return 0
+    elif z == -inf:
+        shi[0] = -inf
+        chi[0] = inf
+        return 0
+    elif zabs(z) < 0.8:
+        # Use the series to avoid cancellation in shi
+        power_series(1, z, shi, chi)
+        if z == 0:
+            sf_error.error("shichi", sf_error.DOMAIN, NULL)
+            chi[0] = zpack(-inf, nan)
+        else:
+            chi[0] += NPY_EULER + zlog(z)
+        return 0
+
+    term1 = zexpi(z)
+    term2 = zexpi(-z)
+    shi[0] = 0.5*(term1 - term2)
+    chi[0] = 0.5*(term1 + term2)
+    if z.imag > 0:
+        shi[0] -= 0.5j*NPY_PI
+        chi[0] += 0.5j*NPY_PI
+    elif z.imag < 0:
+        shi[0] += 0.5j*NPY_PI
+        chi[0] -= 0.5j*NPY_PI
+    elif z.real < 0:
+        chi[0] += 1j*NPY_PI
+
+    return 0

scipy/special/_ufuncs.pyx

@@ -1894,8 +1894,14 @@ cdef extern from "_ufuncs_defs.h":
     cdef double _func_round "round"(double) nogil
 cdef extern from "_ufuncs_defs.h":
     cdef int _func_shichi "shichi"(double, double *, double *) nogil
+from _sici cimport cshichi as _func_cshichi
+ctypedef int _proto_cshichi_t(double complex, double complex *, double complex *) nogil
+cdef _proto_cshichi_t *_proto_cshichi_t_var = &_func_cshichi
 cdef extern from "_ufuncs_defs.h":
     cdef int _func_sici "sici"(double, double *, double *) nogil
+from _sici cimport csici as _func_csici
+ctypedef int _proto_csici_t(double complex, double complex *, double complex *) nogil
+cdef _proto_csici_t *_proto_csici_t_var = &_func_csici
 cdef extern from "_ufuncs_defs.h":
     cdef double _func_sindg "sindg"(double) nogil
 cdef extern from "_ufuncs_defs.h":
@@ -11465,69 +11471,173 @@ ufunc_round_data[0] = &ufunc_round_ptr[2*0]
 ufunc_round_data[1] = &ufunc_round_ptr[2*1]
 round = np.PyUFunc_FromFuncAndData(ufunc_round_loops, ufunc_round_data, ufunc_round_types, 2, 1, 1, 0, "round", ufunc_round_doc, 0)
 
-cdef np.PyUFuncGenericFunction ufunc_shichi_loops[2]
-cdef void *ufunc_shichi_ptr[4]
-cdef void *ufunc_shichi_data[2]
-cdef char ufunc_shichi_types[6]
+cdef np.PyUFuncGenericFunction ufunc_shichi_loops[4]
+cdef void *ufunc_shichi_ptr[8]
+cdef void *ufunc_shichi_data[4]
+cdef char ufunc_shichi_types[12]
 cdef char *ufunc_shichi_doc = (
-    "shichi(x)\n"
+    "shichi(x, out=None)\n"
     "\n"
-    "Hyperbolic sine and cosine integrals\n"
+    "Hyperbolic sine and cosine integrals.\n"
+    "\n"
+    "The hyperbolic sine integral is\n"
+    "\n"
+    ".. math::\n"
+    "\n"
+    "  \\int_0^x \\frac{\\sinh{t}}{t}dt\n"
+    "\n"
+    "and the hyperbolic cosine integral is\n"
+    "\n"
+    ".. math::\n"
+    "\n"
+    "  \\gamma + \\log(x) + \\int_0^x \\frac{\\cosh{t} - 1}{t} dt\n"
+    "\n"
+    "where :math:`\\gamma` is Euler's constant and :math:`\\log` is the\n"
+    "principle branch of the logarithm.\n"
+    "\n"
+    "Parameters\n"
+    "----------\n"
+    "x : array_like\n"
+    "    Real or complex points at which to compute the hyperbolic sine\n"
+    "    and cosine integrals.\n"
     "\n"
     "Returns\n"
     "-------\n"
-    "shi\n"
-    "    ``integral(sinh(t)/t, t=0..x)``\n"
-    "chi\n"
-    "    ``eul + ln x + integral((cosh(t)-1)/t, t=0..x)``\n"
-    "    where ``eul`` is Euler's constant.")
+    "si : ndarray\n"
+    "    Hyperbolic sine integral at ``x``\n"
+    "ci : ndarray\n"
+    "    Hyperbolic cosine integral at ``x``\n"
+    "\n"
+    "Notes\n"
+    "-----\n"
+    "For real arguments with ``x < 0``, ``chi`` is the real part of the\n"
+    "hyperbolic cosine integral. For such points ``chi(x)`` and ``chi(x\n"
+    "+ 0j)`` differ by a factor of ``1j*pi``.\n"
+    "\n"
+    "For real arguments the function is computed by calling Cephes'\n"
+    "[1]_ *shichi* routine. For complex arguments the algorithm is based\n"
+    "on Mpmath's [2]_ *shi* and *chi* routines.\n"
+    "\n"
+    "References\n"
+    "----------\n"
+    ".. [1] Cephes Mathematical Functions Library,\n"
+    "       http://www.netlib.org/cephes/index.html\n"
+    ".. [2] Fredrik Johansson and others.\n"
+    "       \"mpmath: a Python library for arbitrary-precision floating-point arithmetic\"\n"
+    "       (Version 0.19) http://mpmath.org/")
 ufunc_shichi_loops[0] = <np.PyUFuncGenericFunction>loop_i_d_dd_As_f_ff
 ufunc_shichi_loops[1] = <np.PyUFuncGenericFunction>loop_i_d_dd_As_d_dd
+ufunc_shichi_loops[2] = <np.PyUFuncGenericFunction>loop_i_D_DD_As_F_FF
+ufunc_shichi_loops[3] = <np.PyUFuncGenericFunction>loop_i_D_DD_As_D_DD
 ufunc_shichi_types[0] = <char>NPY_FLOAT
 ufunc_shichi_types[1] = <char>NPY_FLOAT
 ufunc_shichi_types[2] = <char>NPY_FLOAT
 ufunc_shichi_types[3] = <char>NPY_DOUBLE
 ufunc_shichi_types[4] = <char>NPY_DOUBLE
 ufunc_shichi_types[5] = <char>NPY_DOUBLE
+ufunc_shichi_types[6] = <char>NPY_CFLOAT
+ufunc_shichi_types[7] = <char>NPY_CFLOAT
+ufunc_shichi_types[8] = <char>NPY_CFLOAT
+ufunc_shichi_types[9] = <char>NPY_CDOUBLE
+ufunc_shichi_types[10] = <char>NPY_CDOUBLE
+ufunc_shichi_types[11] = <char>NPY_CDOUBLE
 ufunc_shichi_ptr[2*0] = <void*>_func_shichi
 ufunc_shichi_ptr[2*0+1] = <void*>(<char*>"shichi")
 ufunc_shichi_ptr[2*1] = <void*>_func_shichi
 ufunc_shichi_ptr[2*1+1] = <void*>(<char*>"shichi")
+ufunc_shichi_ptr[2*2] = <void*>_func_cshichi
+ufunc_shichi_ptr[2*2+1] = <void*>(<char*>"shichi")
+ufunc_shichi_ptr[2*3] = <void*>_func_cshichi
+ufunc_shichi_ptr[2*3+1] = <void*>(<char*>"shichi")
 ufunc_shichi_data[0] = &ufunc_shichi_ptr[2*0]
 ufunc_shichi_data[1] = &ufunc_shichi_ptr[2*1]
-shichi = np.PyUFunc_FromFuncAndData(ufunc_shichi_loops, ufunc_shichi_data, ufunc_shichi_types, 2, 1, 2, 0, "shichi", ufunc_shichi_doc, 0)
-
-cdef np.PyUFuncGenericFunction ufunc_sici_loops[2]
-cdef void *ufunc_sici_ptr[4]
-cdef void *ufunc_sici_data[2]
-cdef char ufunc_sici_types[6]
+ufunc_shichi_data[2] = &ufunc_shichi_ptr[2*2]
+ufunc_shichi_data[3] = &ufunc_shichi_ptr[2*3]
+shichi = np.PyUFunc_FromFuncAndData(ufunc_shichi_loops, ufunc_shichi_data, ufunc_shichi_types, 4, 1, 2, 0, "shichi", ufunc_shichi_doc, 0)
+
+cdef np.PyUFuncGenericFunction ufunc_sici_loops[4]
+cdef void *ufunc_sici_ptr[8]
+cdef void *ufunc_sici_data[4]
+cdef char ufunc_sici_types[12]
 cdef char *ufunc_sici_doc = (
-    "sici(x)\n"
+    "sici(x, out=None)\n"
     "\n"
-    "Sine and cosine integrals\n"
+    "Sine and cosine integrals.\n"
+    "\n"
+    "The sine integral is\n"
+    "\n"
+    ".. math::\n"
+    "\n"
+    "  \\int_0^x \\frac{\\sin{t}}{t}dt\n"
+    "\n"
+    "and the cosine integral is\n"
+    "\n"
+    ".. math::\n"
+    "\n"
+    "  \\gamma + \\log(x) + \\int_0^x \\frac{\\cos{t} - 1}{t}dt\n"
+    "\n"
+    "where :math:`\\gamma` is Euler's constant and :math:`\\log` is the\n"
+    "principle branch of the logarithm.\n"
+    "\n"
+    "Parameters\n"
+    "----------\n"
+    "x : array_like\n"
+    "    Real or complex points at which to compute the sine and cosine\n"
+    "    integrals.\n"
     "\n"
     "Returns\n"
     "-------\n"
-    "si\n"
-    "    ``integral(sin(t)/t, t=0..x)``\n"
-    "ci\n"
-    "    ``eul + ln x + integral((cos(t) - 1)/t, t=0..x)``\n"
-    "    where ``eul`` is Euler's constant.")
+    "si : ndarray\n"
+    "    Sine integral at ``x``\n"
+    "ci : ndarray\n"
+    "    Cosine integral at ``x``\n"
+    "\n"
+    "Notes\n"
+    "-----\n"
+    "For real arguments with ``x < 0``, ``ci`` is the real part of the\n"
+    "cosine integral. For such points ``ci(x)`` and ``ci(x + 0j)``\n"
+    "differ by a factor of ``1j*pi``.\n"
+    "\n"
+    "For real arguments the function is computed by calling Cephes'\n"
+    "[1]_ *sici* routine. For complex arguments the algorithm is based\n"
+    "on Mpmath's [2]_ *si* and *ci* routines.\n"
+    "\n"
+    "References\n"
+    "----------\n"
+    ".. [1] Cephes Mathematical Functions Library,\n"
+    "       http://www.netlib.org/cephes/index.html\n"
+    ".. [2] Fredrik Johansson and others.\n"
+    "       \"mpmath: a Python library for arbitrary-precision floating-point arithmetic\"\n"
+    "       (Version 0.19) http://mpmath.org/")
 ufunc_sici_loops[0] = <np.PyUFuncGenericFunction>loop_i_d_dd_As_f_ff
 ufunc_sici_loops[1] = <np.PyUFuncGenericFunction>loop_i_d_dd_As_d_dd
+ufunc_sici_loops[2] = <np.PyUFuncGenericFunction>loop_i_D_DD_As_F_FF
+ufunc_sici_loops[3] = <np.PyUFuncGenericFunction>loop_i_D_DD_As_D_DD
 ufunc_sici_types[0] = <char>NPY_FLOAT
 ufunc_sici_types[1] = <char>NPY_FLOAT
 ufunc_sici_types[2] = <char>NPY_FLOAT
 ufunc_sici_types[3] = <char>NPY_DOUBLE
 ufunc_sici_types[4] = <char>NPY_DOUBLE
 ufunc_sici_types[5] = <char>NPY_DOUBLE
+ufunc_sici_types[6] = <char>NPY_CFLOAT
+ufunc_sici_types[7] = <char>NPY_CFLOAT
+ufunc_sici_types[8] = <char>NPY_CFLOAT
+ufunc_sici_types[9] = <char>NPY_CDOUBLE
+ufunc_sici_types[10] = <char>NPY_CDOUBLE
+ufunc_sici_types[11] = <char>NPY_CDOUBLE
 ufunc_sici_ptr[2*0] = <void*>_func_sici
 ufunc_sici_ptr[2*0+1] = <void*>(<char*>"sici")
 ufunc_sici_ptr[2*1] = <void*>_func_sici
 ufunc_sici_ptr[2*1+1] = <void*>(<char*>"sici")
+ufunc_sici_ptr[2*2] = <void*>_func_csici
+ufunc_sici_ptr[2*2+1] = <void*>(<char*>"sici")
+ufunc_sici_ptr[2*3] = <void*>_func_csici
+ufunc_sici_ptr[2*3+1] = <void*>(<char*>"sici")
 ufunc_sici_data[0] = &ufunc_sici_ptr[2*0]
 ufunc_sici_data[1] = &ufunc_sici_ptr[2*1]
-sici = np.PyUFunc_FromFuncAndData(ufunc_sici_loops, ufunc_sici_data, ufunc_sici_types, 2, 1, 2, 0, "sici", ufunc_sici_doc, 0)
+ufunc_sici_data[2] = &ufunc_sici_ptr[2*2]
+ufunc_sici_data[3] = &ufunc_sici_ptr[2*3]
+sici = np.PyUFunc_FromFuncAndData(ufunc_sici_loops, ufunc_sici_data, ufunc_sici_types, 4, 1, 2, 0, "sici", ufunc_sici_doc, 0)
 
 cdef np.PyUFuncGenericFunction ufunc_sindg_loops[2]
 cdef void *ufunc_sindg_ptr[4]