

@defun
def shanks(ctx, seq, table=None, randomized=False):
    "\n    Given a list ``seq`` of the first `N` elements of a slowly\n    convergent infinite sequence `(A_k)`, :func:`~mpmath.shanks` computes the iterated\n    Shanks transformation `S(A), S(S(A)), \\ldots, S^{N/2}(A)`. The Shanks\n    transformation often provides strong convergence acceleration,\n    especially if the sequence is oscillating.\n\n    The iterated Shanks transformation is computed using the Wynn\n    epsilon algorithm (see [1]). :func:`~mpmath.shanks` returns the full\n    epsilon table generated by Wynn's algorithm, which can be read\n    off as follows:\n\n    * The table is a list of lists forming a lower triangular matrix,\n      where higher row and column indices correspond to more accurate\n      values.\n    * The columns with even index hold dummy entries (required for the\n      computation) and the columns with odd index hold the actual\n      extrapolates.\n    * The last element in the last row is typically the most\n      accurate estimate of the limit.\n    * The difference to the third last element in the last row\n      provides an estimate of the approximation error.\n    * The magnitude of the second last element provides an estimate\n      of the numerical accuracy lost to cancellation.\n\n    For convenience, so the extrapolation is stopped at an odd index\n    so that ``shanks(seq)[-1][-1]`` always gives an estimate of the\n    limit.\n\n    Optionally, an existing table can be passed to :func:`~mpmath.shanks`.\n    This can be used to efficiently extend a previous computation after\n    new elements have been appended to the sequence. The table will\n    then be updated in-place.\n\n    **The Shanks transformation**\n\n    The Shanks transformation is defined as follows (see [2]): given\n    the input sequence `(A_0, A_1, \\ldots)`, the transformed sequence is\n    given by\n\n    .. math ::\n\n        S(A_k) = \\frac{A_{k+1}A_{k-1}-A_k^2}{A_{k+1}+A_{k-1}-2 A_k}\n\n    The Shanks transformation gives the exact limit `A_{\\infty}` in a\n    single step if `A_k = A + a q^k`. Note in particular that it\n    extrapolates the exact sum of a geometric series in a single step.\n\n    Applying the Shanks transformation once often improves convergence\n    substantially for an arbitrary sequence, but the optimal effect is\n    obtained by applying it iteratively:\n    `S(S(A_k)), S(S(S(A_k))), \\ldots`.\n\n    Wynn's epsilon algorithm provides an efficient way to generate\n    the table of iterated Shanks transformations. It reduces the\n    computation of each element to essentially a single division, at\n    the cost of requiring dummy elements in the table. See [1] for\n    details.\n\n    **Precision issues**\n\n    Due to cancellation effects, the sequence must be typically be\n    computed at a much higher precision than the target accuracy\n    of the extrapolation.\n\n    If the Shanks transformation converges to the exact limit (such\n    as if the sequence is a geometric series), then a division by\n    zero occurs. By default, :func:`~mpmath.shanks` handles this case by\n    terminating the iteration and returning the table it has\n    generated so far. With *randomized=True*, it will instead\n    replace the zero by a pseudorandom number close to zero.\n    (TODO: find a better solution to this problem.)\n\n    **Examples**\n\n    We illustrate by applying Shanks transformation to the Leibniz\n    series for `\\pi`::\n\n        >>> from mpmath import *\n        >>> mp.dps = 50\n        >>> S = [4*sum(mpf(-1)**n/(2*n+1) for n in range(m))\n        ...     for m in range(1,30)]\n        >>>\n        >>> T = shanks(S[:7])\n        >>> for row in T:\n        ...     nprint(row)\n        ...\n        [-0.75]\n        [1.25, 3.16667]\n        [-1.75, 3.13333, -28.75]\n        [2.25, 3.14524, 82.25, 3.14234]\n        [-2.75, 3.13968, -177.75, 3.14139, -969.937]\n        [3.25, 3.14271, 327.25, 3.14166, 3515.06, 3.14161]\n\n    The extrapolated accuracy is about 4 digits, and about 4 digits\n    may have been lost due to cancellation::\n\n        >>> L = T[-1]\n        >>> nprint([abs(L[-1] - pi), abs(L[-1] - L[-3]), abs(L[-2])])\n        [2.22532e-5, 4.78309e-5, 3515.06]\n\n    Now we extend the computation::\n\n        >>> T = shanks(S[:25], T)\n        >>> L = T[-1]\n        >>> nprint([abs(L[-1] - pi), abs(L[-1] - L[-3]), abs(L[-2])])\n        [3.75527e-19, 1.48478e-19, 2.96014e+17]\n\n    The value for pi is now accurate to 18 digits. About 18 digits may\n    also have been lost to cancellation.\n\n    Here is an example with a geometric series, where the convergence\n    is immediate (the sum is exactly 1)::\n\n        >>> mp.dps = 15\n        >>> for row in shanks([0.5, 0.75, 0.875, 0.9375, 0.96875]):\n        ...     nprint(row)\n        [4.0]\n        [8.0, 1.0]\n\n    **References**\n\n    1. [GravesMorris]_\n\n    2. [BenderOrszag]_ pp. 368-375\n\n    "
    if (len(seq) < 2):
        raise ValueError('seq should be of minimum length 2')
    if table:
        START = len(table)
    else:
        START = 0
        table = []
    STOP = (len(seq) - 1)
    if (STOP & 1):
        STOP -= 1
    one = row.one
    eps = (+ ctx.eps)
    if randomized:
        from random import Random
        rnd = Random()
        rnd.seed(START)
    for i in xrange(START, STOP):
        row = []
        for j in xrange((i + 1)):
            if (j == 0):
                (a, b) = (0, (seq[(i + 1)] - seq[i]))
            else:
                if (j == 1):
                    a = seq[i]
                else:
                    a = table[(i - 1)][(j - 2)]
                b = (row[(j - 1)] - table[(i - 1)][(j - 1)])
            if (not b):
                if randomized:
                    b = (rnd.getrandbits(10) * eps)
                elif (i & 1):
                    return table[:(- 1)]
                else:
                    return table
            row.append((a + (one / b)))
        table.append(row)
    return table

dataset/ETHPy150Open fredrik-johansson/mpmath/mpmath/calculus/extrapolation.py shanks/VarMisuse@98/2217 `ctx`->`row`