loosening `tolerances`.
          Alternatively, the callable may not be unimodal, or the limits of
          summation may be too far from the function maximum. Consider
          increasing `maxterms` or breaking the sum into pieces.

    sum : float array
        An estimate of the sum.
    error : float array
        An estimate of the absolute error, assuming all terms are non-negative,
        the function is computed exactly, and direct summation is accurate to
        the precision of the result dtype.
    nfev : int array
        The number of points at which `f` was evaluated.

See Also
--------
mpmath.nsum

Notes
-----
The method implemented for infinite summation is related to the integral
test for convergence of an infinite series: assuming `step` size 1 for
simplicity of exposition, the sum of a monotone decreasing function is bounded by

.. math::

    \int_u^\infty f(x) dx \leq \sum_{k=u}^\infty f(k) \leq \int_u^\infty f(x) dx + f(u)

Let :math:`a` represent  `a`, :math:`n` represent `maxterms`, :math:`\epsilon_a`
represent `atol`, and :math:`\epsilon_r` represent `rtol`.
The implementation first evaluates the integral :math:`S_l=\int_a^\infty f(x) dx`
as a lower bound of the infinite sum. Then, it seeks a value :math:`c > a` such
that :math:`f(c) < \epsilon_a + S_l \epsilon_r`, if it exists; otherwise,
let :math:`c = a + n`. Then the infinite sum is approximated as

.. math::

    \sum_{k=a}^{c-1} f(k) + \int_c^\infty f(x) dx + f(c)/2

and the reported error is :math:`f(c)/2` plus the error estimate of
numerical integration. Note that the integral approximations may require
evaluation of the function at points besides those that appear in the sum,
so `f` must be a continuous and monotonically decreasing function defined
for all reals within the integration interval. However, due to the nature
of the integral approximation, the shape of the function between points
that appear in the sum has little effect. If there is not a natural
extension of the function to all reals, consider using linear interpolation,
which is easy to evaluate and preserves monotonicity.

The approach described above is generalized for non-unit
`step` and finite `b` that is too large for direct evaluation of the sum,
i.e. ``b - a + 1 > maxterms``. It is further generalized to unimodal
functions by directly summing terms surrounding the maximum.
This strategy may fail:

- If the left limit is finite and the maximum is far from it.
- If the right limit is finite and the maximum is far from it.
- If both limits are finite and the maximum is far from the origin.

In these cases, accuracy may be poor, and `nsum` may return status code ``4``.

Although the callable `f` must be non-negative and unimodal,
`nsum` can be used to evaluate more general forms of series. For instance, to
evaluate an alternating series, pass a callable that returns the difference
between pairs of adjacent terms, and adjust `step` accordingly. See Examples.

References
----------
.. [1] Wikipedia. "Integral test for convergence."
       https://en.wikipedia.org/wiki/Integral_test_for_convergence

Examples
--------
Compute the infinite sum of the reciprocals of squared integers.

>>> import numpy as np
>>> from scipy.integrate import nsum
>>> res = nsum(lambda k: 1/k**2, 1, np.inf)
>>> ref = np.pi**2/6  # true value
>>> res.error  # estimated error
np.float64(7.448762306416137e-09)
>>> (res.sum - ref)/ref  # true error
np.float64(-1.839871898894426e-13)
>>> res.nfev  # number of points at which callable was evaluated
np.int32(8561)

Compute the infinite sums of the reciprocals of integers raised to powers ``p``,
where ``p`` is an array.

>>> from scipy import special
>>> p = np.arange(3, 10)
>>> res = nsum(lambda k, p: 1/k**p, 1, np.inf, maxterms=1e3, args=(p,))
>>> ref = special.zeta(p, 1)
>>> np.allclose(res.sum, ref)
True

Evaluate the alternating harmonic series.

>>> res = nsum(lambda x: 1/x - 1/(x+1), 1, np.inf, step=2)
>>> res.sum, res.sum - np.log(2)  # result, difference vs analytical sum
(np.float64(0.6931471805598691), np.float64(-7.616129948928574e-14))