code: 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)
  
imports: from ..libmp.backend import xrange
  from .calculus import defun
label_desc: Use image node_id 8 for calling a global function with example usage: nsum(ctx, f) and returns: unknown, unknown, v