

@defun
def nsum(ctx, f, *intervals, **options):
    '\n    Computes the sum\n\n    .. math :: S = \\sum_{k=a}^b f(k)\n\n    where `(a, b)` = *interval*, and where `a = -\\infty` and/or\n    `b = \\infty` are allowed, or more generally\n\n    .. math :: S = \\sum_{k_1=a_1}^{b_1} \\cdots\n                   \\sum_{k_n=a_n}^{b_n} f(k_1,\\ldots,k_n)\n\n    if multiple intervals are given.\n\n    Two examples of infinite series that can be summed by :func:`~mpmath.nsum`,\n    where the first converges rapidly and the second converges slowly,\n    are::\n\n        >>> from mpmath import *\n        >>> mp.dps = 15; mp.pretty = True\n        >>> nsum(lambda n: 1/fac(n), [0, inf])\n        2.71828182845905\n        >>> nsum(lambda n: 1/n**2, [1, inf])\n        1.64493406684823\n\n    When appropriate, :func:`~mpmath.nsum` applies convergence acceleration to\n    accurately estimate the sums of slowly convergent series. If the series is\n    finite, :func:`~mpmath.nsum` currently does not attempt to perform any\n    extrapolation, and simply calls :func:`~mpmath.fsum`.\n\n    Multidimensional infinite series are reduced to a single-dimensional\n    series over expanding hypercubes; if both infinite and finite dimensions\n    are present, the finite ranges are moved innermost. For more advanced\n    control over the summation order, use nested calls to :func:`~mpmath.nsum`,\n    or manually rewrite the sum as a single-dimensional series.\n\n    **Options**\n\n    *tol*\n        Desired maximum final error. Defaults roughly to the\n        epsilon of the working precision.\n\n    *method*\n        Which summation algorithm to use (described below).\n        Default: ``\'richardson+shanks\'``.\n\n    *maxterms*\n        Cancel after at most this many terms. Default: 10*dps.\n\n    *steps*\n        An iterable giving the number of terms to add between\n        each extrapolation attempt. The default sequence is\n        [10, 20, 30, 40, ...]. For example, if you know that\n        approximately 100 terms will be required, efficiency might be\n        improved by setting this to [100, 10]. Then the first\n        extrapolation will be performed after 100 terms, the second\n        after 110, etc.\n\n    *verbose*\n        Print details about progress.\n\n    *ignore*\n        If enabled, any term that raises ``ArithmeticError``\n        or ``ValueError`` (e.g. through division by zero) is replaced\n        by a zero. This is convenient for lattice sums with\n        a singular term near the origin.\n\n    **Methods**\n\n    Unfortunately, an algorithm that can efficiently sum any infinite\n    series does not exist. :func:`~mpmath.nsum` implements several different\n    algorithms that each work well in different cases. The *method*\n    keyword argument selects a method.\n\n    The default method is ``\'r+s\'``, i.e. both Richardson extrapolation\n    and Shanks transformation is attempted. A slower method that\n    handles more cases is ``\'r+s+e\'``. For very high precision\n    summation, or if the summation needs to be fast (for example if\n    multiple sums need to be evaluated), it is a good idea to\n    investigate which one method works best and only use that.\n\n    ``\'richardson\'`` / ``\'r\'``:\n        Uses Richardson extrapolation. Provides useful extrapolation\n        when `f(k) \\sim P(k)/Q(k)` or when `f(k) \\sim (-1)^k P(k)/Q(k)`\n        for polynomials `P` and `Q`. See :func:`~mpmath.richardson` for\n        additional information.\n\n    ``\'shanks\'`` / ``\'s\'``:\n        Uses Shanks transformation. Typically provides useful\n        extrapolation when `f(k) \\sim c^k` or when successive terms\n        alternate signs. Is able to sum some divergent series.\n        See :func:`~mpmath.shanks` for additional information.\n\n    ``\'levin\'`` / ``\'l\'``:\n        Uses the Levin transformation. It performs better than the Shanks\n        transformation for logarithmic convergent or alternating divergent\n        series. The ``\'levin_variant\'``-keyword selects the variant. Valid\n        choices are "u", "t", "v" and "all" whereby "all" uses all three\n        u,t and v simultanously (This is good for performance comparison in\n        conjunction with "verbose=True"). Instead of the Levin transform one can\n        also use the Sidi-S transform by selecting the method ``\'sidi\'``.\n        See :func:`~mpmath.levin` for additional details.\n\n    ``\'alternating\'`` / ``\'a\'``:\n        This is the convergence acceleration of alternating series developped\n        by Cohen, Villegras and Zagier.\n        See :func:`~mpmath.cohen_alt` for additional details.\n\n    ``\'euler-maclaurin\'`` / ``\'e\'``:\n        Uses the Euler-Maclaurin summation formula to approximate\n        the remainder sum by an integral. This requires high-order\n        numerical derivatives and numerical integration. The advantage\n        of this algorithm is that it works regardless of the\n        decay rate of `f`, as long as `f` is sufficiently smooth.\n        See :func:`~mpmath.sumem` for additional information.\n\n    ``\'direct\'`` / ``\'d\'``:\n        Does not perform any extrapolation. This can be used\n        (and should only be used for) rapidly convergent series.\n        The summation automatically stops when the terms\n        decrease below the target tolerance.\n\n    **Basic examples**\n\n    A finite sum::\n\n        >>> nsum(lambda k: 1/k, [1, 6])\n        2.45\n\n    Summation of a series going to negative infinity and a doubly\n    infinite series::\n\n        >>> nsum(lambda k: 1/k**2, [-inf, -1])\n        1.64493406684823\n        >>> nsum(lambda k: 1/(1+k**2), [-inf, inf])\n        3.15334809493716\n\n    :func:`~mpmath.nsum` handles sums of complex numbers::\n\n        >>> nsum(lambda k: (0.5+0.25j)**k, [0, inf])\n        (1.6 + 0.8j)\n\n    The following sum converges very rapidly, so it is most\n    efficient to sum it by disabling convergence acceleration::\n\n        >>> mp.dps = 1000\n        >>> a = nsum(lambda k: -(-1)**k * k**2 / fac(2*k), [1, inf],\n        ...     method=\'direct\')\n        >>> b = (cos(1)+sin(1))/4\n        >>> abs(a-b) < mpf(\'1e-998\')\n        True\n\n    **Examples with Richardson extrapolation**\n\n    Richardson extrapolation works well for sums over rational\n    functions, as well as their alternating counterparts::\n\n        >>> mp.dps = 50\n        >>> nsum(lambda k: 1 / k**3, [1, inf],\n        ...     method=\'richardson\')\n        1.2020569031595942853997381615114499907649862923405\n        >>> zeta(3)\n        1.2020569031595942853997381615114499907649862923405\n\n        >>> nsum(lambda n: (n + 3)/(n**3 + n**2), [1, inf],\n        ...     method=\'richardson\')\n        2.9348022005446793094172454999380755676568497036204\n        >>> pi**2/2-2\n        2.9348022005446793094172454999380755676568497036204\n\n        >>> nsum(lambda k: (-1)**k / k**3, [1, inf],\n        ...     method=\'richardson\')\n        -0.90154267736969571404980362113358749307373971925537\n        >>> -3*zeta(3)/4\n        -0.90154267736969571404980362113358749307373971925538\n\n    **Examples with Shanks transformation**\n\n    The Shanks transformation works well for geometric series\n    and typically provides excellent acceleration for Taylor\n    series near the border of their disk of convergence.\n    Here we apply it to a series for `\\log(2)`, which can be\n    seen as the Taylor series for `\\log(1+x)` with `x = 1`::\n\n        >>> nsum(lambda k: -(-1)**k/k, [1, inf],\n        ...     method=\'shanks\')\n        0.69314718055994530941723212145817656807550013436025\n        >>> log(2)\n        0.69314718055994530941723212145817656807550013436025\n\n    Here we apply it to a slowly convergent geometric series::\n\n        >>> nsum(lambda k: mpf(\'0.995\')**k, [0, inf],\n        ...     method=\'shanks\')\n        200.0\n\n    Finally, Shanks\' method works very well for alternating series\n    where `f(k) = (-1)^k g(k)`, and often does so regardless of\n    the exact decay rate of `g(k)`::\n\n        >>> mp.dps = 15\n        >>> nsum(lambda k: (-1)**(k+1) / k**1.5, [1, inf],\n        ...     method=\'shanks\')\n        0.765147024625408\n        >>> (2-sqrt(2))*zeta(1.5)/2\n        0.765147024625408\n\n    The following slowly convergent alternating series has no known\n    closed-form value. Evaluating the sum a second time at higher\n    precision indicates that the value is probably correct::\n\n        >>> nsum(lambda k: (-1)**k / log(k), [2, inf],\n        ...     method=\'shanks\')\n        0.924299897222939\n        >>> mp.dps = 30\n        >>> nsum(lambda k: (-1)**k / log(k), [2, inf],\n        ...     method=\'shanks\')\n        0.92429989722293885595957018136\n\n    **Examples with Levin transformation**\n\n    The following example calculates Euler\'s constant as the constant term in\n    the Laurent expansion of zeta(s) at s=1. This sum converges extremly slow\n    because of the logarithmic convergence behaviour of the Dirichlet series\n    for zeta.\n\n      >>> mp.dps = 30\n      >>> z = mp.mpf(10) ** (-10)\n      >>> a = mp.nsum(lambda n: n**(-(1+z)), [1, mp.inf], method = "levin") - 1 / z\n      >>> print(mp.chop(a - mp.euler, tol = 1e-10))\n      0.0\n\n    Now we sum the zeta function outside its range of convergence\n    (attention: This does not work at the negative integers!):\n\n      >>> mp.dps = 15\n      >>> w = mp.nsum(lambda n: n ** (2 + 3j), [1, mp.inf], method = "levin", levin_variant = "v")\n      >>> print(mp.chop(w - mp.zeta(-2-3j)))\n      0.0\n\n    The next example resummates an asymptotic series expansion of an integral\n    related to the exponential integral.\n\n      >>> mp.dps = 15\n      >>> z = mp.mpf(10)\n      >>> # exact = mp.quad(lambda x: mp.exp(-x)/(1+x/z),[0,mp.inf])\n      >>> exact = z * mp.exp(z) * mp.expint(1,z) # this is the symbolic expression for the integral\n      >>> w = mp.nsum(lambda n: (-1) ** n * mp.fac(n) * z ** (-n), [0, mp.inf], method = "sidi", levin_variant = "t")\n      >>> print(mp.chop(w - exact))\n      0.0\n\n    Following highly divergent asymptotic expansion needs some care. Firstly we\n    need copious amount of working precision. Secondly the stepsize must not be\n    chosen to large, otherwise nsum may miss the point where the Levin transform\n    converges and reach the point where only numerical garbage is produced due to\n    numerical cancellation.\n\n      >>> mp.dps = 15\n      >>> z = mp.mpf(2)\n      >>> # exact = mp.quad(lambda x: mp.exp( -x * x / 2 - z * x ** 4), [0,mp.inf]) * 2 / mp.sqrt(2 * mp.pi)\n      >>> 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\n      >>> w = mp.nsum(lambda n: (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n)),\n      ...   [0, mp.inf], method = "levin", levin_variant = "t", workprec = 8*mp.prec, steps = [2] + [1 for x in xrange(1000)])\n      >>> print(mp.chop(w - exact))\n      0.0\n\n    The hypergeoemtric function can also be summed outside its range of convergence:\n\n      >>> mp.dps = 15\n      >>> z = 2 + 1j\n      >>> exact = mp.hyp2f1(2 / mp.mpf(3), 4 / mp.mpf(3), 1 / mp.mpf(3), z)\n      >>> 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))\n      >>> v = mp.nsum(f, [0, mp.inf], method = "levin", steps = [10 for x in xrange(1000)])\n      >>> print(mp.chop(exact-v))\n      0.0\n\n    **Examples with Cohen\'s alternating series resummation**\n\n      The next example sums the alternating zeta function:\n\n      >>> v = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "a")\n      >>> print(mp.chop(v - mp.log(2)))\n      0.0\n\n      The derivate of the alternating zeta function outside its range of\n      convergence:\n\n      >>> v = mp.nsum(lambda n: (-1)**n * mp.log(n) * n, [1, mp.inf], method = "a")\n      >>> print(mp.chop(v - mp.diff(lambda s: mp.altzeta(s), -1)))\n      0.0\n\n    **Examples with Euler-Maclaurin summation**\n\n    The sum in the following example has the wrong rate of convergence\n    for either Richardson or Shanks to be effective.\n\n        >>> f = lambda k: log(k)/k**2.5\n        >>> mp.dps = 15\n        >>> nsum(f, [1, inf], method=\'euler-maclaurin\')\n        0.38734195032621\n        >>> -diff(zeta, 2.5)\n        0.38734195032621\n\n    Increasing ``steps`` improves speed at higher precision::\n\n        >>> mp.dps = 50\n        >>> nsum(f, [1, inf], method=\'euler-maclaurin\', steps=[250])\n        0.38734195032620997271199237593105101319948228874688\n        >>> -diff(zeta, 2.5)\n        0.38734195032620997271199237593105101319948228874688\n\n    **Divergent series**\n\n    The Shanks transformation is able to sum some *divergent*\n    series. In particular, it is often able to sum Taylor series\n    beyond their radius of convergence (this is due to a relation\n    between the Shanks transformation and Pade approximations;\n    see :func:`~mpmath.pade` for an alternative way to evaluate divergent\n    Taylor series). Furthermore the Levin-transform examples above\n    contain some divergent series resummation.\n\n    Here we apply it to `\\log(1+x)` far outside the region of\n    convergence::\n\n        >>> mp.dps = 50\n        >>> nsum(lambda k: -(-9)**k/k, [1, inf],\n        ...     method=\'shanks\')\n        2.3025850929940456840179914546843642076011014886288\n        >>> log(10)\n        2.3025850929940456840179914546843642076011014886288\n\n    A particular type of divergent series that can be summed\n    using the Shanks transformation is geometric series.\n    The result is the same as using the closed-form formula\n    for an infinite geometric series::\n\n        >>> mp.dps = 15\n        >>> for n in range(-8, 8):\n        ...     if n == 1:\n        ...         continue\n        ...     print("%s %s %s" % (mpf(n), mpf(1)/(1-n),\n        ...         nsum(lambda k: n**k, [0, inf], method=\'shanks\')))\n        ...\n        -8.0 0.111111111111111 0.111111111111111\n        -7.0 0.125 0.125\n        -6.0 0.142857142857143 0.142857142857143\n        -5.0 0.166666666666667 0.166666666666667\n        -4.0 0.2 0.2\n        -3.0 0.25 0.25\n        -2.0 0.333333333333333 0.333333333333333\n        -1.0 0.5 0.5\n        0.0 1.0 1.0\n        2.0 -1.0 -1.0\n        3.0 -0.5 -0.5\n        4.0 -0.333333333333333 -0.333333333333333\n        5.0 -0.25 -0.25\n        6.0 -0.2 -0.2\n        7.0 -0.166666666666667 -0.166666666666667\n\n    **Multidimensional sums**\n\n    Any combination of finite and infinite ranges is allowed for the\n    summation indices::\n\n        >>> mp.dps = 15\n        >>> nsum(lambda x,y: x+y, [2,3], [4,5])\n        28.0\n        >>> nsum(lambda x,y: x/2**y, [1,3], [1,inf])\n        6.0\n        >>> nsum(lambda x,y: y/2**x, [1,inf], [1,3])\n        6.0\n        >>> nsum(lambda x,y,z: z/(2**x*2**y), [1,inf], [1,inf], [3,4])\n        7.0\n        >>> nsum(lambda x,y,z: y/(2**x*2**z), [1,inf], [3,4], [1,inf])\n        7.0\n        >>> nsum(lambda x,y,z: x/(2**z*2**y), [3,4], [1,inf], [1,inf])\n        7.0\n\n    Some nice examples of double series with analytic solutions or\n    reductions to single-dimensional series (see [1])::\n\n        >>> nsum(lambda m, n: 1/2**(m*n), [1,inf], [1,inf])\n        1.60669515241529\n        >>> nsum(lambda n: 1/(2**n-1), [1,inf])\n        1.60669515241529\n\n        >>> nsum(lambda i,j: (-1)**(i+j)/(i**2+j**2), [1,inf], [1,inf])\n        0.278070510848213\n        >>> pi*(pi-3*ln2)/12\n        0.278070510848213\n\n        >>> nsum(lambda i,j: (-1)**(i+j)/(i+j)**2, [1,inf], [1,inf])\n        0.129319852864168\n        >>> altzeta(2) - altzeta(1)\n        0.129319852864168\n\n        >>> nsum(lambda i,j: (-1)**(i+j)/(i+j)**3, [1,inf], [1,inf])\n        0.0790756439455825\n        >>> altzeta(3) - altzeta(2)\n        0.0790756439455825\n\n        >>> nsum(lambda m,n: m**2*n/(3**m*(n*3**m+m*3**n)),\n        ...     [1,inf], [1,inf])\n        0.28125\n        >>> mpf(9)/32\n        0.28125\n\n        >>> nsum(lambda i,j: fac(i-1)*fac(j-1)/fac(i+j),\n        ...     [1,inf], [1,inf], workprec=400)\n        1.64493406684823\n        >>> zeta(2)\n        1.64493406684823\n\n    A hard example of a multidimensional sum is the Madelung constant\n    in three dimensions (see [2]). The defining sum converges very\n    slowly and only conditionally, so :func:`~mpmath.nsum` is lucky to\n    obtain an accurate value through convergence acceleration. The\n    second evaluation below uses a much more efficient, rapidly\n    convergent 2D sum::\n\n        >>> nsum(lambda x,y,z: (-1)**(x+y+z)/(x*x+y*y+z*z)**0.5,\n        ...     [-inf,inf], [-inf,inf], [-inf,inf], ignore=True)\n        -1.74756459463318\n        >>> nsum(lambda x,y: -12*pi*sech(0.5*pi * \\\n        ...     sqrt((2*x+1)**2+(2*y+1)**2))**2, [0,inf], [0,inf])\n        -1.74756459463318\n\n    Another example of a lattice sum in 2D::\n\n        >>> nsum(lambda x,y: (-1)**(x+y) / (x**2+y**2), [-inf,inf],\n        ...     [-inf,inf], ignore=True)\n        -2.1775860903036\n        >>> -pi*ln2\n        -2.1775860903036\n\n    An example of an Eisenstein series::\n\n        >>> nsum(lambda m,n: (m+n*1j)**(-4), [-inf,inf], [-inf,inf],\n        ...     ignore=True)\n        (3.1512120021539 + 0.0j)\n\n    **References**\n\n    1. [Weisstein]_ http://mathworld.wolfram.com/DoubleSeries.html,\n    2. [Weisstein]_ http://mathworld.wolfram.com/MadelungConstants.html\n\n    '
    (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 (+ g.adaptive_extrapolation(update, emfun, options))

dataset/ETHPy150Open fredrik-johansson/mpmath/mpmath/calculus/extrapolation.py nsum/VarMisuse@196/7983 `ctx`->`g`