@defun
def nsum(ctx, f, *intervals, **options):
    r"""
    Computes the sum

    .. math :: S = \sum_{k=a}^b f(k)

    where `(a, b)` = *interval*, and where `a = -\infty` and/or
    `b = \infty` are allowed, or more generally

    .. math :: S = \sum_{k_1=a_1}^{b_1} \cdots
                   \sum_{k_n=a_n}^{b_n} f(k_1,\ldots,k_n)

    if multiple intervals are given.

    Two examples of infinite series that can be summed by :func:`~mpmath.nsum`,
    where the first converges rapidly and the second converges slowly,
    are::

        >>> from mpmath import *
        >>> mp.dps = 15; mp.pretty = True
        >>> nsum(lambda n: 1/fac(n), [0, inf])
        2.71828182845905
        >>> nsum(lambda n: 1/n**2, [1, inf])
        1.64493406684823

    When appropriate, :func:`~mpmath.nsum` applies convergence acceleration to
    accurately estimate the sums of slowly convergent series. If the series is
    finite, :func:`~mpmath.nsum` currently does not attempt to perform any
    extrapolation, and simply calls :func:`~mpmath.fsum`.

    Multidimensional infinite series are reduced to a single-dimensional
    series over expanding hypercubes; if both infinite and finite dimensions
    are present, the finite ranges are moved innermost. For more advanced
    control over the summation order, use nested calls to :func:`~mpmath.nsum`,
    or manually rewrite the sum as a single-dimensional series.

    **Options**

    *tol*
        Desired maximum final error. Defaults roughly to the
        epsilon of the working precision.

    *method*
        Which summation algorithm to use (described below).
        Default: ``'richardson+shanks'``.

    *maxterms*
        Cancel after at most this many terms. Default: 10*dps.

    *steps*
        An iterable giving the number of terms to add between
        each extrapolation attempt. The default sequence is
        [10, 20, 30, 40, ...]. For example, if you know that
        approximately 100 terms will be required, efficiency might be
        improved by setting this to [100, 10]. Then the first
        extrapolation will be performed after 100 terms, the second
        after 110, etc.

    *verbose*
        Print details about progress.

    *ignore*
        If enabled, any term that raises ``ArithmeticError``
        or ``ValueError`` (e.g. through division by zero) is replaced
        by a zero. This is convenient for lattice sums with
        a singular term near the origin.

    **Methods**

    Unfortunately, an algorithm that can efficiently sum any infinite
    series does not exist. :func:`~mpmath.nsum` implements several different
    algorithms that each work well in different cases. The *method*
    keyword argument selects a method.

    The default method is ``'r+s'``, i.e. both Richardson extrapolation
    and Shanks transformation is attempted. A slower method that
    handles more cases is ``'r+s+e'``. For very high precision
    summation, or if the summation needs to be fast (for example if
    multiple sums need to be evaluated), it is a good idea to
    investigate which one method works best and only use that.

    ``'richardson'`` / ``'r'``:
        Uses Richardson extrapolation. Provides useful extrapolation
        when `f(k) \sim P(k)/Q(k)` or when `f(k) \sim (-1)^k P(k)/Q(k)`
        for polynomials `P` and `Q`. See :func:`~mpmath.richardson` for
        additional information.

    ``'shanks'`` / ``'s'``:
        Uses Shanks transformation. Typically provides useful
        extrapolation when `f(k) \sim c^k` or when successive terms
        alternate signs. Is able to sum some divergent series.
        See :func:`~mpmath.shanks` for additional information.

    ``'levin'`` / ``'l'``:
        Uses the Levin transformation. It performs better than the Shanks
        transformation for logarithmic convergent or alternating divergent
        series. The ``'levin_variant'``-keyword selects the variant. Valid
        choices are "u", "t", "v" and "all" whereby "all" uses all three
        u,t and v simultanously (This is good for performance comparison in
        conjunction with "verbose=True"). Instead of the Levin transform one can
        also use the Sidi-S transform by selecting the method ``'sidi'``.
        See :func:`~mpmath.levin` for additional details.

    ``'alternating'`` / ``'a'``:
        This is the convergence acceleration of alternating series developped
        by Cohen, Villegras and Zagier.
        See :func:`~mpmath.cohen_alt` for additional details.

    ``'euler-maclaurin'`` / ``'e'``:
        Uses the Euler-Maclaurin summation formula to approximate
        the remainder sum by an integral. This requires high-order
        numerical derivatives and numerical integration. The advantage
        of this algorithm is that it works regardless of the
        decay rate of `f`, as long as `f` is sufficiently smooth.
        See :func:`~mpmath.sumem` for additional information.

    ``'direct'`` / ``'d'``:
        Does not perform any extrapolation. This can be used
        (and should only be used for) rapidly convergent series.
        The summation automatically stops when the terms
        decrease below the target tolerance.

    **Basic examples**

    A finite sum::

        >>> nsum(lambda k: 1/k, [1, 6])
        2.45

    Summation of a series going to negative infinity and a doubly
    infinite series::

        >>> nsum(lambda k: 1/k**2, [-inf, -1])
        1.64493406684823
        >>> nsum(lambda k: 1/(1+k**2), [-inf, inf])
        3.15334809493716

    :func:`~mpmath.nsum` handles sums of complex numbers::

        >>> nsum(lambda k: (0.5+0.25j)**k, [0, inf])
        (1.6 + 0.8j)

    The following sum converges very rapidly, so it is most
    efficient to sum it by disabling convergence acceleration::

        >>> mp.dps = 1000
        >>> a = nsum(lambda k: -(-1)**k * k**2 / fac(2*k), [1, inf],
        ...     method='direct')
        >>> b = (cos(1)+sin(1))/4
        >>> abs(a-b) < mpf('1e-998')
        True

    **Examples with Richardson extrapolation**

    Richardson extrapolation works well for sums over rational
    functions, as well as their alternating counterparts::

        >>> mp.dps = 50
        >>> nsum(lambda k: 1 / k**3, [1, inf],
        ...     method='richardson')
        1.2020569031595942853997381615114499907649862923405
        >>> zeta(3)
        1.2020569031595942853997381615114499907649862923405

        >>> nsum(lambda n: (n + 3)/(n**3 + n**2), [1, inf],
        ...     method='richardson')
        2.9348022005446793094172454999380755676568497036204
        >>> pi**2/2-2
        2.9348022005446793094172454999380755676568497036204

        >>> nsum(lambda k: (-1)**k / k**3, [1, inf],
        ...     method='richardson')
        -0.90154267736969571404980362113358749307373971925537
        >>> -3*zeta(3)/4
        -0.90154267736969571404980362113358749307373971925538

    **Examples with Shanks transformation**

    The Shanks transformation works well for geometric series
    and typically provides excellent acceleration for Taylor
    series near the border of their disk of convergence.
    Here we apply it to a series for `\log(2)`, which can be
    seen as the Taylor series for `\log(1+x)` with `x = 1`::

        >>> nsum(lambda k: -(-1)**k/k, [1, inf],
        ...     method='shanks')
        0.69314718055994530941723212145817656807550013436025
        >>> log(2)
        0.69314718055994530941723212145817656807550013436025

    Here we apply it to a slowly convergent geometric series::

        >>> nsum(lambda k: mpf('0.995')**k, [0, inf],
        ...     method='shanks')
        200.0

    Finally, Shanks' method works very well for alternating series
    where `f(k) = (-1)^k g(k)`, and often does so regardless of
    the exact decay rate of `g(k)`::

        >>> mp.dps = 15
        >>> nsum(lambda k: (-1)**(k+1) / k**1.5, [1, inf],
        ...     method='shanks')
        0.765147024625408
        >>> (2-sqrt(2))*zeta(1.5)/2
        0.765147024625408

    The following slowly convergent alternating series has no known
    closed-form value. Evaluating the sum a second time at higher
    precision indicates that the value is probably correct::

        >>> nsum(lambda k: (-1)**k / log(k), [2, inf],
        ...     method='shanks')
        0.924299897222939
        >>> mp.dps = 30
        >>> nsum(lambda k: (-1)**k / log(k), [2, inf],
        ...     method='shanks')
        0.92429989722293885595957018136

    **Examples with Levin transformation**

    The following example calculates Euler's constant as the constant term in
    the Laurent expansion of zeta(s) at s=1. This sum converges extremly slow
    because of the logarithmic convergence behaviour of the Dirichlet series
    for zeta.

      >>> mp.dps = 30
      >>> z = mp.mpf(10) ** (-10)
      >>> a = mp.nsum(lambda n: n**(-(1+z)), [1, mp.inf], method = "levin") - 1 / z
      >>> print(mp.chop(a - mp.euler, tol = 1e-10))
      0.0

    Now we sum the zeta function outside its range of convergence
    (attention: This does not work at the negative integers!):

      >>> mp.dps = 15
      >>> w = mp.nsum(lambda n: n ** (2 + 3j), [1, mp.inf], method = "levin", levin_variant = "v")
      >>> print(mp.chop(w - mp.zeta(-2-3j)))
      0.0

    The next example resummates an asymptotic series expansion of an integral
    related to the exponential integral.

      >>> mp.dps = 15
      >>> z = mp.mpf(10)
      >>> # exact = mp.quad(lambda x: mp.exp(-x)/(1+x/z),[0,mp.inf])
      >>> exact = z * mp.exp(z) * mp.expint(1,z) # this is the symbolic expression for the integral
      >>> w = mp.nsum(lambda n: (-1) ** n * mp.fac(n) * z ** (-n), [0, mp.inf], method = "sidi", levin_variant = "t")
      >>> print(mp.chop(w - exact))
      0.0

    Following highly divergent asymptotic expansion needs some care. Firstly we
    need copious amount of working precision. Secondly the stepsize must not be
    chosen to large, otherwise nsum may miss the point where the Levin transform
    converges and reach the point where only numerical garbage is produced due to
    numerical cancellation.

      >>> mp.dps = 15
      >>> z = mp.mpf(2)
      >>> # exact = mp.quad(lambda x: mp.exp( -x * x / 2 - z * x ** 4), [0,mp.inf]) * 2 / mp.sqrt(2 * mp.pi)
      >>> exact = mp.exp(mp.one / (32 * z)) * mp.besselk(mp.one / 4, mp.one / (32 * z)) / (4 * mp.sqrt(z * mp.pi)) # this is the symbolic expression for the integral
      >>> w = mp.nsum(lambda n: (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n)),
      ...   [0, mp.inf], method = "levin", levin_variant = "t", workprec = 8*mp.prec, steps = [2] + [1 for x in xrange(1000)])
      >>> print(mp.chop(w - exact))
      0.0

    The hypergeoemtric function can also be summed outside its range of convergence:

      >>> mp.dps = 15
      >>> z = 2 + 1j
      >>> exact = mp.hyp2f1(2 / mp.mpf(3), 4 / mp.mpf(3), 1 / mp.mpf(3), z)
      >>> f = lambda n: mp.rf(2 / mp.mpf(3), n) * mp.rf(4 / mp.mpf(3), n) * z**n / (mp.rf(1 / mp.mpf(3), n) * mp.fac(n))
      >>> v = mp.nsum(f, [0, mp.inf], method = "levin", steps = [10 for x in xrange(1000)])
      >>> print(mp.chop(exact-v))
      0.0

    **Examples with Cohen's alternating series resummation**

      The next example sums the alternating zeta function:

      >>> v = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "a")
      >>> print(mp.chop(v - mp.log(2)))
      0.0

      The derivate of the alternating zeta function outside its range of
      convergence:

      >>> v = mp.nsum(lambda n: (-1)**n * mp.log(n) * n, [1, mp.inf], method = "a")
      >>> print(mp.chop(v - mp.diff(lambda s: mp.altzeta(s), -1)))
      0.0

    **Examples with Euler-Maclaurin summation**

    The sum in the following example has the wrong rate of convergence
    for either Richardson or Shanks to be effective.

        >>> f = lambda k: log(k)/k**2.5
        >>> mp.dps = 15
        >>> nsum(f, [1, inf], method='euler-maclaurin')
        0.38734195032621
        >>> -diff(zeta, 2.5)
        0.38734195032621

    Increasing ``steps`` improves speed at higher precision::

        >>> mp.dps = 50
        >>> nsum(f, [1, inf], method='euler-maclaurin', steps=[250])
        0.38734195032620997271199237593105101319948228874688
        >>> -diff(zeta, 2.5)
        0.38734195032620997271199237593105101319948228874688

    **Divergent series**

    The Shanks transformation is able to sum some *divergent*
    series. In particular, it is often able to sum Taylor series
    beyond their radius of convergence (this is due to a relation
    between the Shanks transformation and Pade approximations;
    see :func:`~mpmath.pade` for an alternative way to evaluate divergent
    Taylor series). Furthermore the Levin-transform examples above
    contain some divergent series resummation.

    Here we apply it to `\log(1+x)` far outside the region of
    convergence::

        >>> mp.dps = 50
        >>> nsum(lambda k: -(-9)**k/k, [1, inf],
        ...     method='shanks')
        2.3025850929940456840179914546843642076011014886288
        >>> log(10)
        2.3025850929940456840179914546843642076011014886288

    A particular type of divergent series that can be summed
    using the Shanks transformation is geometric series.
    The result is the same as using the closed-form formula
    for an infinite geometric series::

        >>> mp.dps = 15
        >>> for n in range(-8, 8):
        ...     if n == 1:
        ...         continue
        ...     print("%s %s %s" % (mpf(n), mpf(1)/(1-n),
        ...         nsum(lambda k: n**k, [0, inf], method='shanks')))
        ...
        -8.0 0.111111111111111 0.111111111111111
        -7.0 0.125 0.125
        -6.0 0.142857142857143 0.142857142857143
        -5.0 0.166666666666667 0.166666666666667
        -4.0 0.2 0.2
        -3.0 0.25 0.25
        -2.0 0.333333333333333 0.333333333333333
        -1.0 0.5 0.5
        0.0 1.0 1.0
        2.0 -1.0 -1.0
        3.0 -0.5 -0.5
        4.0 -0.333333333333333 -0.333333333333333
        5.0 -0.25 -0.25
        6.0 -0.2 -0.2
        7.0 -0.166666666666667 -0.166666666666667

    **Multidimensional sums**

    Any combination of finite and infinite ranges is allowed for the
    summation indices::

        >>> mp.dps = 15
        >>> nsum(lambda x,y: x+y, [2,3], [4,5])
        28.0
        >>> nsum(lambda x,y: x/2**y, [1,3], [1,inf])
        6.0
        >>> nsum(lambda x,y: y/2**x, [1,inf], [1,3])
        6.0
        >>> nsum(lambda x,y,z: z/(2**x*2**y), [1,inf], [1,inf], [3,4])
        7.0
        >>> nsum(lambda x,y,z: y/(2**x*2**z), [1,inf], [3,4], [1,inf])
        7.0
        >>> nsum(lambda x,y,z: x/(2**z*2**y), [3,4], [1,inf], [1,inf])
        7.0

    Some nice examples of double series with analytic solutions or
    reductions to single-dimensional series (see [1])::

        >>> nsum(lambda m, n: 1/2**(m*n), [1,inf], [1,inf])
        1.60669515241529
        >>> nsum(lambda n: 1/(2**n-1), [1,inf])
        1.60669515241529

        >>> nsum(lambda i,j: (-1)**(i+j)/(i**2+j**2), [1,inf], [1,inf])
        0.278070510848213
        >>> pi*(pi-3*ln2)/12
        0.278070510848213

        >>> nsum(lambda i,j: (-1)**(i+j)/(i+j)**2, [1,inf], [1,inf])
        0.129319852864168
        >>> altzeta(2) - altzeta(1)
        0.129319852864168

        >>> nsum(lambda i,j: (-1)**(i+j)/(i+j)**3, [1,inf], [1,inf])
        0.0790756439455825
        >>> altzeta(3) - altzeta(2)
        0.0790756439455825

        >>> nsum(lambda m,n: m**2*n/(3**m*(n*3**m+m*3**n)),
        ...     [1,inf], [1,inf])
        0.28125
        >>> mpf(9)/32
        0.28125

        >>> nsum(lambda i,j: fac(i-1)*fac(j-1)/fac(i+j),
        ...     [1,inf], [1,inf], workprec=400)
        1.64493406684823
        >>> zeta(2)
        1.64493406684823

    A hard example of a multidimensional sum is the Madelung constant
    in three dimensions (see [2]). The defining sum converges very
    slowly and only conditionally, so :func:`~mpmath.nsum` is lucky to
    obtain an accurate value through convergence acceleration. The
    second evaluation below uses a much more efficient, rapidly
    convergent 2D sum::

        >>> nsum(lambda x,y,z: (-1)**(x+y+z)/(x*x+y*y+z*z)**0.5,
        ...     [-inf,inf], [-inf,inf], [-inf,inf], ignore=True)
        -1.74756459463318
        >>> nsum(lambda x,y: -12*pi*sech(0.5*pi * \
        ...     sqrt((2*x+1)**2+(2*y+1)**2))**2, [0,inf], [0,inf])
        -1.74756459463318

    Another example of a lattice sum in 2D::

        >>> nsum(lambda x,y: (-1)**(x+y) / (x**2+y**2), [-inf,inf],
        ...     [-inf,inf], ignore=True)
        -2.1775860903036
        >>> -pi*ln2
        -2.1775860903036

    An example of an Eisenstein series::

        >>> nsum(lambda m,n: (m+n*1j)**(-4), [-inf,inf], [-inf,inf],
        ...     ignore=True)
        (3.1512120021539 + 0.0j)

    **References**

    1. [Weisstein]_ http://mathworld.wolfram.com/DoubleSeries.html,
    2. [Weisstein]_ http://mathworld.wolfram.com/MadelungConstants.html

    """
    infinite, g = standardize(ctx, f, intervals, options)
    if not infinite:
        return +g()

    def update(partial_sums, indices):
        if partial_sums:
            psum = partial_sums[-1]
        else:
            psum = ctx.zero
        for k in indices:
            psum = psum + g(ctx.mpf(k))
            partial_sums.append(psum)

    prec = ctx.prec

    def emfun(point, tol):
        workprec = ctx.prec
        ctx.prec = prec * 10
        v = ctx.sumem(g, [point, ctx.inf], tol, error=1)
        ctx.prec = workprec
        return v

    return +ctx.adaptive_extrapolation(update, emfun, options)
dataset/ETHPy150Open fredrik-johansson/mpmath/mpmath/calculus/extrapolation.py/nsum/operator `+` -> `*`