matical Tables. New York: Dover, 1972.
           (See Section 5.2.)
    .. [2] Cephes Mathematical Functions Library,
           http://www.netlib.org/cephes/
    .. [3] Fredrik Johansson and others.
           "mpmath: a Python library for arbitrary-precision floating-point
           arithmetic" (Version 0.19) http://mpmath.org/

    Examples
    --------
    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> from scipy.special import sici, exp1

    `sici` accepts real or complex input:

    >>> sici(2.5)
    (1.7785201734438267, 0.2858711963653835)
    >>> sici(2.5 + 3j)
    ((4.505735874563953+0.06863305018999577j),
    (0.0793644206906966-2.935510262937543j))

    For z in the right half plane, the sine and cosine integrals are
    related to the exponential integral E1 (implemented in SciPy as
    `scipy.special.exp1`) by

    * Si(z) = (E1(i*z) - E1(-i*z))/2i + pi/2
    * Ci(z) = -(E1(i*z) + E1(-i*z))/2

    See [1]_ (equations 5.2.21 and 5.2.23).

    We can verify these relations:

    >>> z = 2 - 3j
    >>> sici(z)
    ((4.54751388956229-1.3991965806460565j),
    (1.408292501520851+2.9836177420296055j))

    >>> (exp1(1j*z) - exp1(-1j*z))/2j + np.pi/2  # Same as sine integral
    (4.54751388956229-1.3991965806460565j)

    >>> -(exp1(1j*z) + exp1(-1j*z))/2            # Same as cosine integral
    (1.408292501520851+2.9836177420296055j)

    Plot the functions evaluated on the real axis; the dotted horizontal
    lines are at pi/2 and -pi/2:

    >>> x = np.linspace(-16, 16, 150)
    >>> si, ci = sici(x)

    >>> fig, ax = plt.subplots()
    >>> ax.plot(x, si, label='Si(x)')
    >>> ax.plot(x, ci, '--', label='Ci(x)')
    >>> ax.legend(shadow=True, framealpha=1, loc='upper left')
    >>> ax.set_xlabel('x')
    >>> ax.set_title('Sine and Cosine Integrals')
    >>> ax.axhline(np.pi/2, linestyle=':', alpha=0.5, color='k')
    >>> ax.axhline(-np.pi/2, linestyle=':', alpha=0.5, color='k')
    >>> ax.grid(True)
    >>> plt.show()

    