text: try:
      from itertools import izip
  except ImportError:
      izip = zip
  
  from ..libmp.backend import xrange
  from .calculus import defun
  
  try:
      next = next
  except NameError:
      next = lambda _: _.next()
  
  @defun
  def richardson(ctx, seq):
      r"""
      Given a list ``seq`` of the first `N` elements of a slowly convergent
      infinite sequence, :func:`~mpmath.richardson` computes the `N`-term
      Richardson extrapolate for the limit.
  
      :func:`~mpmath.richardson` returns `(v, c)` where `v` is the estimated
      limit and `c` is the magnitude of the largest weight used during the
      computation. The weight provides an estimate of the precision
      lost to cancellation. Due to cancellation effects, the sequence must
      be typically be computed at a much higher precision than the target
      accuracy of the extrapolation.
  
      **Applicability and issues**
  
      The `N`-step Richardson extrapolation algorithm used by
      :func:`~mpmath.richardson` is described in [1].
  
      Richardson extrapolation only works for a specific type of sequence,
      namely one converging like partial sums of
      `P(1)/Q(1) + P(2)/Q(2) + \ldots` where `P` and `Q` are polynomials.
      When the sequence does not convergence at such a rate
      :func:`~mpmath.richardson` generally produces garbage.
  
      Richardson extrapolation has the advantage of being fast: the `N`-term
      extrapolate requires only `O(N)` arithmetic operations, and usually
      produces an estimate that is accurate to `O(N)` digits. Contrast with
      the Shanks transformation (see :func:`~mpmath.shanks`), which requires
      `O(N^2)` operations.
  
      :func:`~mpmath.richardson` is unable to produce an estimate for the
      approximation error. One way to estimate the error is to perform
      two extrapolations with slightly different `N` and comparing the
      results.
  
      Richardson extrapolation does not work for oscillating sequences.
      As a simple workaround, :func:`~mpmath.richardson` detects if the last
      three elements do not differ monotonically, and in that case
      applies extrapolation only to the even-index elements.
  
      **Example**
  
      Applying Richardson extrapolation to the Leibniz series for `\pi`::
  
          >>> from mpmath import *
          >>> mp.dps = 30; mp.pretty = True
          >>> S = [4*sum(mpf(-1)**n/(2*n+1) for n in range(m))
          ...     for m in range(1,30)]
          >>> v, c = richardson(S[:10])
          >>> v
          3.2126984126984126984126984127
          >>> nprint([v-pi, c])
          [0.0711058, 2.0]
  
          >>> v, c = richardson(S[:30])
          >>> v
          3.14159265468624052829954206226
          >>> nprint([v-pi, c])
          [1.09645e-9, 20833.3]
  
      **References**
  
      1. [BenderOrszag]_ pp. 375-376
  
      """
      if len(seq) < 3:
          raise ValueError("seq should be of minimum length 3")
      if ctx.sign(seq[-1]-seq[-2]) != ctx.sign(seq[-2]-seq[-3]):
          seq = seq[::2]
      N = len(seq)//2-1
      s = ctx.zero
      # The general weight is c[k] = (N+k)**N * (-1)**(k+N) / k! / (N-k)!
      # To avoid repeated factorials, we simplify the quotient
      # of successive weights to obtain a recurrence relation
      c = (-1)**N * N**N / ctx.mpf(ctx._ifac(N))
      maxc = 1
      for k in xrange(N+1):
          s += c * seq[N+k]
          maxc = max(abs(c), maxc)
          c *= (k-N)*ctx.mpf(k+N+1)**N
          c /= ((1+k)*ctx.mpf(k+N)**N)
      return s, maxc
  
  @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
  
  
  class levin_class:
      # levin: Copyright 2013 Timo Hartmann (thartmann15 at gmail.com)
      r"""
      This interface implements Levin's (nonlinear) sequence transformation for
      convergence acceleration and summation of divergent series. It performs
      better than the Shanks/Wynn-epsilon algorithm for logarithmic convergent
      or alternating divergent series.
  
      Let *A* be the series we want to sum:
  
      .. math ::
  
          A = \sum_{k=0}^{\infty} a_k
  
      Attention: all `a_k` must be non-zero!
  
      Let `s_n` be the partial sums of this series:
  
      .. math ::
  
          s_n = \sum_{k=0}^n a_k.
  
      **Methods**
  
      Calling ``levin`` returns an object with the following methods.
  
      ``update(...)`` works with the list of individual terms `a_k` of *A*, and
      ``update_step(...)`` works with the list of partial sums `s_k` of *A*:
  
      .. code ::
  
          v, e = ...update([a_0, a_1,..., a_k])
          v, e = ...update_psum([s_0, s_1,..., s_k])
  
      ``step(...)`` works with the individual terms `a_k` and ``step_psum(...)``
      works with the partial sums `s_k`:
  
      .. code ::
  
          v, e = ...step(a_k)
          v, e = ...step_psum(s_k)
  
      *v* is the current estimate for *A*, and *e* is an error estimate which is
      simply the difference between the current estimate and the last estimate.
      One should not mix ``update``, ``update_psum``, ``step`` and ``step_psum``.
  
      **A word of caution**
  
      One can only hope for good results (i.e. convergence acceleration or
      resummation) if the `s_n` have some well defind asymptotic behavior for
      large `n` and are not erratic or random. Furthermore one usually needs very
      high working precision because of the numerical cancellation. If the working
      precision is insufficient, levin may produce silently numerical garbage.
      Furthermore even if the Levin-transformation converges, in the general case
      there is no proof that the result is mathematically sound. Only for very
      special classes of problems one can prove that the Levin-transformation
      converges to the expected result (for example Stieltjes-type integrals).
      Furthermore the Levin-transform is quite expensive (i.e. slow) in comparison
      to Shanks/Wynn-epsilon, Richardson & co.
      In summary one can say that the Levin-transformation is powerful but
      unreliable and that it may need a copious amount of working precision.
  
      The Levin transform has several variants differing in the choice of weights.
      Some variants are better suited for the possible flavours of convergence
      behaviour of *A* than other variants:
  
      .. code ::
  
         convergence behaviour   levin-u   levin-t   levin-v   shanks/wynn-epsilon
  
         logarithmic               +         -         +           -
         linear                    +         +         +           +
         alternating divergent     +         +         +           +
  
           "+" means the variant is suitable,"-" means the variant is not suitable;
           for comparison the Shanks/Wynn-epsilon transform is listed, too.
  
      The variant is controlled though the variant keyword (i.e. ``variant="u"``,
      ``variant="t"`` or ``variant="v"``). Overall "u" is probably the best choice.
  
      Finally it is possible to use the Sidi-S transform instead of the Levin transform
      by using the keyword ``method='sidi'``. The Sidi-S transform works better than the
      Levin transformation for some divergent series (see the examples).
  
      Parameters:
  
      .. code ::
  
         method      "levin" or "sidi" chooses either the Levin or the Sidi-S transformation
         variant     "u","t" or "v" chooses the weight variant.
  
      The Levin transform is also accessible through the nsum interface.
      ``method="l"`` or ``method="levin"`` select the normal Levin transform while
      ``method="sidi"``
      selects the Sidi-S transform. The variant is in both cases selected through the
      levin_variant keyword. The stepsize in :func:`~mpmath.nsum` must not be chosen too large, otherwise
      it will miss the point where the Levin transform converges resulting in numerical
      overflow/garbage. For highly divergent series a copious amount of working precision
      must be chosen.
  
      **Examples**
  
      First we sum the zeta function::
  
          >>> from mpmath import mp
          >>> mp.prec = 53
          >>> eps = mp.mpf(mp.eps)
          >>> with mp.extraprec(2 * mp.prec): # levin needs a high working precision
          ...     L = mp.levin(method = "levin", variant = "u")
          ...     S, s, n = [], 0, 1
          ...     while 1:
          ...         s += mp.one / (n * n)
          ...         n += 1
          ...         S.append(s)
          ...         v, e = L.update_psum(S)
          ...         if e < eps:
          ...             break
          ...         if n > 1000: raise RuntimeError("iteration limit exceeded")
          >>> print(mp.chop(v - mp.pi ** 2 / 6))
          0.0
          >>> w = mp.nsum(lambda n: 1 / (n*n), [1, mp.inf], method = "levin", levin_variant = "u")
          >>> print(mp.chop(v - w))
          0.0
  
      Now we sum the zeta function outside its range of convergence
      (attention: This does not work at the negative integers!)::
  
          >>> eps = mp.mpf(mp.eps)
          >>> with mp.extraprec(2 * mp.prec): # levin needs a high working precision
          ...     L = mp.levin(method = "levin", variant = "v")
          ...     A, n = [], 1
          ...     while 1:
          ...         s = mp.mpf(n) ** (2 + 3j)
          ...         n += 1
          ...         A.append(s)
          ...         v, e = L.update(A)
          ...         if e < eps:
          ...             break
          ...         if n > 1000: raise RuntimeError("iteration limit exceeded")
          >>> print(mp.chop(v - mp.zeta(-2-3j)))
          0.0
          >>> w = mp.nsum(lambda n: n ** (2 + 3j), [1, mp.inf], method = "levin", levin_variant = "v")
          >>> print(mp.chop(v - w))
          0.0
  
      Now we sum the divergent asymptotic expansion of an integral related to the
      exponential integral (see also [2] p.373). The Sidi-S transform works best here::
  
          >>> 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
          >>> eps = mp.mpf(mp.eps)
          >>> with mp.extraprec(2 * mp.prec): # high working precisions are mandatory for divergent resummation
          ...     L = mp.levin(method = "sidi", variant = "t")
          ...     n = 0
          ...     while 1:
          ...         s = (-1)**n * mp.fac(n) * z ** (-n)
          ...         v, e = L.step(s)
          ...         n += 1
          ...         if e < eps:
          ...             break
          ...         if n > 1000: raise RuntimeError("iteration limit exceeded")
          >>> print(mp.chop(v - exact))
          0.0
          >>> w = mp.nsum(lambda n: (-1) ** n * mp.fac(n) * z ** (-n), [0, mp.inf], method = "sidi", levin_variant = "t")
          >>> print(mp.chop(v - w))
          0.0
  
      Another highly divergent integral is also summable::
  
          >>> z = mp.mpf(2)
          >>> eps = mp.mpf(mp.eps)
          >>> 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
          >>> with mp.extraprec(7 * mp.prec):  # we need copious amount of precision to sum this highly divergent series
          ...     L = mp.levin(method = "levin", variant = "t")
          ...     n, s = 0, 0
          ...     while 1:
          ...         s += (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n))
          ...         n += 1
          ...         v, e = L.step_psum(s)
          ...         if e < eps:
          ...             break
          ...         if n > 1000: raise RuntimeError("iteration limit exceeded")
          >>> print(mp.chop(v - exact))
          0.0
          >>> 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(v - w))
          0.0
  
      These examples run with 15-20 decimal digits precision. For higher precision the
      working precision must be raised.
  
      **Examples for nsum**
  
      Here we calculate Euler's constant as the constant term in the Laurent
      expansion of `\zeta(s)` at `s=1`. This sum converges extremly slowly 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 = "l") - 1 / z
          >>> print(mp.chop(a - mp.euler, tol = 1e-10))
          0.0
  
      The Sidi-S transform performs excellently for the alternating series of `\log(2)`::
  
          >>> a = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "sidi")
          >>> print(mp.chop(a - mp.log(2)))
          0.0
  
      Hypergeometric series can also be summed outside their range of convergence.
      The stepsize in :func:`~mpmath.nsum` must not be chosen too large, otherwise it will miss the
      point where the Levin transform converges resulting in numerical overflow/garbage::
  
          >>> 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
  
      References:
  
        [1] E.J. Weniger - "Nonlinear Sequence Transformations for the Acceleration of
            Convergence and the Summation of Divergent Series" arXiv:math/0306302
  
        [2] A. Sidi - "Pratical Extrapolation Methods"
  
        [3] H.H.H. Homeier - "Scalar Levin-Type Sequence Transformations" arXiv:math/0005209
  
      """
  
      def __init__(self, method = "levin", variant = "u"):
          self.variant = variant
          self.n = 0
          self.a0 = 0
          self.theta = 1
          self.A = []
          self.B = []
          self.last = 0
          self.last_s = False
  
          if method == "levin":
              self.factor = self.factor_levin
          elif method == "sidi":
              self.factor = self.factor_sidi
          else:
              raise ValueError("levin: unknown method \"%s\"" % method)
  
      def factor_levin(self, i):
          # original levin
          # [1] p.50,e.7.5-7 (with n-j replaced by i)
          return (self.theta + i) * (self.theta + self.n - 1) ** (self.n - i - 2) / self.ctx.mpf(self.theta + self.n) ** (self.n - i - 1)
  
      def factor_sidi(self, i):
          # sidi analogon to levin (factorial series)
          # [1] p.59,e.8.3-16 (with n-j replaced by i)
          return (self.theta + self.n - 1) * (self.theta + self.n - 2) / self.ctx.mpf((self.theta + 2 * self.n - i - 2) * (self.theta + 2 * self.n - i - 3))
  
      def run(self, s, a0, a1 = 0):
          if self.variant=="t":
              # levin t
              w=a0
          elif self.variant=="u":
              # levin u
              w=a0*(self.theta+self.n)
          elif self.variant=="v":
              # levin v
              w=a0*a1/(a0-a1)
          else:
              assert False, "unknown variant"
  
          if w==0:
              raise ValueError("levin: zero weight")
  
          self.A.append(s/w)
          self.B.append(1/w)
  
          for i in range(self.n-1,-1,-1):
              if i==self.n-1:
                  f=1
              else:
                  f=self.factor(i)
  
              self.A[i]=self.A[i+1]-f*self.A[i]
              self.B[i]=self.B[i+1]-f*self.B[i]
  
          self.n+=1
  
      ###########################################################################
  
      def update_psum(self,S):
          """
          This routine applies the convergence acceleration to the list of partial sums.
  
          A   = sum(a_k, k = 0..infinity)
          s_n = sum(a_k, k = 0..n)
  
          v, e = ...update_psum([s_0, s_1,..., s_k])
  
          output:
            v      current estimate of the series A
            e      an error estimate which is simply the difference between the current
                   estimate and the last estimate.
          """
  
          if self.variant!="v":
              if self.n==0:
                  self.run(S[0],S[0])
              while self.n<len(S):
                  self.run(S[self.n],S[self.n]-S[self.n-1])
          else:
              if len(S)==1:
                  self.last=0
                  return S[0],abs(S[0])
  
              if self.n==0:
                  self.a1=S[1]-S[0]
                  self.run(S[0],S[0],self.a1)
  
              while self.n<len(S)-1:
                  na1=S[self.n+1]-S[self.n]
                  self.run(S[self.n],self.a1,na1)
                  self.a1=na1
  
          value=self.A[0]/self.B[0]
          err=abs(value-self.last)
          self.last=value
  
          return value,err
  
      def update(self,X):
          """
          This routine applies the convergence acceleration to the list of individual terms.
  
          A = sum(a_k, k = 0..infinity)
  
          v, e = ...update([a_0, a_1,..., a_k])
  
          output:
            v      current estimate of the series A
            e      an error estimate which is simply the difference between the current
                   estimate and the last estimate.
          """
  
          if self.variant!="v":
              if self.n==0:
                  self.s=X[0]
                  self.run(self.s,X[0])
              while self.n<len(X):
                  self.s+=X[self.n]
                  self.run(self.s,X[self.n])
          else:
              if len(X)==1:
                  self.last=0
                  return X[0],abs(X[0])
  
              if self.n==0:
                  self.s=X[0]
                  self.run(self.s,X[0],X[1])
  
              while self.n<len(X)-1:
                  self.s+=X[self.n]
                  self.run(self.s,X[self.n],X[self.n+1])
  
          value=self.A[0]/self.B[0]
          err=abs(value-self.last)
          self.last=value
  
          return value,err
  
      ###########################################################################
  
      def step_psum(self,s):
          """
          This routine applies the convergence acceleration to the partial sums.
  
          A   = sum(a_k, k = 0..infinity)
          s_n = sum(a_k, k = 0..n)
  
          v, e = ...step_psum(s_k)
  
          output:
            v      current estimate of the series A
            e      an error estimate which is simply the difference between the current
                   estimate and the last estimate.
          """
  
          if self.variant!="v":
              if self.n==0:
                  self.last_s=s
                  self.run(s,s)
              else:
                  self.run(s,s-self.last_s)
                  self.last_s=s
          else:
              if isinstance(self.last_s,bool):
                  self.last_s=s
                  self.last_w=s
                  self.last=0
                  return s,abs(s)
  
              na1=s-self.last_s
              self.run(self.last_s,self.last_w,na1)
              self.last_w=na1
              self.last_s=s
  
          value=self.A[0]/self.B[0]
          err=abs(value-self.last)
          self.last=value
  
          return value,err
  
      def step(self,x):
          """
          This routine applies the convergence acceleration to the individual terms.
  
          A = sum(a_k, k = 0..infinity)
  
          v, e = ...step(a_k)
  
          output:
            v      current estimate of the series A
            e      an error estimate which is simply the difference between the current
                   estimate and the last estimate.
          """
  
          if self.variant!="v":
              if self.n==0:
                  self.s=x
                  self.run(self.s,x)
              else:
                  self.s+=x
                  self.run(self.s,x)
          else:
              if isinstance(self.last_s,bool):
                  self.last_s=x
                  self.s=0
                  self.last=0
                  return x,abs(x)
  
              self.s+=self.last_s
              self.run(self.s,self.last_s,x)
              self.last_s=x
  
          value=self.A[0]/self.B[0]
          err=abs(value-self.last)
          self.last=value
  
          return value,err
  
  def levin(ctx, method = "levin", variant = "u"):
      L = levin_class(method = method, variant = variant)
      L.ctx = ctx
      return L
  
  levin.__doc__ = levin_class.__doc__
  defun(levin)
  
  
  class cohen_alt_class:
      # cohen_alt: Copyright 2013 Timo Hartmann (thartmann15 at gmail.com)
      """
      This interface implements the convergence acceleration of alternating series
      as described in H. Cohen, F.R. Villegas, D. Zagier - "Convergence Acceleration
      of Alternating Series". This series transformation works only well if the
      individual terms of the series have an alternating sign. It belongs to the
      class of linear series transformations (in contrast to the Shanks/Wynn-epsilon
      or Levin transform). This series transformation is also able to sum some types
      of divergent series. See the paper under which conditions this resummation is
      mathematical sound.
  
      Let *A* be the series we want to sum:
  
      .. math ::
  
          A = \sum_{k=0}^{\infty} a_k
  
      Let `s_n` be the partial sums of this series:
  
      .. math ::
  
          s_n = \sum_{k=0}^n a_k.
  
  
      **Interface**
  
      Calling ``cohen_alt`` returns an object with the following methods.
  
      Then ``update(...)`` works with the list of individual terms `a_k` and
      ``update_psum(...)`` works with the list of partial sums `s_k`:
  
      .. code ::
  
          v, e = ...update([a_0, a_1,..., a_k])
          v, e = ...update_psum([s_0, s_1,..., s_k])
  
      *v* is the current estimate for *A*, and *e* is an error estimate which is
      simply the difference between the current estimate and the last estimate.
  
      **Examples**
  
      Here we compute the alternating zeta function using ``update_psum``::
  
          >>> from mpmath import mp
          >>> AC = mp.cohen_alt()
          >>> S, s, n = [], 0, 1
          >>> while 1:
          ...     s += -((-1) ** n) * mp.one / (n * n)
          ...     n += 1
          ...     S.append(s)
          ...     v, e = AC.update_psum(S)
          ...     if e < mp.eps:
          ...         break
          ...     if n > 1000: raise RuntimeError("iteration limit exceeded")
          >>> print(mp.chop(v - mp.pi ** 2 / 12))
          0.0
  
      Here we compute the product `\prod_{n=1}^{\infty} \Gamma(1+1/(2n-1)) / \Gamma(1+1/(2n))`::
  
          >>> A = []
          >>> AC = mp.cohen_alt()
          >>> n = 1
          >>> while 1:
          ...     A.append( mp.loggamma(1 + mp.one / (2 * n - 1)))
          ...     A.append(-mp.loggamma(1 + mp.one / (2 * n)))
          ...     n += 1
          ...     v, e = AC.update(A)
          ...     if e < mp.eps:
          ...         break
          ...     if n > 1000: raise RuntimeError("iteration limit exceeded")
          >>> v = mp.exp(v)
          >>> print(mp.chop(v - 1.06215090557106, tol = 1e-12))
          0.0
  
      ``cohen_alt`` is also accessible through the :func:`~mpmath.nsum` interface::
  
          >>> v = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "a")
          >>> print(mp.chop(v - mp.log(2)))
          0.0
          >>> v = mp.nsum(lambda n: (-1)**n / (2 * n + 1), [0, mp.inf], method = "a")
          >>> print(mp.chop(v - mp.pi / 4))
          0.0
          >>> 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
  
      """
  
      def __init__(self):
          self.last=0
  
      def update(self, A):
          """
          This routine applies the convergence acceleration to the list of individual terms.
  
          A    = sum(a_k, k = 0..infinity)
  
          v, e = ...update([a_0, a_1,..., a_k])
  
          output:
            v      current estimate of the series A
            e      an error estimate which is simply the difference between the current
                   estimate and the last estimate.
          """
  
          n = len(A)
          d = (3 + self.ctx.sqrt(8)) ** n
          d = (d + 1 / d) / 2
          b = -self.ctx.one
          c = -d
          s = 0
  
          for k in xrange(n):
              c = b - c
              if k % 2 == 0:
                  s = s + c * A[k]
              else:
                  s = s - c * A[k]
              b = 2 * (k + n) * (k - n) * b / ((2 * k + 1) * (k + self.ctx.one))
  
          value = s / d
  
          err = abs(value - self.last)
          self.last = value
  
          return value, err
  
      def update_psum(self, S):
          """
          This routine applies the convergence acceleration to the list of partial sums.
  
          A   = sum(a_k, k = 0..infinity)
          s_n = sum(a_k ,k = 0..n)
  
          v, e = ...update_psum([s_0, s_1,..., s_k])
  
          output:
            v      current estimate of the series A
            e      an error estimate which is simply the difference between the current
                   estimate and the last estimate.
          """
  
          n = len(S)
          d = (3 + self.ctx.sqrt(8)) ** n
          d = (d + 1 / d) / 2
          b = self.ctx.one
          s = 0
  
          for k in xrange(n):
              b = 2 * (n + k) * (n - k) * b / ((2 * k + 1) * (k + self.ctx.one))
              s += b * S[k]
  
          value = s / d
  
          err = abs(value - self.last)
          self.last = value
  
          return value, err
  
  def cohen_alt(ctx):
      L = cohen_alt_class()
      L.ctx = ctx
      return L
  
  cohen_alt.__doc__ = cohen_alt_class.__doc__
  defun(cohen_alt)
  
  
  @defun
  def sumap(ctx, f, interval, integral=None, error=False):
      r"""
      Evaluates an infinite series of an analytic summand *f* using the
      Abel-Plana formula
  
      .. math ::
  
          \sum_{k=0}^{\infty} f(k) = \int_0^{\infty} f(t) dt + \frac{1}{2} f(0) +
              i \int_0^{\infty} \frac{f(it)-f(-it)}{e^{2\pi t}-1} dt.
  
      Unlike the Euler-Maclaurin formula (see :func:`~mpmath.sumem`),
      the Abel-Plana formula does not require derivatives. However,
      it only works when `|f(it)-f(-it)|` does not
      increase too rapidly with `t`.
  
      **Examples**
  
      The Abel-Plana formula is particularly useful when the summand
      decreases like a power of `k`; for example when the sum is a pure
      zeta function::
  
          >>> from mpmath import *
          >>> mp.dps = 25; mp.pretty = True
          >>> sumap(lambda k: 1/k**2.5, [1,inf])
          1.34148725725091717975677
          >>> zeta(2.5)
          1.34148725725091717975677
          >>> sumap(lambda k: 1/(k+1j)**(2.5+2.5j), [1,inf])
          (-3.385361068546473342286084 - 0.7432082105196321803869551j)
          >>> zeta(2.5+2.5j, 1+1j)
          (-3.385361068546473342286084 - 0.7432082105196321803869551j)
  
      If the series is alternating, numerical quadrature along the real
      line is likely to give poor results, so it is better to evaluate
      the first term symbolically whenever possible:
  
          >>> n=3; z=-0.75
          >>> I = expint(n,-log(z))
          >>> chop(sumap(lambda k: z**k / k**n, [1,inf], integral=I))
          -0.6917036036904594510141448
          >>> polylog(n,z)
          -0.6917036036904594510141448
  
      """
      prec = ctx.prec
      try:
          ctx.prec += 10
          a, b = interval
          if  b != ctx.inf:
              raise ValueError("b should be equal to ctx.inf")
          g = lambda x: f(x+a)
          if integral is None:
              i1, err1 = ctx.quad(g, [0,ctx.inf], error=True)
          else:
              i1, err1 = integral, 0
          j = ctx.j
          p = ctx.pi * 2
          if ctx._is_real_type(i1):
              h = lambda t: -2 * ctx.im(g(j*t)) / ctx.expm1(p*t)
          else:
              h = lambda t: j*(g(j*t)-g(-j*t)) / ctx.expm1(p*t)
          i2, err2 = ctx.quad(h, [0,ctx.inf], error=True)
          err = err1+err2
          v = i1+i2+0.5*g(ctx.mpf(0))
      finally:
          ctx.prec = prec
      if error:
          return +v, err
      return +v
  
  
  @defun
  def sumem(ctx, f, interval, tol=None, reject=10, integral=None,
      adiffs=None, bdiffs=None, verbose=False, error=False,
      _fast_abort=False):
      r"""
      Uses the Euler-Maclaurin formula to compute an approximation accurate
      to within ``tol`` (which defaults to the present epsilon) of the sum
  
      .. math ::
  
          S = \sum_{k=a}^b f(k)
  
      where `(a,b)` are given by ``interval`` and `a` or `b` may be
      infinite. The approximation is
  
      .. math ::
  
          S \sim \int_a^b f(x) \,dx + \frac{f(a)+f(b)}{2} +
          \sum_{k=1}^{\infty} \frac{B_{2k}}{(2k)!}
          \left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right).
  
      The last sum in the Euler-Maclaurin formula is not generally
      convergent (a notable exception is if `f` is a polynomial, in
      which case Euler-Maclaurin actually gives an exact result).
  
      The summation is stopped as soon as the quotient between two
      consecutive terms falls below *reject*. That is, by default
      (*reject* = 10), the summation is continued as long as each
      term adds at least one decimal.
  
      Although not convergent, convergence to a given tolerance can
      often be "forced" if `b = \infty` by summing up to `a+N` and then
      applying the Euler-Maclaurin formula to the sum over the range
      `(a+N+1, \ldots, \infty)`. This procedure is implemented by
      :func:`~mpmath.nsum`.
  
      By default numerical quadrature and differentiation is used.
      If the symbolic values of the integral and endpoint derivatives
      are known, it is more efficient to pass the value of the
      integral explicitly as ``integral`` and the derivatives
      explicitly as ``adiffs`` and ``bdiffs``. The derivatives
      should be given as iterables that yield
      `f(a), f'(a), f''(a), \ldots` (and the equivalent for `b`).
  
      **Examples**
  
      Summation of an infinite series, with automatic and symbolic
      integral and derivative values (the second should be much faster)::
  
          >>> from mpmath import *
          >>> mp.dps = 50; mp.pretty = True
          >>> sumem(lambda n: 1/n**2, [32, inf])
          0.03174336652030209012658168043874142714132886413417
          >>> I = mpf(1)/32
          >>> D = adiffs=((-1)**n*fac(n+1)*32**(-2-n) for n in range(999))
          >>> sumem(lambda n: 1/n**2, [32, inf], integral=I, adiffs=D)
          0.03174336652030209012658168043874142714132886413417
  
      An exact evaluation of a finite polynomial sum::
  
          >>> sumem(lambda n: n**5-12*n**2+3*n, [-100000, 200000])
          10500155000624963999742499550000.0
          >>> print(sum(n**5-12*n**2+3*n for n in range(-100000, 200001)))
          10500155000624963999742499550000
  
      """
      tol = tol or +ctx.eps
      interval = ctx._as_points(interval)
      a = ctx.convert(interval[0])
      b = ctx.convert(interval[-1])
      err = ctx.zero
      prev = 0
      M = 10000
      if a == ctx.ninf: adiffs = (0 for n in xrange(M))
      else:             adiffs = adiffs or ctx.diffs(f, a)
      if b == ctx.inf:  bdiffs = (0 for n in xrange(M))
      else:             bdiffs = bdiffs or ctx.diffs(f, b)
      orig = ctx.prec
      #verbose = 1
      try:
          ctx.prec += 10
          s = ctx.zero
          for k, (da, db) in enumerate(izip(adiffs, bdiffs)):
              if k & 1:
                  term = (db-da) * ctx.bernoulli(k+1) / ctx.factorial(k+1)
                  mag = abs(term)
                  if verbose:
                      print("term", k, "magnitude =", ctx.nstr(mag))
                  if k > 4 and mag < tol:
                      s += term
                      break
                  elif k > 4 and abs(prev) / mag < reject:
                      err += mag
                      if _fast_abort:
                          return [s, (s, err)][error]
                      if verbose:
                          print("Failed to converge")
                      break
                  else:
                      s += term
                  prev = term
          # Endpoint correction
          if a != ctx.ninf: s += f(a)/2
          if b != ctx.inf: s += f(b)/2
          # Tail integral
          if verbose:
              print("Integrating f(x) from x = %s to %s" % (ctx.nstr(a), ctx.nstr(b)))
          if integral:
              s += integral
          else:
              integral, ierr = ctx.quad(f, interval, error=True)
              if verbose:
                  print("Integration error:", ierr)
              s += integral
              err += ierr
      finally:
          ctx.prec = orig
      if error:
          return s, err
      else:
          return s
  
  @defun
  def adaptive_extrapolation(ctx, update, emfun, kwargs):
      option = kwargs.get
      if ctx._fixed_precision:
          tol = option('tol', ctx.eps*2**10)
      else:
          tol = option('tol', ctx.eps/2**10)
      verbose = option('verbose', False)
      maxterms = option('maxterms', ctx.dps*10)
      method = set(option('method', 'r+s').split('+'))
      skip = option('skip', 0)
      steps = iter(option('steps', xrange(10, 10**9, 10)))
      strict = option('strict')
      #steps = (10 for i in xrange(1000))
      summer=[]
      if 'd' in method or 'direct' in method:
          TRY_RICHARDSON = TRY_SHANKS = TRY_EULER_MACLAURIN = False
      else:
          TRY_RICHARDSON = ('r' in method) or ('richardson' in method)
          TRY_SHANKS = ('s' in method) or ('shanks' in method)
          TRY_EULER_MACLAURIN = ('e' in method) or \
              ('euler-maclaurin' in method)
  
          def init_levin(m):
              variant = kwargs.get("levin_variant", "u")
              if isinstance(variant, str):
                  if variant == "all":
                      variant = ["u", "v", "t"]
                  else:
                      variant = [variant]
              for s in variant:
                  L = levin_class(method = m, variant = s)
                  L.ctx = ctx
                  L.name = m + "(" + s + ")"
                  summer.append(L)
  
          if ('l' in method) or ('levin' in method):
              init_levin("levin")
  
          if ('sidi' in method):
              init_levin("sidi")
  
          if ('a' in method) or ('alternating' in method):
              L = cohen_alt_class()
              L.ctx = ctx
              L.name = "alternating"
              summer.append(L)
  
      last_richardson_value = 0
      shanks_table = []
      index = 0
      step = 10
      partial = []
      best = ctx.zero
      orig = ctx.prec
      try:
          if 'workprec' in kwargs:
              ctx.prec = kwargs['workprec']
          elif TRY_RICHARDSON or TRY_SHANKS or len(summer)!=0:
              ctx.prec = (ctx.prec+10) * 4
          else:
              ctx.prec += 30
          while 1:
              if index >= maxterms:
                  break
  
              # Get new batch of terms
              try:
                  step = next(steps)
              except StopIteration:
                  pass
              if verbose:
                  print("-"*70)
                  print("Adding terms #%i-#%i" % (index, index+step))
              update(partial, xrange(index, index+step))
              index += step
  
              # Check direct error
              best = partial[-1]
              error = abs(best - partial[-2])
              if verbose:
                  print("Direct error: %s" % ctx.nstr(error))
              if error <= tol:
                  return best
  
              # Check each extrapolation method
              if TRY_RICHARDSON:
                  value, maxc = ctx.richardson(partial)
                  # Convergence
                  richardson_error = abs(value - last_richardson_value)
                  if verbose:
                      print("Richardson error: %s" % ctx.nstr(richardson_error))
                  # Convergence
                  if richardson_error <= tol:
                      return value
                  last_richardson_value = value
                  # Unreliable due to cancellation
                  if ctx.eps*maxc > tol:
                      if verbose:
                          print("Ran out of precision for Richardson")
                      TRY_RICHARDSON = False
                  if richardson_error < error:
                      error = richardson_error
                      best = value
              if TRY_SHANKS:
                  shanks_table = ctx.shanks(partial, shanks_table, randomized=True)
                  row = shanks_table[-1]
                  if len(row) == 2:
                      est1 = row[-1]
                      shanks_error = 0
                  else:
                      est1, maxc, est2 = row[-1], abs(row[-2]), row[-3]
                      shanks_error = abs(est1-est2)
                  if verbose:
                      print("Shanks error: %s" % ctx.nstr(shanks_error))
                  if shanks_error <= tol:
                      return est1
                  if ctx.eps*maxc > tol:
                      if verbose:
                          print("Ran out of precision for Shanks")
                      TRY_SHANKS = False
                  if shanks_error < error:
                      error = shanks_error
                      best = est1
              for L in summer:
                  est, lerror = L.update_psum(partial)
                  if verbose:
                      print("%s error: %s" % (L.name, ctx.nstr(lerror)))
                  if lerror <= tol:
                      return est
                  if lerror < error:
                      error = lerror
                      best = est
              if TRY_EULER_MACLAURIN:
                  if ctx.mpc(ctx.sign(partial[-1]) / ctx.sign(partial[-2])).ae(-1):
                      if verbose:
                          print ("NOT using Euler-Maclaurin: the series appears"
                              " to be alternating, so numerical\n quadrature"
                              " will most likely fail")
                      TRY_EULER_MACLAURIN = False
                  else:
                      value, em_error = emfun(index, tol)
                      value += partial[-1]
                      if verbose:
                          print("Euler-Maclaurin error: %s" % ctx.nstr(em_error))
                      if em_error <= tol:
                          return value
                      if em_error < error:
                          best = value
      finally:
          ctx.prec = orig
      if strict:
          raise ctx.NoConvergence
      if verbose:
          print("Warning: failed to converge to target accuracy")
      return best
  
  @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)
  
  
  def wrapsafe(f):
      def g(*args):
          try:
              return f(*args)
          except (ArithmeticError, ValueError):
              return 0
      return g
  
  def standardize(ctx, f, intervals, options):
      if options.get("ignore"):
          f = wrapsafe(f)
      finite = []
      infinite = []
      for k, points in enumerate(intervals):
          a, b = ctx._as_points(points)
          if b < a:
              return False, (lambda: ctx.zero)
          if a == ctx.ninf or b == ctx.inf:
              infinite.append((k, (a,b)))
          else:
              finite.append((k, (int(a), int(b))))
      if finite:
          f = fold_finite(ctx, f, finite)
          if not infinite:
              return False, lambda: f(*([0]*len(intervals)))
      if infinite:
          f = standardize_infinite(ctx, f, infinite)
          f = fold_infinite(ctx, f, infinite)
          args = [0] * len(intervals)
          d = infinite[0][0]
          def g(k):
              args[d] = k
              return f(*args)
          return True, g
  
  # backwards compatible itertools.product
  def cartesian_product(args):
      pools = map(tuple, args)
      result = [[]]
      for pool in pools:
          result = [x+[y] for x in result for y in pool]
      for prod in result:
          yield tuple(prod)
  
  def fold_finite(ctx, f, intervals):
      if not intervals:
          return f
      indices = [v[0] for v in intervals]
      points = [v[1] for v in intervals]
      ranges = [xrange(a, b+1) for (a,b) in points]
      def g(*args):
          args = list(args)
          s = ctx.zero
          for xs in cartesian_product(ranges):
              for dim, x in zip(indices, xs):
                  args[dim] = ctx.mpf(x)
              s += f(*args)
          return s
      #print "Folded finite", indices
      return g
  
  # Standardize each interval to [0,inf]
  def standardize_infinite(ctx, f, intervals):
      if not intervals:
          return f
      dim, [a,b] = intervals[-1]
      if a == ctx.ninf:
          if b == ctx.inf:
              def g(*args):
                  args = list(args)
                  k = args[dim]
                  if k:
                      s = f(*args)
                      args[dim] = -k
                      s += f(*args)
                      return s
                  else:
                      return f(*args)
          else:
              def g(*args):
                  args = list(args)
                  args[dim] = b - args[dim]
                  return f(*args)
      else:
          def g(*args):
              args = list(args)
              args[dim] += a
              return f(*args)
      #print "Standardized infinity along dimension", dim, a, b
      return standardize_infinite(ctx, g, intervals[:-1])
  
  def fold_infinite(ctx, f, intervals):
      if len(intervals) < 2:
          return f
      dim1 = intervals[-2][0]
      dim2 = intervals[-1][0]
      # Assume intervals are [0,inf] x [0,inf] x ...
      def g(*args):
          args = list(args)
          #args.insert(dim2, None)
          n = int(args[dim1])
          s = ctx.zero
          #y = ctx.mpf(n)
          args[dim2] = ctx.mpf(n) #y
          for x in xrange(n+1):
              args[dim1] = ctx.mpf(x)
              s += f(*args)
          args[dim1] = ctx.mpf(n) #ctx.mpf(n)
          for y in xrange(n):
              args[dim2] = ctx.mpf(y)
              s += f(*args)
          return s
      #print "Folded infinite from", len(intervals), "to", (len(intervals)-1)
      return fold_infinite(ctx, g, intervals[:-1])
  
  @defun
  def nprod(ctx, f, interval, nsum=False, **kwargs):
      r"""
      Computes the product
  
      .. math ::
  
          P = \prod_{k=a}^b f(k)
  
      where `(a, b)` = *interval*, and where `a = -\infty` and/or
      `b = \infty` are allowed.
  
      By default, :func:`~mpmath.nprod` uses the same extrapolation methods as
      :func:`~mpmath.nsum`, except applied to the partial products rather than
      partial sums, and the same keyword options as for :func:`~mpmath.nsum` are
      supported. If ``nsum=True``, the product is instead computed via
      :func:`~mpmath.nsum` as
  
      .. math ::
  
          P = \exp\left( \sum_{k=a}^b \log(f(k)) \right).
  
      This is slower, but can sometimes yield better results. It is
      also required (and used automatically) when Euler-Maclaurin
      summation is requested.
  
      **Examples**
  
      A simple finite product::
  
          >>> from mpmath import *
          >>> mp.dps = 25; mp.pretty = True
          >>> nprod(lambda k: k, [1, 4])
          24.0
  
      A large number of infinite products have known exact values,
      and can therefore be used as a reference. Most of the following
      examples are taken from MathWorld [1].
  
      A few infinite products with simple values are::
  
          >>> 2*nprod(lambda k: (4*k**2)/(4*k**2-1), [1, inf])
          3.141592653589793238462643
          >>> nprod(lambda k: (1+1/k)**2/(1+2/k), [1, inf])
          2.0
          >>> nprod(lambda k: (k**3-1)/(k**3+1), [2, inf])
          0.6666666666666666666666667
          >>> nprod(lambda k: (1-1/k**2), [2, inf])
          0.5
  
      Next, several more infinite products with more complicated
      values::
  
          >>> nprod(lambda k: exp(1/k**2), [1, inf]); exp(pi**2/6)
          5.180668317897115748416626
          5.180668317897115748416626
  
          >>> nprod(lambda k: (k**2-1)/(k**2+1), [2, inf]); pi*csch(pi)
          0.2720290549821331629502366
          0.2720290549821331629502366
  
          >>> nprod(lambda k: (k**4-1)/(k**4+1), [2, inf])
          0.8480540493529003921296502
          >>> pi*sinh(pi)/(cosh(sqrt(2)*pi)-cos(sqrt(2)*pi))
          0.8480540493529003921296502
  
          >>> nprod(lambda k: (1+1/k+1/k**2)**2/(1+2/k+3/k**2), [1, inf])
          1.848936182858244485224927
          >>> 3*sqrt(2)*cosh(pi*sqrt(3)/2)**2*csch(pi*sqrt(2))/pi
          1.848936182858244485224927
  
          >>> nprod(lambda k: (1-1/k**4), [2, inf]); sinh(pi)/(4*pi)
          0.9190194775937444301739244
          0.9190194775937444301739244
  
          >>> nprod(lambda k: (1-1/k**6), [2, inf])
          0.9826842777421925183244759
          >>> (1+cosh(pi*sqrt(3)))/(12*pi**2)
          0.9826842777421925183244759
  
          >>> nprod(lambda k: (1+1/k**2), [2, inf]); sinh(pi)/(2*pi)
          1.838038955187488860347849
          1.838038955187488860347849
  
          >>> nprod(lambda n: (1+1/n)**n * exp(1/(2*n)-1), [1, inf])
          1.447255926890365298959138
          >>> exp(1+euler/2)/sqrt(2*pi)
          1.447255926890365298959138
  
      The following two products are equivalent and can be evaluated in
      terms of a Jacobi theta function. Pi can be replaced by any value
      (as long as convergence is preserved)::
  
          >>> nprod(lambda k: (1-pi**-k)/(1+pi**-k), [1, inf])
          0.3838451207481672404778686
          >>> nprod(lambda k: tanh(k*log(pi)/2), [1, inf])
          0.3838451207481672404778686
          >>> jtheta(4,0,1/pi)
          0.3838451207481672404778686
  
      This product does not have a known closed form value::
  
          >>> nprod(lambda k: (1-1/2**k), [1, inf])
          0.2887880950866024212788997
  
      A product taken from `-\infty`::
  
          >>> nprod(lambda k: 1-k**(-3), [-inf,-2])
          0.8093965973662901095786805
          >>> cosh(pi*sqrt(3)/2)/(3*pi)
          0.8093965973662901095786805
  
      A doubly infinite product::
  
          >>> nprod(lambda k: exp(1/(1+k**2)), [-inf, inf])
          23.41432688231864337420035
          >>> exp(pi/tanh(pi))
          23.41432688231864337420035
  
      A product requiring the use of Euler-Maclaurin summation to compute
      an accurate value::
  
          >>> nprod(lambda k: (1-1/k**2.5), [2, inf], method='e')
          0.696155111336231052898125
  
      **References**
  
      1. [Weisstein]_ http://mathworld.wolfram.com/InfiniteProduct.html
  
      """
      if nsum or ('e' in kwargs.get('method', '')):
          orig = ctx.prec
          try:
              # TODO: we are evaluating log(1+eps) -> eps, which is
              # inaccurate. This currently works because nsum greatly
              # increases the working precision. But we should be
              # more intelligent and handle the precision here.
              ctx.prec += 10
              v = ctx.nsum(lambda n: ctx.ln(f(n)), interval, **kwargs)
          finally:
              ctx.prec = orig
          return +ctx.exp(v)
  
      a, b = ctx._as_points(interval)
      if a == ctx.ninf:
          if b == ctx.inf:
              return f(0) * ctx.nprod(lambda k: f(-k) * f(k), [1, ctx.inf], **kwargs)
          return ctx.nprod(f, [-b, ctx.inf], **kwargs)
      elif b != ctx.inf:
          return ctx.fprod(f(ctx.mpf(k)) for k in xrange(int(a), int(b)+1))
  
      a = int(a)
  
      def update(partial_products, indices):
          if partial_products:
              pprod = partial_products[-1]
          else:
              pprod = ctx.one
          for k in indices:
              pprod = pprod * f(a + ctx.mpf(k))
              partial_products.append(pprod)
  
      return +ctx.adaptive_extrapolation(update, None, kwargs)
  
  
  @defun
  def limit(ctx, f, x, direction=1, exp=False, **kwargs):
      r"""
      Computes an estimate of the limit
  
      .. math ::
  
          \lim_{t \to x} f(t)
  
      where `x` may be finite or infinite.
  
      For finite `x`, :func:`~mpmath.limit` evaluates `f(x + d/n)` for
      consecutive integer values of `n`, where the approach direction
      `d` may be specified using the *direction* keyword argument.
      For infinite `x`, :func:`~mpmath.limit` evaluates values of
      `f(\mathrm{sign}(x) \cdot n)`.
  
      If the approach to the limit is not sufficiently fast to give
      an accurate estimate directly, :func:`~mpmath.limit` attempts to find
      the limit using Richardson extrapolation or the Shanks
      transformation. You can select between these methods using
      the *method* keyword (see documentation of :func:`~mpmath.nsum` for
      more information).
  
      **Options**
  
      The following options are available with essentially the
      same meaning as for :func:`~mpmath.nsum`: *tol*, *method*, *maxterms*,
      *steps*, *verbose*.
  
      If the option *exp=True* is set, `f` will be
      sampled at exponentially spaced points `n = 2^1, 2^2, 2^3, \ldots`
      instead of the linearly spaced points `n = 1, 2, 3, \ldots`.
      This can sometimes improve the rate of convergence so that
      :func:`~mpmath.limit` may return a more accurate answer (and faster).
      However, do note that this can only be used if `f`
      supports fast and accurate evaluation for arguments that
      are extremely close to the limit point (or if infinite,
      very large arguments).
  
      **Examples**
  
      A basic evaluation of a removable singularity::
  
          >>> from mpmath import *
          >>> mp.dps = 30; mp.pretty = True
          >>> limit(lambda x: (x-sin(x))/x**3, 0)
          0.166666666666666666666666666667
  
      Computing the exponential function using its limit definition::
  
          >>> limit(lambda n: (1+3/n)**n, inf)
          20.0855369231876677409285296546
          >>> exp(3)
          20.0855369231876677409285296546
  
      A limit for `\pi`::
  
          >>> f = lambda n: 2**(4*n+1)*fac(n)**4/(2*n+1)/fac(2*n)**2
          >>> limit(f, inf)
          3.14159265358979323846264338328
  
      Calculating the coefficient in Stirling's formula::
  
          >>> limit(lambda n: fac(n) / (sqrt(n)*(n/e)**n), inf)
          2.50662827463100050241576528481
          >>> sqrt(2*pi)
          2.50662827463100050241576528481
  
      Evaluating Euler's constant `\gamma` using the limit representation
  
      .. math ::
  
          \gamma = \lim_{n \rightarrow \infty } \left[ \left(
          \sum_{k=1}^n \frac{1}{k} \right) - \log(n) \right]
  
      (which converges notoriously slowly)::
  
          >>> f = lambda n: sum([mpf(1)/k for k in range(1,int(n)+1)]) - log(n)
          >>> limit(f, inf)
          0.577215664901532860606512090082
          >>> +euler
          0.577215664901532860606512090082
  
      With default settings, the following limit converges too slowly
      to be evaluated accurately. Changing to exponential sampling
      however gives a perfect result::
  
          >>> f = lambda x: sqrt(x**3+x**2)/(sqrt(x**3)+x)
          >>> limit(f, inf)
          0.992831158558330281129249686491
          >>> limit(f, inf, exp=True)
          1.0
  
      """
  
      if ctx.isinf(x):
          direction = ctx.sign(x)
          g = lambda k: f(ctx.mpf(k+1)*direction)
      else:
          direction *= ctx.one
          g = lambda k: f(x + direction/(k+1))
      if exp:
          h = g
          g = lambda k: h(2**k)
  
      def update(values, indices):
          for k in indices:
              values.append(g(k+1))
  
      # XXX: steps used by nsum don't work well
      if not 'steps' in kwargs:
          kwargs['steps'] = [10]
  
      return +ctx.adaptive_extrapolation(update, None, kwargs)
  