@defun
def shanks(ctx, seq, table=None, randomized=False):
    r"""
    Given a list ``seq`` of the first `N` elements of a slowly
    convergent infinite sequence `(A_k)`, :func:`~mpmath.shanks` computes the iterated
    Shanks transformation `S(A), S(S(A)), \ldots, S^{N/2}(A)`. The Shanks
    transformation often provides strong convergence acceleration,
    especially if the sequence is oscillating.

    The iterated Shanks transformation is computed using the Wynn
    epsilon algorithm (see [1]). :func:`~mpmath.shanks` returns the full
    epsilon table generated by Wynn's algorithm, which can be read
    off as follows:

    * The table is a list of lists forming a lower triangular matrix,
      where higher row and column indices correspond to more accurate
      values.
    * The columns with even index hold dummy entries (required for the
      computation) and the columns with odd index hold the actual
      extrapolates.
    * The last element in the last row is typically the most
      accurate estimate of the limit.
    * The difference to the third last element in the last row
      provides an estimate of the approximation error.
    * The magnitude of the second last element provides an estimate
      of the numerical accuracy lost to cancellation.

    For convenience, so the extrapolation is stopped at an odd index
    so that ``shanks(seq)[-1][-1]`` always gives an estimate of the
    limit.

    Optionally, an existing table can be passed to :func:`~mpmath.shanks`.
    This can be used to efficiently extend a previous computation after
    new elements have been appended to the sequence. The table will
    then be updated in-place.

    **The Shanks transformation**

    The Shanks transformation is defined as follows (see [2]): given
    the input sequence `(A_0, A_1, \ldots)`, the transformed sequence is
    given by

    .. math ::

        S(A_k) = \frac{A_{k+1}A_{k-1}-A_k^2}{A_{k+1}+A_{k-1}-2 A_k}

    The Shanks transformation gives the exact limit `A_{\infty}` in a
    single step if `A_k = A + a q^k`. Note in particular that it
    extrapolates the exact sum of a geometric series in a single step.

    Applying the Shanks transformation once often improves convergence
    substantially for an arbitrary sequence, but the optimal effect is
    obtained by applying it iteratively:
    `S(S(A_k)), S(S(S(A_k))), \ldots`.

    Wynn's epsilon algorithm provides an efficient way to generate
    the table of iterated Shanks transformations. It reduces the
    computation of each element to essentially a single division, at
    the cost of requiring dummy elements in the table. See [1] for
    details.

    **Precision issues**

    Due to cancellation effects, the sequence must be typically be
    computed at a much higher precision than the target accuracy
    of the extrapolation.

    If the Shanks transformation converges to the exact limit (such
    as if the sequence is a geometric series), then a division by
    zero occurs. By default, :func:`~mpmath.shanks` handles this case by
    terminating the iteration and returning the table it has
    generated so far. With *randomized=True*, it will instead
    replace the zero by a pseudorandom number close to zero.
    (TODO: find a better solution to this problem.)

    **Examples**

    We illustrate by applying Shanks transformation to the Leibniz
    series for `\pi`::

        >>> from mpmath import *
        >>> mp.dps = 50
        >>> S = [4*sum(mpf(-1)**n/(2*n+1) for n in range(m))
        ...     for m in range(1,30)]
        >>>
        >>> T = shanks(S[:7])
        >>> for row in T:
        ...     nprint(row)
        ...
        [-0.75]
        [1.25, 3.16667]
        [-1.75, 3.13333, -28.75]
        [2.25, 3.14524, 82.25, 3.14234]
        [-2.75, 3.13968, -177.75, 3.14139, -969.937]
        [3.25, 3.14271, 327.25, 3.14166, 3515.06, 3.14161]

    The extrapolated accuracy is about 4 digits, and about 4 digits
    may have been lost due to cancellation::

        >>> L = T[-1]
        >>> nprint([abs(L[-1] - pi), abs(L[-1] - L[-3]), abs(L[-2])])
        [2.22532e-5, 4.78309e-5, 3515.06]

    Now we extend the computation::

        >>> T = shanks(S[:25], T)
        >>> L = T[-1]
        >>> nprint([abs(L[-1] - pi), abs(L[-1] - L[-3]), abs(L[-2])])
        [3.75527e-19, 1.48478e-19, 2.96014e+17]

    The value for pi is now accurate to 18 digits. About 18 digits may
    also have been lost to cancellation.

    Here is an example with a geometric series, where the convergence
    is immediate (the sum is exactly 1)::

        >>> mp.dps = 15
        >>> for row in shanks([0.5, 0.75, 0.875, 0.9375, 0.96875]):
        ...     nprint(row)
        [4.0]
        [8.0, 1.0]

    **References**

    1. [GravesMorris]_

    2. [BenderOrszag]_ pp. 368-375

    """
    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 = ctx.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/original