Coverage for /usr/lib/python3/dist-packages/mpmath/calculus/extrapolation.py: 8%
587 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1try:
2 from itertools import izip
3except ImportError:
4 izip = zip
6from ..libmp.backend import xrange
7from .calculus import defun
9try:
10 next = next
11except NameError:
12 next = lambda _: _.next()
14@defun
15def richardson(ctx, seq):
16 r"""
17 Given a list ``seq`` of the first `N` elements of a slowly convergent
18 infinite sequence, :func:`~mpmath.richardson` computes the `N`-term
19 Richardson extrapolate for the limit.
21 :func:`~mpmath.richardson` returns `(v, c)` where `v` is the estimated
22 limit and `c` is the magnitude of the largest weight used during the
23 computation. The weight provides an estimate of the precision
24 lost to cancellation. Due to cancellation effects, the sequence must
25 be typically be computed at a much higher precision than the target
26 accuracy of the extrapolation.
28 **Applicability and issues**
30 The `N`-step Richardson extrapolation algorithm used by
31 :func:`~mpmath.richardson` is described in [1].
33 Richardson extrapolation only works for a specific type of sequence,
34 namely one converging like partial sums of
35 `P(1)/Q(1) + P(2)/Q(2) + \ldots` where `P` and `Q` are polynomials.
36 When the sequence does not convergence at such a rate
37 :func:`~mpmath.richardson` generally produces garbage.
39 Richardson extrapolation has the advantage of being fast: the `N`-term
40 extrapolate requires only `O(N)` arithmetic operations, and usually
41 produces an estimate that is accurate to `O(N)` digits. Contrast with
42 the Shanks transformation (see :func:`~mpmath.shanks`), which requires
43 `O(N^2)` operations.
45 :func:`~mpmath.richardson` is unable to produce an estimate for the
46 approximation error. One way to estimate the error is to perform
47 two extrapolations with slightly different `N` and comparing the
48 results.
50 Richardson extrapolation does not work for oscillating sequences.
51 As a simple workaround, :func:`~mpmath.richardson` detects if the last
52 three elements do not differ monotonically, and in that case
53 applies extrapolation only to the even-index elements.
55 **Example**
57 Applying Richardson extrapolation to the Leibniz series for `\pi`::
59 >>> from mpmath import *
60 >>> mp.dps = 30; mp.pretty = True
61 >>> S = [4*sum(mpf(-1)**n/(2*n+1) for n in range(m))
62 ... for m in range(1,30)]
63 >>> v, c = richardson(S[:10])
64 >>> v
65 3.2126984126984126984126984127
66 >>> nprint([v-pi, c])
67 [0.0711058, 2.0]
69 >>> v, c = richardson(S[:30])
70 >>> v
71 3.14159265468624052829954206226
72 >>> nprint([v-pi, c])
73 [1.09645e-9, 20833.3]
75 **References**
77 1. [BenderOrszag]_ pp. 375-376
79 """
80 if len(seq) < 3:
81 raise ValueError("seq should be of minimum length 3")
82 if ctx.sign(seq[-1]-seq[-2]) != ctx.sign(seq[-2]-seq[-3]):
83 seq = seq[::2]
84 N = len(seq)//2-1
85 s = ctx.zero
86 # The general weight is c[k] = (N+k)**N * (-1)**(k+N) / k! / (N-k)!
87 # To avoid repeated factorials, we simplify the quotient
88 # of successive weights to obtain a recurrence relation
89 c = (-1)**N * N**N / ctx.mpf(ctx._ifac(N))
90 maxc = 1
91 for k in xrange(N+1):
92 s += c * seq[N+k]
93 maxc = max(abs(c), maxc)
94 c *= (k-N)*ctx.mpf(k+N+1)**N
95 c /= ((1+k)*ctx.mpf(k+N)**N)
96 return s, maxc
98@defun
99def shanks(ctx, seq, table=None, randomized=False):
100 r"""
101 Given a list ``seq`` of the first `N` elements of a slowly
102 convergent infinite sequence `(A_k)`, :func:`~mpmath.shanks` computes the iterated
103 Shanks transformation `S(A), S(S(A)), \ldots, S^{N/2}(A)`. The Shanks
104 transformation often provides strong convergence acceleration,
105 especially if the sequence is oscillating.
107 The iterated Shanks transformation is computed using the Wynn
108 epsilon algorithm (see [1]). :func:`~mpmath.shanks` returns the full
109 epsilon table generated by Wynn's algorithm, which can be read
110 off as follows:
112 * The table is a list of lists forming a lower triangular matrix,
113 where higher row and column indices correspond to more accurate
114 values.
115 * The columns with even index hold dummy entries (required for the
116 computation) and the columns with odd index hold the actual
117 extrapolates.
118 * The last element in the last row is typically the most
119 accurate estimate of the limit.
120 * The difference to the third last element in the last row
121 provides an estimate of the approximation error.
122 * The magnitude of the second last element provides an estimate
123 of the numerical accuracy lost to cancellation.
125 For convenience, so the extrapolation is stopped at an odd index
126 so that ``shanks(seq)[-1][-1]`` always gives an estimate of the
127 limit.
129 Optionally, an existing table can be passed to :func:`~mpmath.shanks`.
130 This can be used to efficiently extend a previous computation after
131 new elements have been appended to the sequence. The table will
132 then be updated in-place.
134 **The Shanks transformation**
136 The Shanks transformation is defined as follows (see [2]): given
137 the input sequence `(A_0, A_1, \ldots)`, the transformed sequence is
138 given by
140 .. math ::
142 S(A_k) = \frac{A_{k+1}A_{k-1}-A_k^2}{A_{k+1}+A_{k-1}-2 A_k}
144 The Shanks transformation gives the exact limit `A_{\infty}` in a
145 single step if `A_k = A + a q^k`. Note in particular that it
146 extrapolates the exact sum of a geometric series in a single step.
148 Applying the Shanks transformation once often improves convergence
149 substantially for an arbitrary sequence, but the optimal effect is
150 obtained by applying it iteratively:
151 `S(S(A_k)), S(S(S(A_k))), \ldots`.
153 Wynn's epsilon algorithm provides an efficient way to generate
154 the table of iterated Shanks transformations. It reduces the
155 computation of each element to essentially a single division, at
156 the cost of requiring dummy elements in the table. See [1] for
157 details.
159 **Precision issues**
161 Due to cancellation effects, the sequence must be typically be
162 computed at a much higher precision than the target accuracy
163 of the extrapolation.
165 If the Shanks transformation converges to the exact limit (such
166 as if the sequence is a geometric series), then a division by
167 zero occurs. By default, :func:`~mpmath.shanks` handles this case by
168 terminating the iteration and returning the table it has
169 generated so far. With *randomized=True*, it will instead
170 replace the zero by a pseudorandom number close to zero.
171 (TODO: find a better solution to this problem.)
173 **Examples**
175 We illustrate by applying Shanks transformation to the Leibniz
176 series for `\pi`::
178 >>> from mpmath import *
179 >>> mp.dps = 50
180 >>> S = [4*sum(mpf(-1)**n/(2*n+1) for n in range(m))
181 ... for m in range(1,30)]
182 >>>
183 >>> T = shanks(S[:7])
184 >>> for row in T:
185 ... nprint(row)
186 ...
187 [-0.75]
188 [1.25, 3.16667]
189 [-1.75, 3.13333, -28.75]
190 [2.25, 3.14524, 82.25, 3.14234]
191 [-2.75, 3.13968, -177.75, 3.14139, -969.937]
192 [3.25, 3.14271, 327.25, 3.14166, 3515.06, 3.14161]
194 The extrapolated accuracy is about 4 digits, and about 4 digits
195 may have been lost due to cancellation::
197 >>> L = T[-1]
198 >>> nprint([abs(L[-1] - pi), abs(L[-1] - L[-3]), abs(L[-2])])
199 [2.22532e-5, 4.78309e-5, 3515.06]
201 Now we extend the computation::
203 >>> T = shanks(S[:25], T)
204 >>> L = T[-1]
205 >>> nprint([abs(L[-1] - pi), abs(L[-1] - L[-3]), abs(L[-2])])
206 [3.75527e-19, 1.48478e-19, 2.96014e+17]
208 The value for pi is now accurate to 18 digits. About 18 digits may
209 also have been lost to cancellation.
211 Here is an example with a geometric series, where the convergence
212 is immediate (the sum is exactly 1)::
214 >>> mp.dps = 15
215 >>> for row in shanks([0.5, 0.75, 0.875, 0.9375, 0.96875]):
216 ... nprint(row)
217 [4.0]
218 [8.0, 1.0]
220 **References**
222 1. [GravesMorris]_
224 2. [BenderOrszag]_ pp. 368-375
226 """
227 if len(seq) < 2:
228 raise ValueError("seq should be of minimum length 2")
229 if table:
230 START = len(table)
231 else:
232 START = 0
233 table = []
234 STOP = len(seq) - 1
235 if STOP & 1:
236 STOP -= 1
237 one = ctx.one
238 eps = +ctx.eps
239 if randomized:
240 from random import Random
241 rnd = Random()
242 rnd.seed(START)
243 for i in xrange(START, STOP):
244 row = []
245 for j in xrange(i+1):
246 if j == 0:
247 a, b = 0, seq[i+1]-seq[i]
248 else:
249 if j == 1:
250 a = seq[i]
251 else:
252 a = table[i-1][j-2]
253 b = row[j-1] - table[i-1][j-1]
254 if not b:
255 if randomized:
256 b = (1 + rnd.getrandbits(10))*eps
257 elif i & 1:
258 return table[:-1]
259 else:
260 return table
261 row.append(a + one/b)
262 table.append(row)
263 return table
266class levin_class:
267 # levin: Copyright 2013 Timo Hartmann (thartmann15 at gmail.com)
268 r"""
269 This interface implements Levin's (nonlinear) sequence transformation for
270 convergence acceleration and summation of divergent series. It performs
271 better than the Shanks/Wynn-epsilon algorithm for logarithmic convergent
272 or alternating divergent series.
274 Let *A* be the series we want to sum:
276 .. math ::
278 A = \sum_{k=0}^{\infty} a_k
280 Attention: all `a_k` must be non-zero!
282 Let `s_n` be the partial sums of this series:
284 .. math ::
286 s_n = \sum_{k=0}^n a_k.
288 **Methods**
290 Calling ``levin`` returns an object with the following methods.
292 ``update(...)`` works with the list of individual terms `a_k` of *A*, and
293 ``update_step(...)`` works with the list of partial sums `s_k` of *A*:
295 .. code ::
297 v, e = ...update([a_0, a_1,..., a_k])
298 v, e = ...update_psum([s_0, s_1,..., s_k])
300 ``step(...)`` works with the individual terms `a_k` and ``step_psum(...)``
301 works with the partial sums `s_k`:
303 .. code ::
305 v, e = ...step(a_k)
306 v, e = ...step_psum(s_k)
308 *v* is the current estimate for *A*, and *e* is an error estimate which is
309 simply the difference between the current estimate and the last estimate.
310 One should not mix ``update``, ``update_psum``, ``step`` and ``step_psum``.
312 **A word of caution**
314 One can only hope for good results (i.e. convergence acceleration or
315 resummation) if the `s_n` have some well defind asymptotic behavior for
316 large `n` and are not erratic or random. Furthermore one usually needs very
317 high working precision because of the numerical cancellation. If the working
318 precision is insufficient, levin may produce silently numerical garbage.
319 Furthermore even if the Levin-transformation converges, in the general case
320 there is no proof that the result is mathematically sound. Only for very
321 special classes of problems one can prove that the Levin-transformation
322 converges to the expected result (for example Stieltjes-type integrals).
323 Furthermore the Levin-transform is quite expensive (i.e. slow) in comparison
324 to Shanks/Wynn-epsilon, Richardson & co.
325 In summary one can say that the Levin-transformation is powerful but
326 unreliable and that it may need a copious amount of working precision.
328 The Levin transform has several variants differing in the choice of weights.
329 Some variants are better suited for the possible flavours of convergence
330 behaviour of *A* than other variants:
332 .. code ::
334 convergence behaviour levin-u levin-t levin-v shanks/wynn-epsilon
336 logarithmic + - + -
337 linear + + + +
338 alternating divergent + + + +
340 "+" means the variant is suitable,"-" means the variant is not suitable;
341 for comparison the Shanks/Wynn-epsilon transform is listed, too.
343 The variant is controlled though the variant keyword (i.e. ``variant="u"``,
344 ``variant="t"`` or ``variant="v"``). Overall "u" is probably the best choice.
346 Finally it is possible to use the Sidi-S transform instead of the Levin transform
347 by using the keyword ``method='sidi'``. The Sidi-S transform works better than the
348 Levin transformation for some divergent series (see the examples).
350 Parameters:
352 .. code ::
354 method "levin" or "sidi" chooses either the Levin or the Sidi-S transformation
355 variant "u","t" or "v" chooses the weight variant.
357 The Levin transform is also accessible through the nsum interface.
358 ``method="l"`` or ``method="levin"`` select the normal Levin transform while
359 ``method="sidi"``
360 selects the Sidi-S transform. The variant is in both cases selected through the
361 levin_variant keyword. The stepsize in :func:`~mpmath.nsum` must not be chosen too large, otherwise
362 it will miss the point where the Levin transform converges resulting in numerical
363 overflow/garbage. For highly divergent series a copious amount of working precision
364 must be chosen.
366 **Examples**
368 First we sum the zeta function::
370 >>> from mpmath import mp
371 >>> mp.prec = 53
372 >>> eps = mp.mpf(mp.eps)
373 >>> with mp.extraprec(2 * mp.prec): # levin needs a high working precision
374 ... L = mp.levin(method = "levin", variant = "u")
375 ... S, s, n = [], 0, 1
376 ... while 1:
377 ... s += mp.one / (n * n)
378 ... n += 1
379 ... S.append(s)
380 ... v, e = L.update_psum(S)
381 ... if e < eps:
382 ... break
383 ... if n > 1000: raise RuntimeError("iteration limit exceeded")
384 >>> print(mp.chop(v - mp.pi ** 2 / 6))
385 0.0
386 >>> w = mp.nsum(lambda n: 1 / (n*n), [1, mp.inf], method = "levin", levin_variant = "u")
387 >>> print(mp.chop(v - w))
388 0.0
390 Now we sum the zeta function outside its range of convergence
391 (attention: This does not work at the negative integers!)::
393 >>> eps = mp.mpf(mp.eps)
394 >>> with mp.extraprec(2 * mp.prec): # levin needs a high working precision
395 ... L = mp.levin(method = "levin", variant = "v")
396 ... A, n = [], 1
397 ... while 1:
398 ... s = mp.mpf(n) ** (2 + 3j)
399 ... n += 1
400 ... A.append(s)
401 ... v, e = L.update(A)
402 ... if e < eps:
403 ... break
404 ... if n > 1000: raise RuntimeError("iteration limit exceeded")
405 >>> print(mp.chop(v - mp.zeta(-2-3j)))
406 0.0
407 >>> w = mp.nsum(lambda n: n ** (2 + 3j), [1, mp.inf], method = "levin", levin_variant = "v")
408 >>> print(mp.chop(v - w))
409 0.0
411 Now we sum the divergent asymptotic expansion of an integral related to the
412 exponential integral (see also [2] p.373). The Sidi-S transform works best here::
414 >>> z = mp.mpf(10)
415 >>> exact = mp.quad(lambda x: mp.exp(-x)/(1+x/z),[0,mp.inf])
416 >>> # exact = z * mp.exp(z) * mp.expint(1,z) # this is the symbolic expression for the integral
417 >>> eps = mp.mpf(mp.eps)
418 >>> with mp.extraprec(2 * mp.prec): # high working precisions are mandatory for divergent resummation
419 ... L = mp.levin(method = "sidi", variant = "t")
420 ... n = 0
421 ... while 1:
422 ... s = (-1)**n * mp.fac(n) * z ** (-n)
423 ... v, e = L.step(s)
424 ... n += 1
425 ... if e < eps:
426 ... break
427 ... if n > 1000: raise RuntimeError("iteration limit exceeded")
428 >>> print(mp.chop(v - exact))
429 0.0
430 >>> w = mp.nsum(lambda n: (-1) ** n * mp.fac(n) * z ** (-n), [0, mp.inf], method = "sidi", levin_variant = "t")
431 >>> print(mp.chop(v - w))
432 0.0
434 Another highly divergent integral is also summable::
436 >>> z = mp.mpf(2)
437 >>> eps = mp.mpf(mp.eps)
438 >>> exact = mp.quad(lambda x: mp.exp( -x * x / 2 - z * x ** 4), [0,mp.inf]) * 2 / mp.sqrt(2 * mp.pi)
439 >>> # 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
440 >>> with mp.extraprec(7 * mp.prec): # we need copious amount of precision to sum this highly divergent series
441 ... L = mp.levin(method = "levin", variant = "t")
442 ... n, s = 0, 0
443 ... while 1:
444 ... s += (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n))
445 ... n += 1
446 ... v, e = L.step_psum(s)
447 ... if e < eps:
448 ... break
449 ... if n > 1000: raise RuntimeError("iteration limit exceeded")
450 >>> print(mp.chop(v - exact))
451 0.0
452 >>> w = mp.nsum(lambda n: (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n)),
453 ... [0, mp.inf], method = "levin", levin_variant = "t", workprec = 8*mp.prec, steps = [2] + [1 for x in xrange(1000)])
454 >>> print(mp.chop(v - w))
455 0.0
457 These examples run with 15-20 decimal digits precision. For higher precision the
458 working precision must be raised.
460 **Examples for nsum**
462 Here we calculate Euler's constant as the constant term in the Laurent
463 expansion of `\zeta(s)` at `s=1`. This sum converges extremly slowly because of
464 the logarithmic convergence behaviour of the Dirichlet series for zeta::
466 >>> mp.dps = 30
467 >>> z = mp.mpf(10) ** (-10)
468 >>> a = mp.nsum(lambda n: n**(-(1+z)), [1, mp.inf], method = "l") - 1 / z
469 >>> print(mp.chop(a - mp.euler, tol = 1e-10))
470 0.0
472 The Sidi-S transform performs excellently for the alternating series of `\log(2)`::
474 >>> a = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "sidi")
475 >>> print(mp.chop(a - mp.log(2)))
476 0.0
478 Hypergeometric series can also be summed outside their range of convergence.
479 The stepsize in :func:`~mpmath.nsum` must not be chosen too large, otherwise it will miss the
480 point where the Levin transform converges resulting in numerical overflow/garbage::
482 >>> z = 2 + 1j
483 >>> exact = mp.hyp2f1(2 / mp.mpf(3), 4 / mp.mpf(3), 1 / mp.mpf(3), z)
484 >>> 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))
485 >>> v = mp.nsum(f, [0, mp.inf], method = "levin", steps = [10 for x in xrange(1000)])
486 >>> print(mp.chop(exact-v))
487 0.0
489 References:
491 [1] E.J. Weniger - "Nonlinear Sequence Transformations for the Acceleration of
492 Convergence and the Summation of Divergent Series" arXiv:math/0306302
494 [2] A. Sidi - "Pratical Extrapolation Methods"
496 [3] H.H.H. Homeier - "Scalar Levin-Type Sequence Transformations" arXiv:math/0005209
498 """
500 def __init__(self, method = "levin", variant = "u"):
501 self.variant = variant
502 self.n = 0
503 self.a0 = 0
504 self.theta = 1
505 self.A = []
506 self.B = []
507 self.last = 0
508 self.last_s = False
510 if method == "levin":
511 self.factor = self.factor_levin
512 elif method == "sidi":
513 self.factor = self.factor_sidi
514 else:
515 raise ValueError("levin: unknown method \"%s\"" % method)
517 def factor_levin(self, i):
518 # original levin
519 # [1] p.50,e.7.5-7 (with n-j replaced by i)
520 return (self.theta + i) * (self.theta + self.n - 1) ** (self.n - i - 2) / self.ctx.mpf(self.theta + self.n) ** (self.n - i - 1)
522 def factor_sidi(self, i):
523 # sidi analogon to levin (factorial series)
524 # [1] p.59,e.8.3-16 (with n-j replaced by i)
525 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))
527 def run(self, s, a0, a1 = 0):
528 if self.variant=="t":
529 # levin t
530 w=a0
531 elif self.variant=="u":
532 # levin u
533 w=a0*(self.theta+self.n)
534 elif self.variant=="v":
535 # levin v
536 w=a0*a1/(a0-a1)
537 else:
538 assert False, "unknown variant"
540 if w==0:
541 raise ValueError("levin: zero weight")
543 self.A.append(s/w)
544 self.B.append(1/w)
546 for i in range(self.n-1,-1,-1):
547 if i==self.n-1:
548 f=1
549 else:
550 f=self.factor(i)
552 self.A[i]=self.A[i+1]-f*self.A[i]
553 self.B[i]=self.B[i+1]-f*self.B[i]
555 self.n+=1
557 ###########################################################################
559 def update_psum(self,S):
560 """
561 This routine applies the convergence acceleration to the list of partial sums.
563 A = sum(a_k, k = 0..infinity)
564 s_n = sum(a_k, k = 0..n)
566 v, e = ...update_psum([s_0, s_1,..., s_k])
568 output:
569 v current estimate of the series A
570 e an error estimate which is simply the difference between the current
571 estimate and the last estimate.
572 """
574 if self.variant!="v":
575 if self.n==0:
576 self.run(S[0],S[0])
577 while self.n<len(S):
578 self.run(S[self.n],S[self.n]-S[self.n-1])
579 else:
580 if len(S)==1:
581 self.last=0
582 return S[0],abs(S[0])
584 if self.n==0:
585 self.a1=S[1]-S[0]
586 self.run(S[0],S[0],self.a1)
588 while self.n<len(S)-1:
589 na1=S[self.n+1]-S[self.n]
590 self.run(S[self.n],self.a1,na1)
591 self.a1=na1
593 value=self.A[0]/self.B[0]
594 err=abs(value-self.last)
595 self.last=value
597 return value,err
599 def update(self,X):
600 """
601 This routine applies the convergence acceleration to the list of individual terms.
603 A = sum(a_k, k = 0..infinity)
605 v, e = ...update([a_0, a_1,..., a_k])
607 output:
608 v current estimate of the series A
609 e an error estimate which is simply the difference between the current
610 estimate and the last estimate.
611 """
613 if self.variant!="v":
614 if self.n==0:
615 self.s=X[0]
616 self.run(self.s,X[0])
617 while self.n<len(X):
618 self.s+=X[self.n]
619 self.run(self.s,X[self.n])
620 else:
621 if len(X)==1:
622 self.last=0
623 return X[0],abs(X[0])
625 if self.n==0:
626 self.s=X[0]
627 self.run(self.s,X[0],X[1])
629 while self.n<len(X)-1:
630 self.s+=X[self.n]
631 self.run(self.s,X[self.n],X[self.n+1])
633 value=self.A[0]/self.B[0]
634 err=abs(value-self.last)
635 self.last=value
637 return value,err
639 ###########################################################################
641 def step_psum(self,s):
642 """
643 This routine applies the convergence acceleration to the partial sums.
645 A = sum(a_k, k = 0..infinity)
646 s_n = sum(a_k, k = 0..n)
648 v, e = ...step_psum(s_k)
650 output:
651 v current estimate of the series A
652 e an error estimate which is simply the difference between the current
653 estimate and the last estimate.
654 """
656 if self.variant!="v":
657 if self.n==0:
658 self.last_s=s
659 self.run(s,s)
660 else:
661 self.run(s,s-self.last_s)
662 self.last_s=s
663 else:
664 if isinstance(self.last_s,bool):
665 self.last_s=s
666 self.last_w=s
667 self.last=0
668 return s,abs(s)
670 na1=s-self.last_s
671 self.run(self.last_s,self.last_w,na1)
672 self.last_w=na1
673 self.last_s=s
675 value=self.A[0]/self.B[0]
676 err=abs(value-self.last)
677 self.last=value
679 return value,err
681 def step(self,x):
682 """
683 This routine applies the convergence acceleration to the individual terms.
685 A = sum(a_k, k = 0..infinity)
687 v, e = ...step(a_k)
689 output:
690 v current estimate of the series A
691 e an error estimate which is simply the difference between the current
692 estimate and the last estimate.
693 """
695 if self.variant!="v":
696 if self.n==0:
697 self.s=x
698 self.run(self.s,x)
699 else:
700 self.s+=x
701 self.run(self.s,x)
702 else:
703 if isinstance(self.last_s,bool):
704 self.last_s=x
705 self.s=0
706 self.last=0
707 return x,abs(x)
709 self.s+=self.last_s
710 self.run(self.s,self.last_s,x)
711 self.last_s=x
713 value=self.A[0]/self.B[0]
714 err=abs(value-self.last)
715 self.last=value
717 return value,err
719def levin(ctx, method = "levin", variant = "u"):
720 L = levin_class(method = method, variant = variant)
721 L.ctx = ctx
722 return L
724levin.__doc__ = levin_class.__doc__
725defun(levin)
728class cohen_alt_class:
729 # cohen_alt: Copyright 2013 Timo Hartmann (thartmann15 at gmail.com)
730 r"""
731 This interface implements the convergence acceleration of alternating series
732 as described in H. Cohen, F.R. Villegas, D. Zagier - "Convergence Acceleration
733 of Alternating Series". This series transformation works only well if the
734 individual terms of the series have an alternating sign. It belongs to the
735 class of linear series transformations (in contrast to the Shanks/Wynn-epsilon
736 or Levin transform). This series transformation is also able to sum some types
737 of divergent series. See the paper under which conditions this resummation is
738 mathematical sound.
740 Let *A* be the series we want to sum:
742 .. math ::
744 A = \sum_{k=0}^{\infty} a_k
746 Let `s_n` be the partial sums of this series:
748 .. math ::
750 s_n = \sum_{k=0}^n a_k.
753 **Interface**
755 Calling ``cohen_alt`` returns an object with the following methods.
757 Then ``update(...)`` works with the list of individual terms `a_k` and
758 ``update_psum(...)`` works with the list of partial sums `s_k`:
760 .. code ::
762 v, e = ...update([a_0, a_1,..., a_k])
763 v, e = ...update_psum([s_0, s_1,..., s_k])
765 *v* is the current estimate for *A*, and *e* is an error estimate which is
766 simply the difference between the current estimate and the last estimate.
768 **Examples**
770 Here we compute the alternating zeta function using ``update_psum``::
772 >>> from mpmath import mp
773 >>> AC = mp.cohen_alt()
774 >>> S, s, n = [], 0, 1
775 >>> while 1:
776 ... s += -((-1) ** n) * mp.one / (n * n)
777 ... n += 1
778 ... S.append(s)
779 ... v, e = AC.update_psum(S)
780 ... if e < mp.eps:
781 ... break
782 ... if n > 1000: raise RuntimeError("iteration limit exceeded")
783 >>> print(mp.chop(v - mp.pi ** 2 / 12))
784 0.0
786 Here we compute the product `\prod_{n=1}^{\infty} \Gamma(1+1/(2n-1)) / \Gamma(1+1/(2n))`::
788 >>> A = []
789 >>> AC = mp.cohen_alt()
790 >>> n = 1
791 >>> while 1:
792 ... A.append( mp.loggamma(1 + mp.one / (2 * n - 1)))
793 ... A.append(-mp.loggamma(1 + mp.one / (2 * n)))
794 ... n += 1
795 ... v, e = AC.update(A)
796 ... if e < mp.eps:
797 ... break
798 ... if n > 1000: raise RuntimeError("iteration limit exceeded")
799 >>> v = mp.exp(v)
800 >>> print(mp.chop(v - 1.06215090557106, tol = 1e-12))
801 0.0
803 ``cohen_alt`` is also accessible through the :func:`~mpmath.nsum` interface::
805 >>> v = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "a")
806 >>> print(mp.chop(v - mp.log(2)))
807 0.0
808 >>> v = mp.nsum(lambda n: (-1)**n / (2 * n + 1), [0, mp.inf], method = "a")
809 >>> print(mp.chop(v - mp.pi / 4))
810 0.0
811 >>> v = mp.nsum(lambda n: (-1)**n * mp.log(n) * n, [1, mp.inf], method = "a")
812 >>> print(mp.chop(v - mp.diff(lambda s: mp.altzeta(s), -1)))
813 0.0
815 """
817 def __init__(self):
818 self.last=0
820 def update(self, A):
821 """
822 This routine applies the convergence acceleration to the list of individual terms.
824 A = sum(a_k, k = 0..infinity)
826 v, e = ...update([a_0, a_1,..., a_k])
828 output:
829 v current estimate of the series A
830 e an error estimate which is simply the difference between the current
831 estimate and the last estimate.
832 """
834 n = len(A)
835 d = (3 + self.ctx.sqrt(8)) ** n
836 d = (d + 1 / d) / 2
837 b = -self.ctx.one
838 c = -d
839 s = 0
841 for k in xrange(n):
842 c = b - c
843 if k % 2 == 0:
844 s = s + c * A[k]
845 else:
846 s = s - c * A[k]
847 b = 2 * (k + n) * (k - n) * b / ((2 * k + 1) * (k + self.ctx.one))
849 value = s / d
851 err = abs(value - self.last)
852 self.last = value
854 return value, err
856 def update_psum(self, S):
857 """
858 This routine applies the convergence acceleration to the list of partial sums.
860 A = sum(a_k, k = 0..infinity)
861 s_n = sum(a_k ,k = 0..n)
863 v, e = ...update_psum([s_0, s_1,..., s_k])
865 output:
866 v current estimate of the series A
867 e an error estimate which is simply the difference between the current
868 estimate and the last estimate.
869 """
871 n = len(S)
872 d = (3 + self.ctx.sqrt(8)) ** n
873 d = (d + 1 / d) / 2
874 b = self.ctx.one
875 s = 0
877 for k in xrange(n):
878 b = 2 * (n + k) * (n - k) * b / ((2 * k + 1) * (k + self.ctx.one))
879 s += b * S[k]
881 value = s / d
883 err = abs(value - self.last)
884 self.last = value
886 return value, err
888def cohen_alt(ctx):
889 L = cohen_alt_class()
890 L.ctx = ctx
891 return L
893cohen_alt.__doc__ = cohen_alt_class.__doc__
894defun(cohen_alt)
897@defun
898def sumap(ctx, f, interval, integral=None, error=False):
899 r"""
900 Evaluates an infinite series of an analytic summand *f* using the
901 Abel-Plana formula
903 .. math ::
905 \sum_{k=0}^{\infty} f(k) = \int_0^{\infty} f(t) dt + \frac{1}{2} f(0) +
906 i \int_0^{\infty} \frac{f(it)-f(-it)}{e^{2\pi t}-1} dt.
908 Unlike the Euler-Maclaurin formula (see :func:`~mpmath.sumem`),
909 the Abel-Plana formula does not require derivatives. However,
910 it only works when `|f(it)-f(-it)|` does not
911 increase too rapidly with `t`.
913 **Examples**
915 The Abel-Plana formula is particularly useful when the summand
916 decreases like a power of `k`; for example when the sum is a pure
917 zeta function::
919 >>> from mpmath import *
920 >>> mp.dps = 25; mp.pretty = True
921 >>> sumap(lambda k: 1/k**2.5, [1,inf])
922 1.34148725725091717975677
923 >>> zeta(2.5)
924 1.34148725725091717975677
925 >>> sumap(lambda k: 1/(k+1j)**(2.5+2.5j), [1,inf])
926 (-3.385361068546473342286084 - 0.7432082105196321803869551j)
927 >>> zeta(2.5+2.5j, 1+1j)
928 (-3.385361068546473342286084 - 0.7432082105196321803869551j)
930 If the series is alternating, numerical quadrature along the real
931 line is likely to give poor results, so it is better to evaluate
932 the first term symbolically whenever possible:
934 >>> n=3; z=-0.75
935 >>> I = expint(n,-log(z))
936 >>> chop(sumap(lambda k: z**k / k**n, [1,inf], integral=I))
937 -0.6917036036904594510141448
938 >>> polylog(n,z)
939 -0.6917036036904594510141448
941 """
942 prec = ctx.prec
943 try:
944 ctx.prec += 10
945 a, b = interval
946 if b != ctx.inf:
947 raise ValueError("b should be equal to ctx.inf")
948 g = lambda x: f(x+a)
949 if integral is None:
950 i1, err1 = ctx.quad(g, [0,ctx.inf], error=True)
951 else:
952 i1, err1 = integral, 0
953 j = ctx.j
954 p = ctx.pi * 2
955 if ctx._is_real_type(i1):
956 h = lambda t: -2 * ctx.im(g(j*t)) / ctx.expm1(p*t)
957 else:
958 h = lambda t: j*(g(j*t)-g(-j*t)) / ctx.expm1(p*t)
959 i2, err2 = ctx.quad(h, [0,ctx.inf], error=True)
960 err = err1+err2
961 v = i1+i2+0.5*g(ctx.mpf(0))
962 finally:
963 ctx.prec = prec
964 if error:
965 return +v, err
966 return +v
969@defun
970def sumem(ctx, f, interval, tol=None, reject=10, integral=None,
971 adiffs=None, bdiffs=None, verbose=False, error=False,
972 _fast_abort=False):
973 r"""
974 Uses the Euler-Maclaurin formula to compute an approximation accurate
975 to within ``tol`` (which defaults to the present epsilon) of the sum
977 .. math ::
979 S = \sum_{k=a}^b f(k)
981 where `(a,b)` are given by ``interval`` and `a` or `b` may be
982 infinite. The approximation is
984 .. math ::
986 S \sim \int_a^b f(x) \,dx + \frac{f(a)+f(b)}{2} +
987 \sum_{k=1}^{\infty} \frac{B_{2k}}{(2k)!}
988 \left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right).
990 The last sum in the Euler-Maclaurin formula is not generally
991 convergent (a notable exception is if `f` is a polynomial, in
992 which case Euler-Maclaurin actually gives an exact result).
994 The summation is stopped as soon as the quotient between two
995 consecutive terms falls below *reject*. That is, by default
996 (*reject* = 10), the summation is continued as long as each
997 term adds at least one decimal.
999 Although not convergent, convergence to a given tolerance can
1000 often be "forced" if `b = \infty` by summing up to `a+N` and then
1001 applying the Euler-Maclaurin formula to the sum over the range
1002 `(a+N+1, \ldots, \infty)`. This procedure is implemented by
1003 :func:`~mpmath.nsum`.
1005 By default numerical quadrature and differentiation is used.
1006 If the symbolic values of the integral and endpoint derivatives
1007 are known, it is more efficient to pass the value of the
1008 integral explicitly as ``integral`` and the derivatives
1009 explicitly as ``adiffs`` and ``bdiffs``. The derivatives
1010 should be given as iterables that yield
1011 `f(a), f'(a), f''(a), \ldots` (and the equivalent for `b`).
1013 **Examples**
1015 Summation of an infinite series, with automatic and symbolic
1016 integral and derivative values (the second should be much faster)::
1018 >>> from mpmath import *
1019 >>> mp.dps = 50; mp.pretty = True
1020 >>> sumem(lambda n: 1/n**2, [32, inf])
1021 0.03174336652030209012658168043874142714132886413417
1022 >>> I = mpf(1)/32
1023 >>> D = adiffs=((-1)**n*fac(n+1)*32**(-2-n) for n in range(999))
1024 >>> sumem(lambda n: 1/n**2, [32, inf], integral=I, adiffs=D)
1025 0.03174336652030209012658168043874142714132886413417
1027 An exact evaluation of a finite polynomial sum::
1029 >>> sumem(lambda n: n**5-12*n**2+3*n, [-100000, 200000])
1030 10500155000624963999742499550000.0
1031 >>> print(sum(n**5-12*n**2+3*n for n in range(-100000, 200001)))
1032 10500155000624963999742499550000
1034 """
1035 tol = tol or +ctx.eps
1036 interval = ctx._as_points(interval)
1037 a = ctx.convert(interval[0])
1038 b = ctx.convert(interval[-1])
1039 err = ctx.zero
1040 prev = 0
1041 M = 10000
1042 if a == ctx.ninf: adiffs = (0 for n in xrange(M))
1043 else: adiffs = adiffs or ctx.diffs(f, a)
1044 if b == ctx.inf: bdiffs = (0 for n in xrange(M))
1045 else: bdiffs = bdiffs or ctx.diffs(f, b)
1046 orig = ctx.prec
1047 #verbose = 1
1048 try:
1049 ctx.prec += 10
1050 s = ctx.zero
1051 for k, (da, db) in enumerate(izip(adiffs, bdiffs)):
1052 if k & 1:
1053 term = (db-da) * ctx.bernoulli(k+1) / ctx.factorial(k+1)
1054 mag = abs(term)
1055 if verbose:
1056 print("term", k, "magnitude =", ctx.nstr(mag))
1057 if k > 4 and mag < tol:
1058 s += term
1059 break
1060 elif k > 4 and abs(prev) / mag < reject:
1061 err += mag
1062 if _fast_abort:
1063 return [s, (s, err)][error]
1064 if verbose:
1065 print("Failed to converge")
1066 break
1067 else:
1068 s += term
1069 prev = term
1070 # Endpoint correction
1071 if a != ctx.ninf: s += f(a)/2
1072 if b != ctx.inf: s += f(b)/2
1073 # Tail integral
1074 if verbose:
1075 print("Integrating f(x) from x = %s to %s" % (ctx.nstr(a), ctx.nstr(b)))
1076 if integral:
1077 s += integral
1078 else:
1079 integral, ierr = ctx.quad(f, interval, error=True)
1080 if verbose:
1081 print("Integration error:", ierr)
1082 s += integral
1083 err += ierr
1084 finally:
1085 ctx.prec = orig
1086 if error:
1087 return s, err
1088 else:
1089 return s
1091@defun
1092def adaptive_extrapolation(ctx, update, emfun, kwargs):
1093 option = kwargs.get
1094 if ctx._fixed_precision:
1095 tol = option('tol', ctx.eps*2**10)
1096 else:
1097 tol = option('tol', ctx.eps/2**10)
1098 verbose = option('verbose', False)
1099 maxterms = option('maxterms', ctx.dps*10)
1100 method = set(option('method', 'r+s').split('+'))
1101 skip = option('skip', 0)
1102 steps = iter(option('steps', xrange(10, 10**9, 10)))
1103 strict = option('strict')
1104 #steps = (10 for i in xrange(1000))
1105 summer=[]
1106 if 'd' in method or 'direct' in method:
1107 TRY_RICHARDSON = TRY_SHANKS = TRY_EULER_MACLAURIN = False
1108 else:
1109 TRY_RICHARDSON = ('r' in method) or ('richardson' in method)
1110 TRY_SHANKS = ('s' in method) or ('shanks' in method)
1111 TRY_EULER_MACLAURIN = ('e' in method) or \
1112 ('euler-maclaurin' in method)
1114 def init_levin(m):
1115 variant = kwargs.get("levin_variant", "u")
1116 if isinstance(variant, str):
1117 if variant == "all":
1118 variant = ["u", "v", "t"]
1119 else:
1120 variant = [variant]
1121 for s in variant:
1122 L = levin_class(method = m, variant = s)
1123 L.ctx = ctx
1124 L.name = m + "(" + s + ")"
1125 summer.append(L)
1127 if ('l' in method) or ('levin' in method):
1128 init_levin("levin")
1130 if ('sidi' in method):
1131 init_levin("sidi")
1133 if ('a' in method) or ('alternating' in method):
1134 L = cohen_alt_class()
1135 L.ctx = ctx
1136 L.name = "alternating"
1137 summer.append(L)
1139 last_richardson_value = 0
1140 shanks_table = []
1141 index = 0
1142 step = 10
1143 partial = []
1144 best = ctx.zero
1145 orig = ctx.prec
1146 try:
1147 if 'workprec' in kwargs:
1148 ctx.prec = kwargs['workprec']
1149 elif TRY_RICHARDSON or TRY_SHANKS or len(summer)!=0:
1150 ctx.prec = (ctx.prec+10) * 4
1151 else:
1152 ctx.prec += 30
1153 while 1:
1154 if index >= maxterms:
1155 break
1157 # Get new batch of terms
1158 try:
1159 step = next(steps)
1160 except StopIteration:
1161 pass
1162 if verbose:
1163 print("-"*70)
1164 print("Adding terms #%i-#%i" % (index, index+step))
1165 update(partial, xrange(index, index+step))
1166 index += step
1168 # Check direct error
1169 best = partial[-1]
1170 error = abs(best - partial[-2])
1171 if verbose:
1172 print("Direct error: %s" % ctx.nstr(error))
1173 if error <= tol:
1174 return best
1176 # Check each extrapolation method
1177 if TRY_RICHARDSON:
1178 value, maxc = ctx.richardson(partial)
1179 # Convergence
1180 richardson_error = abs(value - last_richardson_value)
1181 if verbose:
1182 print("Richardson error: %s" % ctx.nstr(richardson_error))
1183 # Convergence
1184 if richardson_error <= tol:
1185 return value
1186 last_richardson_value = value
1187 # Unreliable due to cancellation
1188 if ctx.eps*maxc > tol:
1189 if verbose:
1190 print("Ran out of precision for Richardson")
1191 TRY_RICHARDSON = False
1192 if richardson_error < error:
1193 error = richardson_error
1194 best = value
1195 if TRY_SHANKS:
1196 shanks_table = ctx.shanks(partial, shanks_table, randomized=True)
1197 row = shanks_table[-1]
1198 if len(row) == 2:
1199 est1 = row[-1]
1200 shanks_error = 0
1201 else:
1202 est1, maxc, est2 = row[-1], abs(row[-2]), row[-3]
1203 shanks_error = abs(est1-est2)
1204 if verbose:
1205 print("Shanks error: %s" % ctx.nstr(shanks_error))
1206 if shanks_error <= tol:
1207 return est1
1208 if ctx.eps*maxc > tol:
1209 if verbose:
1210 print("Ran out of precision for Shanks")
1211 TRY_SHANKS = False
1212 if shanks_error < error:
1213 error = shanks_error
1214 best = est1
1215 for L in summer:
1216 est, lerror = L.update_psum(partial)
1217 if verbose:
1218 print("%s error: %s" % (L.name, ctx.nstr(lerror)))
1219 if lerror <= tol:
1220 return est
1221 if lerror < error:
1222 error = lerror
1223 best = est
1224 if TRY_EULER_MACLAURIN:
1225 if ctx.mpc(ctx.sign(partial[-1]) / ctx.sign(partial[-2])).ae(-1):
1226 if verbose:
1227 print ("NOT using Euler-Maclaurin: the series appears"
1228 " to be alternating, so numerical\n quadrature"
1229 " will most likely fail")
1230 TRY_EULER_MACLAURIN = False
1231 else:
1232 value, em_error = emfun(index, tol)
1233 value += partial[-1]
1234 if verbose:
1235 print("Euler-Maclaurin error: %s" % ctx.nstr(em_error))
1236 if em_error <= tol:
1237 return value
1238 if em_error < error:
1239 best = value
1240 finally:
1241 ctx.prec = orig
1242 if strict:
1243 raise ctx.NoConvergence
1244 if verbose:
1245 print("Warning: failed to converge to target accuracy")
1246 return best
1248@defun
1249def nsum(ctx, f, *intervals, **options):
1250 r"""
1251 Computes the sum
1253 .. math :: S = \sum_{k=a}^b f(k)
1255 where `(a, b)` = *interval*, and where `a = -\infty` and/or
1256 `b = \infty` are allowed, or more generally
1258 .. math :: S = \sum_{k_1=a_1}^{b_1} \cdots
1259 \sum_{k_n=a_n}^{b_n} f(k_1,\ldots,k_n)
1261 if multiple intervals are given.
1263 Two examples of infinite series that can be summed by :func:`~mpmath.nsum`,
1264 where the first converges rapidly and the second converges slowly,
1265 are::
1267 >>> from mpmath import *
1268 >>> mp.dps = 15; mp.pretty = True
1269 >>> nsum(lambda n: 1/fac(n), [0, inf])
1270 2.71828182845905
1271 >>> nsum(lambda n: 1/n**2, [1, inf])
1272 1.64493406684823
1274 When appropriate, :func:`~mpmath.nsum` applies convergence acceleration to
1275 accurately estimate the sums of slowly convergent series. If the series is
1276 finite, :func:`~mpmath.nsum` currently does not attempt to perform any
1277 extrapolation, and simply calls :func:`~mpmath.fsum`.
1279 Multidimensional infinite series are reduced to a single-dimensional
1280 series over expanding hypercubes; if both infinite and finite dimensions
1281 are present, the finite ranges are moved innermost. For more advanced
1282 control over the summation order, use nested calls to :func:`~mpmath.nsum`,
1283 or manually rewrite the sum as a single-dimensional series.
1285 **Options**
1287 *tol*
1288 Desired maximum final error. Defaults roughly to the
1289 epsilon of the working precision.
1291 *method*
1292 Which summation algorithm to use (described below).
1293 Default: ``'richardson+shanks'``.
1295 *maxterms*
1296 Cancel after at most this many terms. Default: 10*dps.
1298 *steps*
1299 An iterable giving the number of terms to add between
1300 each extrapolation attempt. The default sequence is
1301 [10, 20, 30, 40, ...]. For example, if you know that
1302 approximately 100 terms will be required, efficiency might be
1303 improved by setting this to [100, 10]. Then the first
1304 extrapolation will be performed after 100 terms, the second
1305 after 110, etc.
1307 *verbose*
1308 Print details about progress.
1310 *ignore*
1311 If enabled, any term that raises ``ArithmeticError``
1312 or ``ValueError`` (e.g. through division by zero) is replaced
1313 by a zero. This is convenient for lattice sums with
1314 a singular term near the origin.
1316 **Methods**
1318 Unfortunately, an algorithm that can efficiently sum any infinite
1319 series does not exist. :func:`~mpmath.nsum` implements several different
1320 algorithms that each work well in different cases. The *method*
1321 keyword argument selects a method.
1323 The default method is ``'r+s'``, i.e. both Richardson extrapolation
1324 and Shanks transformation is attempted. A slower method that
1325 handles more cases is ``'r+s+e'``. For very high precision
1326 summation, or if the summation needs to be fast (for example if
1327 multiple sums need to be evaluated), it is a good idea to
1328 investigate which one method works best and only use that.
1330 ``'richardson'`` / ``'r'``:
1331 Uses Richardson extrapolation. Provides useful extrapolation
1332 when `f(k) \sim P(k)/Q(k)` or when `f(k) \sim (-1)^k P(k)/Q(k)`
1333 for polynomials `P` and `Q`. See :func:`~mpmath.richardson` for
1334 additional information.
1336 ``'shanks'`` / ``'s'``:
1337 Uses Shanks transformation. Typically provides useful
1338 extrapolation when `f(k) \sim c^k` or when successive terms
1339 alternate signs. Is able to sum some divergent series.
1340 See :func:`~mpmath.shanks` for additional information.
1342 ``'levin'`` / ``'l'``:
1343 Uses the Levin transformation. It performs better than the Shanks
1344 transformation for logarithmic convergent or alternating divergent
1345 series. The ``'levin_variant'``-keyword selects the variant. Valid
1346 choices are "u", "t", "v" and "all" whereby "all" uses all three
1347 u,t and v simultanously (This is good for performance comparison in
1348 conjunction with "verbose=True"). Instead of the Levin transform one can
1349 also use the Sidi-S transform by selecting the method ``'sidi'``.
1350 See :func:`~mpmath.levin` for additional details.
1352 ``'alternating'`` / ``'a'``:
1353 This is the convergence acceleration of alternating series developped
1354 by Cohen, Villegras and Zagier.
1355 See :func:`~mpmath.cohen_alt` for additional details.
1357 ``'euler-maclaurin'`` / ``'e'``:
1358 Uses the Euler-Maclaurin summation formula to approximate
1359 the remainder sum by an integral. This requires high-order
1360 numerical derivatives and numerical integration. The advantage
1361 of this algorithm is that it works regardless of the
1362 decay rate of `f`, as long as `f` is sufficiently smooth.
1363 See :func:`~mpmath.sumem` for additional information.
1365 ``'direct'`` / ``'d'``:
1366 Does not perform any extrapolation. This can be used
1367 (and should only be used for) rapidly convergent series.
1368 The summation automatically stops when the terms
1369 decrease below the target tolerance.
1371 **Basic examples**
1373 A finite sum::
1375 >>> nsum(lambda k: 1/k, [1, 6])
1376 2.45
1378 Summation of a series going to negative infinity and a doubly
1379 infinite series::
1381 >>> nsum(lambda k: 1/k**2, [-inf, -1])
1382 1.64493406684823
1383 >>> nsum(lambda k: 1/(1+k**2), [-inf, inf])
1384 3.15334809493716
1386 :func:`~mpmath.nsum` handles sums of complex numbers::
1388 >>> nsum(lambda k: (0.5+0.25j)**k, [0, inf])
1389 (1.6 + 0.8j)
1391 The following sum converges very rapidly, so it is most
1392 efficient to sum it by disabling convergence acceleration::
1394 >>> mp.dps = 1000
1395 >>> a = nsum(lambda k: -(-1)**k * k**2 / fac(2*k), [1, inf],
1396 ... method='direct')
1397 >>> b = (cos(1)+sin(1))/4
1398 >>> abs(a-b) < mpf('1e-998')
1399 True
1401 **Examples with Richardson extrapolation**
1403 Richardson extrapolation works well for sums over rational
1404 functions, as well as their alternating counterparts::
1406 >>> mp.dps = 50
1407 >>> nsum(lambda k: 1 / k**3, [1, inf],
1408 ... method='richardson')
1409 1.2020569031595942853997381615114499907649862923405
1410 >>> zeta(3)
1411 1.2020569031595942853997381615114499907649862923405
1413 >>> nsum(lambda n: (n + 3)/(n**3 + n**2), [1, inf],
1414 ... method='richardson')
1415 2.9348022005446793094172454999380755676568497036204
1416 >>> pi**2/2-2
1417 2.9348022005446793094172454999380755676568497036204
1419 >>> nsum(lambda k: (-1)**k / k**3, [1, inf],
1420 ... method='richardson')
1421 -0.90154267736969571404980362113358749307373971925537
1422 >>> -3*zeta(3)/4
1423 -0.90154267736969571404980362113358749307373971925538
1425 **Examples with Shanks transformation**
1427 The Shanks transformation works well for geometric series
1428 and typically provides excellent acceleration for Taylor
1429 series near the border of their disk of convergence.
1430 Here we apply it to a series for `\log(2)`, which can be
1431 seen as the Taylor series for `\log(1+x)` with `x = 1`::
1433 >>> nsum(lambda k: -(-1)**k/k, [1, inf],
1434 ... method='shanks')
1435 0.69314718055994530941723212145817656807550013436025
1436 >>> log(2)
1437 0.69314718055994530941723212145817656807550013436025
1439 Here we apply it to a slowly convergent geometric series::
1441 >>> nsum(lambda k: mpf('0.995')**k, [0, inf],
1442 ... method='shanks')
1443 200.0
1445 Finally, Shanks' method works very well for alternating series
1446 where `f(k) = (-1)^k g(k)`, and often does so regardless of
1447 the exact decay rate of `g(k)`::
1449 >>> mp.dps = 15
1450 >>> nsum(lambda k: (-1)**(k+1) / k**1.5, [1, inf],
1451 ... method='shanks')
1452 0.765147024625408
1453 >>> (2-sqrt(2))*zeta(1.5)/2
1454 0.765147024625408
1456 The following slowly convergent alternating series has no known
1457 closed-form value. Evaluating the sum a second time at higher
1458 precision indicates that the value is probably correct::
1460 >>> nsum(lambda k: (-1)**k / log(k), [2, inf],
1461 ... method='shanks')
1462 0.924299897222939
1463 >>> mp.dps = 30
1464 >>> nsum(lambda k: (-1)**k / log(k), [2, inf],
1465 ... method='shanks')
1466 0.92429989722293885595957018136
1468 **Examples with Levin transformation**
1470 The following example calculates Euler's constant as the constant term in
1471 the Laurent expansion of zeta(s) at s=1. This sum converges extremly slow
1472 because of the logarithmic convergence behaviour of the Dirichlet series
1473 for zeta.
1475 >>> mp.dps = 30
1476 >>> z = mp.mpf(10) ** (-10)
1477 >>> a = mp.nsum(lambda n: n**(-(1+z)), [1, mp.inf], method = "levin") - 1 / z
1478 >>> print(mp.chop(a - mp.euler, tol = 1e-10))
1479 0.0
1481 Now we sum the zeta function outside its range of convergence
1482 (attention: This does not work at the negative integers!):
1484 >>> mp.dps = 15
1485 >>> w = mp.nsum(lambda n: n ** (2 + 3j), [1, mp.inf], method = "levin", levin_variant = "v")
1486 >>> print(mp.chop(w - mp.zeta(-2-3j)))
1487 0.0
1489 The next example resummates an asymptotic series expansion of an integral
1490 related to the exponential integral.
1492 >>> mp.dps = 15
1493 >>> z = mp.mpf(10)
1494 >>> # exact = mp.quad(lambda x: mp.exp(-x)/(1+x/z),[0,mp.inf])
1495 >>> exact = z * mp.exp(z) * mp.expint(1,z) # this is the symbolic expression for the integral
1496 >>> w = mp.nsum(lambda n: (-1) ** n * mp.fac(n) * z ** (-n), [0, mp.inf], method = "sidi", levin_variant = "t")
1497 >>> print(mp.chop(w - exact))
1498 0.0
1500 Following highly divergent asymptotic expansion needs some care. Firstly we
1501 need copious amount of working precision. Secondly the stepsize must not be
1502 chosen to large, otherwise nsum may miss the point where the Levin transform
1503 converges and reach the point where only numerical garbage is produced due to
1504 numerical cancellation.
1506 >>> mp.dps = 15
1507 >>> z = mp.mpf(2)
1508 >>> # exact = mp.quad(lambda x: mp.exp( -x * x / 2 - z * x ** 4), [0,mp.inf]) * 2 / mp.sqrt(2 * mp.pi)
1509 >>> 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
1510 >>> w = mp.nsum(lambda n: (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n)),
1511 ... [0, mp.inf], method = "levin", levin_variant = "t", workprec = 8*mp.prec, steps = [2] + [1 for x in xrange(1000)])
1512 >>> print(mp.chop(w - exact))
1513 0.0
1515 The hypergeoemtric function can also be summed outside its range of convergence:
1517 >>> mp.dps = 15
1518 >>> z = 2 + 1j
1519 >>> exact = mp.hyp2f1(2 / mp.mpf(3), 4 / mp.mpf(3), 1 / mp.mpf(3), z)
1520 >>> 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))
1521 >>> v = mp.nsum(f, [0, mp.inf], method = "levin", steps = [10 for x in xrange(1000)])
1522 >>> print(mp.chop(exact-v))
1523 0.0
1525 **Examples with Cohen's alternating series resummation**
1527 The next example sums the alternating zeta function:
1529 >>> v = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "a")
1530 >>> print(mp.chop(v - mp.log(2)))
1531 0.0
1533 The derivate of the alternating zeta function outside its range of
1534 convergence:
1536 >>> v = mp.nsum(lambda n: (-1)**n * mp.log(n) * n, [1, mp.inf], method = "a")
1537 >>> print(mp.chop(v - mp.diff(lambda s: mp.altzeta(s), -1)))
1538 0.0
1540 **Examples with Euler-Maclaurin summation**
1542 The sum in the following example has the wrong rate of convergence
1543 for either Richardson or Shanks to be effective.
1545 >>> f = lambda k: log(k)/k**2.5
1546 >>> mp.dps = 15
1547 >>> nsum(f, [1, inf], method='euler-maclaurin')
1548 0.38734195032621
1549 >>> -diff(zeta, 2.5)
1550 0.38734195032621
1552 Increasing ``steps`` improves speed at higher precision::
1554 >>> mp.dps = 50
1555 >>> nsum(f, [1, inf], method='euler-maclaurin', steps=[250])
1556 0.38734195032620997271199237593105101319948228874688
1557 >>> -diff(zeta, 2.5)
1558 0.38734195032620997271199237593105101319948228874688
1560 **Divergent series**
1562 The Shanks transformation is able to sum some *divergent*
1563 series. In particular, it is often able to sum Taylor series
1564 beyond their radius of convergence (this is due to a relation
1565 between the Shanks transformation and Pade approximations;
1566 see :func:`~mpmath.pade` for an alternative way to evaluate divergent
1567 Taylor series). Furthermore the Levin-transform examples above
1568 contain some divergent series resummation.
1570 Here we apply it to `\log(1+x)` far outside the region of
1571 convergence::
1573 >>> mp.dps = 50
1574 >>> nsum(lambda k: -(-9)**k/k, [1, inf],
1575 ... method='shanks')
1576 2.3025850929940456840179914546843642076011014886288
1577 >>> log(10)
1578 2.3025850929940456840179914546843642076011014886288
1580 A particular type of divergent series that can be summed
1581 using the Shanks transformation is geometric series.
1582 The result is the same as using the closed-form formula
1583 for an infinite geometric series::
1585 >>> mp.dps = 15
1586 >>> for n in range(-8, 8):
1587 ... if n == 1:
1588 ... continue
1589 ... print("%s %s %s" % (mpf(n), mpf(1)/(1-n),
1590 ... nsum(lambda k: n**k, [0, inf], method='shanks')))
1591 ...
1592 -8.0 0.111111111111111 0.111111111111111
1593 -7.0 0.125 0.125
1594 -6.0 0.142857142857143 0.142857142857143
1595 -5.0 0.166666666666667 0.166666666666667
1596 -4.0 0.2 0.2
1597 -3.0 0.25 0.25
1598 -2.0 0.333333333333333 0.333333333333333
1599 -1.0 0.5 0.5
1600 0.0 1.0 1.0
1601 2.0 -1.0 -1.0
1602 3.0 -0.5 -0.5
1603 4.0 -0.333333333333333 -0.333333333333333
1604 5.0 -0.25 -0.25
1605 6.0 -0.2 -0.2
1606 7.0 -0.166666666666667 -0.166666666666667
1608 **Multidimensional sums**
1610 Any combination of finite and infinite ranges is allowed for the
1611 summation indices::
1613 >>> mp.dps = 15
1614 >>> nsum(lambda x,y: x+y, [2,3], [4,5])
1615 28.0
1616 >>> nsum(lambda x,y: x/2**y, [1,3], [1,inf])
1617 6.0
1618 >>> nsum(lambda x,y: y/2**x, [1,inf], [1,3])
1619 6.0
1620 >>> nsum(lambda x,y,z: z/(2**x*2**y), [1,inf], [1,inf], [3,4])
1621 7.0
1622 >>> nsum(lambda x,y,z: y/(2**x*2**z), [1,inf], [3,4], [1,inf])
1623 7.0
1624 >>> nsum(lambda x,y,z: x/(2**z*2**y), [3,4], [1,inf], [1,inf])
1625 7.0
1627 Some nice examples of double series with analytic solutions or
1628 reductions to single-dimensional series (see [1])::
1630 >>> nsum(lambda m, n: 1/2**(m*n), [1,inf], [1,inf])
1631 1.60669515241529
1632 >>> nsum(lambda n: 1/(2**n-1), [1,inf])
1633 1.60669515241529
1635 >>> nsum(lambda i,j: (-1)**(i+j)/(i**2+j**2), [1,inf], [1,inf])
1636 0.278070510848213
1637 >>> pi*(pi-3*ln2)/12
1638 0.278070510848213
1640 >>> nsum(lambda i,j: (-1)**(i+j)/(i+j)**2, [1,inf], [1,inf])
1641 0.129319852864168
1642 >>> altzeta(2) - altzeta(1)
1643 0.129319852864168
1645 >>> nsum(lambda i,j: (-1)**(i+j)/(i+j)**3, [1,inf], [1,inf])
1646 0.0790756439455825
1647 >>> altzeta(3) - altzeta(2)
1648 0.0790756439455825
1650 >>> nsum(lambda m,n: m**2*n/(3**m*(n*3**m+m*3**n)),
1651 ... [1,inf], [1,inf])
1652 0.28125
1653 >>> mpf(9)/32
1654 0.28125
1656 >>> nsum(lambda i,j: fac(i-1)*fac(j-1)/fac(i+j),
1657 ... [1,inf], [1,inf], workprec=400)
1658 1.64493406684823
1659 >>> zeta(2)
1660 1.64493406684823
1662 A hard example of a multidimensional sum is the Madelung constant
1663 in three dimensions (see [2]). The defining sum converges very
1664 slowly and only conditionally, so :func:`~mpmath.nsum` is lucky to
1665 obtain an accurate value through convergence acceleration. The
1666 second evaluation below uses a much more efficient, rapidly
1667 convergent 2D sum::
1669 >>> nsum(lambda x,y,z: (-1)**(x+y+z)/(x*x+y*y+z*z)**0.5,
1670 ... [-inf,inf], [-inf,inf], [-inf,inf], ignore=True)
1671 -1.74756459463318
1672 >>> nsum(lambda x,y: -12*pi*sech(0.5*pi * \
1673 ... sqrt((2*x+1)**2+(2*y+1)**2))**2, [0,inf], [0,inf])
1674 -1.74756459463318
1676 Another example of a lattice sum in 2D::
1678 >>> nsum(lambda x,y: (-1)**(x+y) / (x**2+y**2), [-inf,inf],
1679 ... [-inf,inf], ignore=True)
1680 -2.1775860903036
1681 >>> -pi*ln2
1682 -2.1775860903036
1684 An example of an Eisenstein series::
1686 >>> nsum(lambda m,n: (m+n*1j)**(-4), [-inf,inf], [-inf,inf],
1687 ... ignore=True)
1688 (3.1512120021539 + 0.0j)
1690 **References**
1692 1. [Weisstein]_ http://mathworld.wolfram.com/DoubleSeries.html,
1693 2. [Weisstein]_ http://mathworld.wolfram.com/MadelungConstants.html
1695 """
1696 infinite, g = standardize(ctx, f, intervals, options)
1697 if not infinite:
1698 return +g()
1700 def update(partial_sums, indices):
1701 if partial_sums:
1702 psum = partial_sums[-1]
1703 else:
1704 psum = ctx.zero
1705 for k in indices:
1706 psum = psum + g(ctx.mpf(k))
1707 partial_sums.append(psum)
1709 prec = ctx.prec
1711 def emfun(point, tol):
1712 workprec = ctx.prec
1713 ctx.prec = prec + 10
1714 v = ctx.sumem(g, [point, ctx.inf], tol, error=1)
1715 ctx.prec = workprec
1716 return v
1718 return +ctx.adaptive_extrapolation(update, emfun, options)
1721def wrapsafe(f):
1722 def g(*args):
1723 try:
1724 return f(*args)
1725 except (ArithmeticError, ValueError):
1726 return 0
1727 return g
1729def standardize(ctx, f, intervals, options):
1730 if options.get("ignore"):
1731 f = wrapsafe(f)
1732 finite = []
1733 infinite = []
1734 for k, points in enumerate(intervals):
1735 a, b = ctx._as_points(points)
1736 if b < a:
1737 return False, (lambda: ctx.zero)
1738 if a == ctx.ninf or b == ctx.inf:
1739 infinite.append((k, (a,b)))
1740 else:
1741 finite.append((k, (int(a), int(b))))
1742 if finite:
1743 f = fold_finite(ctx, f, finite)
1744 if not infinite:
1745 return False, lambda: f(*([0]*len(intervals)))
1746 if infinite:
1747 f = standardize_infinite(ctx, f, infinite)
1748 f = fold_infinite(ctx, f, infinite)
1749 args = [0] * len(intervals)
1750 d = infinite[0][0]
1751 def g(k):
1752 args[d] = k
1753 return f(*args)
1754 return True, g
1756# backwards compatible itertools.product
1757def cartesian_product(args):
1758 pools = map(tuple, args)
1759 result = [[]]
1760 for pool in pools:
1761 result = [x+[y] for x in result for y in pool]
1762 for prod in result:
1763 yield tuple(prod)
1765def fold_finite(ctx, f, intervals):
1766 if not intervals:
1767 return f
1768 indices = [v[0] for v in intervals]
1769 points = [v[1] for v in intervals]
1770 ranges = [xrange(a, b+1) for (a,b) in points]
1771 def g(*args):
1772 args = list(args)
1773 s = ctx.zero
1774 for xs in cartesian_product(ranges):
1775 for dim, x in zip(indices, xs):
1776 args[dim] = ctx.mpf(x)
1777 s += f(*args)
1778 return s
1779 #print "Folded finite", indices
1780 return g
1782# Standardize each interval to [0,inf]
1783def standardize_infinite(ctx, f, intervals):
1784 if not intervals:
1785 return f
1786 dim, [a,b] = intervals[-1]
1787 if a == ctx.ninf:
1788 if b == ctx.inf:
1789 def g(*args):
1790 args = list(args)
1791 k = args[dim]
1792 if k:
1793 s = f(*args)
1794 args[dim] = -k
1795 s += f(*args)
1796 return s
1797 else:
1798 return f(*args)
1799 else:
1800 def g(*args):
1801 args = list(args)
1802 args[dim] = b - args[dim]
1803 return f(*args)
1804 else:
1805 def g(*args):
1806 args = list(args)
1807 args[dim] += a
1808 return f(*args)
1809 #print "Standardized infinity along dimension", dim, a, b
1810 return standardize_infinite(ctx, g, intervals[:-1])
1812def fold_infinite(ctx, f, intervals):
1813 if len(intervals) < 2:
1814 return f
1815 dim1 = intervals[-2][0]
1816 dim2 = intervals[-1][0]
1817 # Assume intervals are [0,inf] x [0,inf] x ...
1818 def g(*args):
1819 args = list(args)
1820 #args.insert(dim2, None)
1821 n = int(args[dim1])
1822 s = ctx.zero
1823 #y = ctx.mpf(n)
1824 args[dim2] = ctx.mpf(n) #y
1825 for x in xrange(n+1):
1826 args[dim1] = ctx.mpf(x)
1827 s += f(*args)
1828 args[dim1] = ctx.mpf(n) #ctx.mpf(n)
1829 for y in xrange(n):
1830 args[dim2] = ctx.mpf(y)
1831 s += f(*args)
1832 return s
1833 #print "Folded infinite from", len(intervals), "to", (len(intervals)-1)
1834 return fold_infinite(ctx, g, intervals[:-1])
1836@defun
1837def nprod(ctx, f, interval, nsum=False, **kwargs):
1838 r"""
1839 Computes the product
1841 .. math ::
1843 P = \prod_{k=a}^b f(k)
1845 where `(a, b)` = *interval*, and where `a = -\infty` and/or
1846 `b = \infty` are allowed.
1848 By default, :func:`~mpmath.nprod` uses the same extrapolation methods as
1849 :func:`~mpmath.nsum`, except applied to the partial products rather than
1850 partial sums, and the same keyword options as for :func:`~mpmath.nsum` are
1851 supported. If ``nsum=True``, the product is instead computed via
1852 :func:`~mpmath.nsum` as
1854 .. math ::
1856 P = \exp\left( \sum_{k=a}^b \log(f(k)) \right).
1858 This is slower, but can sometimes yield better results. It is
1859 also required (and used automatically) when Euler-Maclaurin
1860 summation is requested.
1862 **Examples**
1864 A simple finite product::
1866 >>> from mpmath import *
1867 >>> mp.dps = 25; mp.pretty = True
1868 >>> nprod(lambda k: k, [1, 4])
1869 24.0
1871 A large number of infinite products have known exact values,
1872 and can therefore be used as a reference. Most of the following
1873 examples are taken from MathWorld [1].
1875 A few infinite products with simple values are::
1877 >>> 2*nprod(lambda k: (4*k**2)/(4*k**2-1), [1, inf])
1878 3.141592653589793238462643
1879 >>> nprod(lambda k: (1+1/k)**2/(1+2/k), [1, inf])
1880 2.0
1881 >>> nprod(lambda k: (k**3-1)/(k**3+1), [2, inf])
1882 0.6666666666666666666666667
1883 >>> nprod(lambda k: (1-1/k**2), [2, inf])
1884 0.5
1886 Next, several more infinite products with more complicated
1887 values::
1889 >>> nprod(lambda k: exp(1/k**2), [1, inf]); exp(pi**2/6)
1890 5.180668317897115748416626
1891 5.180668317897115748416626
1893 >>> nprod(lambda k: (k**2-1)/(k**2+1), [2, inf]); pi*csch(pi)
1894 0.2720290549821331629502366
1895 0.2720290549821331629502366
1897 >>> nprod(lambda k: (k**4-1)/(k**4+1), [2, inf])
1898 0.8480540493529003921296502
1899 >>> pi*sinh(pi)/(cosh(sqrt(2)*pi)-cos(sqrt(2)*pi))
1900 0.8480540493529003921296502
1902 >>> nprod(lambda k: (1+1/k+1/k**2)**2/(1+2/k+3/k**2), [1, inf])
1903 1.848936182858244485224927
1904 >>> 3*sqrt(2)*cosh(pi*sqrt(3)/2)**2*csch(pi*sqrt(2))/pi
1905 1.848936182858244485224927
1907 >>> nprod(lambda k: (1-1/k**4), [2, inf]); sinh(pi)/(4*pi)
1908 0.9190194775937444301739244
1909 0.9190194775937444301739244
1911 >>> nprod(lambda k: (1-1/k**6), [2, inf])
1912 0.9826842777421925183244759
1913 >>> (1+cosh(pi*sqrt(3)))/(12*pi**2)
1914 0.9826842777421925183244759
1916 >>> nprod(lambda k: (1+1/k**2), [2, inf]); sinh(pi)/(2*pi)
1917 1.838038955187488860347849
1918 1.838038955187488860347849
1920 >>> nprod(lambda n: (1+1/n)**n * exp(1/(2*n)-1), [1, inf])
1921 1.447255926890365298959138
1922 >>> exp(1+euler/2)/sqrt(2*pi)
1923 1.447255926890365298959138
1925 The following two products are equivalent and can be evaluated in
1926 terms of a Jacobi theta function. Pi can be replaced by any value
1927 (as long as convergence is preserved)::
1929 >>> nprod(lambda k: (1-pi**-k)/(1+pi**-k), [1, inf])
1930 0.3838451207481672404778686
1931 >>> nprod(lambda k: tanh(k*log(pi)/2), [1, inf])
1932 0.3838451207481672404778686
1933 >>> jtheta(4,0,1/pi)
1934 0.3838451207481672404778686
1936 This product does not have a known closed form value::
1938 >>> nprod(lambda k: (1-1/2**k), [1, inf])
1939 0.2887880950866024212788997
1941 A product taken from `-\infty`::
1943 >>> nprod(lambda k: 1-k**(-3), [-inf,-2])
1944 0.8093965973662901095786805
1945 >>> cosh(pi*sqrt(3)/2)/(3*pi)
1946 0.8093965973662901095786805
1948 A doubly infinite product::
1950 >>> nprod(lambda k: exp(1/(1+k**2)), [-inf, inf])
1951 23.41432688231864337420035
1952 >>> exp(pi/tanh(pi))
1953 23.41432688231864337420035
1955 A product requiring the use of Euler-Maclaurin summation to compute
1956 an accurate value::
1958 >>> nprod(lambda k: (1-1/k**2.5), [2, inf], method='e')
1959 0.696155111336231052898125
1961 **References**
1963 1. [Weisstein]_ http://mathworld.wolfram.com/InfiniteProduct.html
1965 """
1966 if nsum or ('e' in kwargs.get('method', '')):
1967 orig = ctx.prec
1968 try:
1969 # TODO: we are evaluating log(1+eps) -> eps, which is
1970 # inaccurate. This currently works because nsum greatly
1971 # increases the working precision. But we should be
1972 # more intelligent and handle the precision here.
1973 ctx.prec += 10
1974 v = ctx.nsum(lambda n: ctx.ln(f(n)), interval, **kwargs)
1975 finally:
1976 ctx.prec = orig
1977 return +ctx.exp(v)
1979 a, b = ctx._as_points(interval)
1980 if a == ctx.ninf:
1981 if b == ctx.inf:
1982 return f(0) * ctx.nprod(lambda k: f(-k) * f(k), [1, ctx.inf], **kwargs)
1983 return ctx.nprod(f, [-b, ctx.inf], **kwargs)
1984 elif b != ctx.inf:
1985 return ctx.fprod(f(ctx.mpf(k)) for k in xrange(int(a), int(b)+1))
1987 a = int(a)
1989 def update(partial_products, indices):
1990 if partial_products:
1991 pprod = partial_products[-1]
1992 else:
1993 pprod = ctx.one
1994 for k in indices:
1995 pprod = pprod * f(a + ctx.mpf(k))
1996 partial_products.append(pprod)
1998 return +ctx.adaptive_extrapolation(update, None, kwargs)
2001@defun
2002def limit(ctx, f, x, direction=1, exp=False, **kwargs):
2003 r"""
2004 Computes an estimate of the limit
2006 .. math ::
2008 \lim_{t \to x} f(t)
2010 where `x` may be finite or infinite.
2012 For finite `x`, :func:`~mpmath.limit` evaluates `f(x + d/n)` for
2013 consecutive integer values of `n`, where the approach direction
2014 `d` may be specified using the *direction* keyword argument.
2015 For infinite `x`, :func:`~mpmath.limit` evaluates values of
2016 `f(\mathrm{sign}(x) \cdot n)`.
2018 If the approach to the limit is not sufficiently fast to give
2019 an accurate estimate directly, :func:`~mpmath.limit` attempts to find
2020 the limit using Richardson extrapolation or the Shanks
2021 transformation. You can select between these methods using
2022 the *method* keyword (see documentation of :func:`~mpmath.nsum` for
2023 more information).
2025 **Options**
2027 The following options are available with essentially the
2028 same meaning as for :func:`~mpmath.nsum`: *tol*, *method*, *maxterms*,
2029 *steps*, *verbose*.
2031 If the option *exp=True* is set, `f` will be
2032 sampled at exponentially spaced points `n = 2^1, 2^2, 2^3, \ldots`
2033 instead of the linearly spaced points `n = 1, 2, 3, \ldots`.
2034 This can sometimes improve the rate of convergence so that
2035 :func:`~mpmath.limit` may return a more accurate answer (and faster).
2036 However, do note that this can only be used if `f`
2037 supports fast and accurate evaluation for arguments that
2038 are extremely close to the limit point (or if infinite,
2039 very large arguments).
2041 **Examples**
2043 A basic evaluation of a removable singularity::
2045 >>> from mpmath import *
2046 >>> mp.dps = 30; mp.pretty = True
2047 >>> limit(lambda x: (x-sin(x))/x**3, 0)
2048 0.166666666666666666666666666667
2050 Computing the exponential function using its limit definition::
2052 >>> limit(lambda n: (1+3/n)**n, inf)
2053 20.0855369231876677409285296546
2054 >>> exp(3)
2055 20.0855369231876677409285296546
2057 A limit for `\pi`::
2059 >>> f = lambda n: 2**(4*n+1)*fac(n)**4/(2*n+1)/fac(2*n)**2
2060 >>> limit(f, inf)
2061 3.14159265358979323846264338328
2063 Calculating the coefficient in Stirling's formula::
2065 >>> limit(lambda n: fac(n) / (sqrt(n)*(n/e)**n), inf)
2066 2.50662827463100050241576528481
2067 >>> sqrt(2*pi)
2068 2.50662827463100050241576528481
2070 Evaluating Euler's constant `\gamma` using the limit representation
2072 .. math ::
2074 \gamma = \lim_{n \rightarrow \infty } \left[ \left(
2075 \sum_{k=1}^n \frac{1}{k} \right) - \log(n) \right]
2077 (which converges notoriously slowly)::
2079 >>> f = lambda n: sum([mpf(1)/k for k in range(1,int(n)+1)]) - log(n)
2080 >>> limit(f, inf)
2081 0.577215664901532860606512090082
2082 >>> +euler
2083 0.577215664901532860606512090082
2085 With default settings, the following limit converges too slowly
2086 to be evaluated accurately. Changing to exponential sampling
2087 however gives a perfect result::
2089 >>> f = lambda x: sqrt(x**3+x**2)/(sqrt(x**3)+x)
2090 >>> limit(f, inf)
2091 0.992831158558330281129249686491
2092 >>> limit(f, inf, exp=True)
2093 1.0
2095 """
2097 if ctx.isinf(x):
2098 direction = ctx.sign(x)
2099 g = lambda k: f(ctx.mpf(k+1)*direction)
2100 else:
2101 direction *= ctx.one
2102 g = lambda k: f(x + direction/(k+1))
2103 if exp:
2104 h = g
2105 g = lambda k: h(2**k)
2107 def update(values, indices):
2108 for k in indices:
2109 values.append(g(k+1))
2111 # XXX: steps used by nsum don't work well
2112 if not 'steps' in kwargs:
2113 kwargs['steps'] = [10]
2115 return +ctx.adaptive_extrapolation(update, None, kwargs)