Coverage for /usr/lib/python3/dist-packages/sympy/concrete/summations.py: 8%
808 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from typing import Tuple as tTuple
3from sympy.calculus.singularities import is_decreasing
4from sympy.calculus.accumulationbounds import AccumulationBounds
5from .expr_with_intlimits import ExprWithIntLimits
6from .expr_with_limits import AddWithLimits
7from .gosper import gosper_sum
8from sympy.core.expr import Expr
9from sympy.core.add import Add
10from sympy.core.containers import Tuple
11from sympy.core.function import Derivative, expand
12from sympy.core.mul import Mul
13from sympy.core.numbers import Float, _illegal
14from sympy.core.relational import Eq
15from sympy.core.singleton import S
16from sympy.core.sorting import ordered
17from sympy.core.symbol import Dummy, Wild, Symbol, symbols
18from sympy.functions.combinatorial.factorials import factorial
19from sympy.functions.combinatorial.numbers import bernoulli, harmonic
20from sympy.functions.elementary.exponential import exp, log
21from sympy.functions.elementary.piecewise import Piecewise
22from sympy.functions.elementary.trigonometric import cot, csc
23from sympy.functions.special.hyper import hyper
24from sympy.functions.special.tensor_functions import KroneckerDelta
25from sympy.functions.special.zeta_functions import zeta
26from sympy.integrals.integrals import Integral
27from sympy.logic.boolalg import And
28from sympy.polys.partfrac import apart
29from sympy.polys.polyerrors import PolynomialError, PolificationFailed
30from sympy.polys.polytools import parallel_poly_from_expr, Poly, factor
31from sympy.polys.rationaltools import together
32from sympy.series.limitseq import limit_seq
33from sympy.series.order import O
34from sympy.series.residues import residue
35from sympy.sets.sets import FiniteSet, Interval
36from sympy.utilities.iterables import sift
37import itertools
40class Sum(AddWithLimits, ExprWithIntLimits):
41 r"""
42 Represents unevaluated summation.
44 Explanation
45 ===========
47 ``Sum`` represents a finite or infinite series, with the first argument
48 being the general form of terms in the series, and the second argument
49 being ``(dummy_variable, start, end)``, with ``dummy_variable`` taking
50 all integer values from ``start`` through ``end``. In accordance with
51 long-standing mathematical convention, the end term is included in the
52 summation.
54 Finite sums
55 ===========
57 For finite sums (and sums with symbolic limits assumed to be finite) we
58 follow the summation convention described by Karr [1], especially
59 definition 3 of section 1.4. The sum:
61 .. math::
63 \sum_{m \leq i < n} f(i)
65 has *the obvious meaning* for `m < n`, namely:
67 .. math::
69 \sum_{m \leq i < n} f(i) = f(m) + f(m+1) + \ldots + f(n-2) + f(n-1)
71 with the upper limit value `f(n)` excluded. The sum over an empty set is
72 zero if and only if `m = n`:
74 .. math::
76 \sum_{m \leq i < n} f(i) = 0 \quad \mathrm{for} \quad m = n
78 Finally, for all other sums over empty sets we assume the following
79 definition:
81 .. math::
83 \sum_{m \leq i < n} f(i) = - \sum_{n \leq i < m} f(i) \quad \mathrm{for} \quad m > n
85 It is important to note that Karr defines all sums with the upper
86 limit being exclusive. This is in contrast to the usual mathematical notation,
87 but does not affect the summation convention. Indeed we have:
89 .. math::
91 \sum_{m \leq i < n} f(i) = \sum_{i = m}^{n - 1} f(i)
93 where the difference in notation is intentional to emphasize the meaning,
94 with limits typeset on the top being inclusive.
96 Examples
97 ========
99 >>> from sympy.abc import i, k, m, n, x
100 >>> from sympy import Sum, factorial, oo, IndexedBase, Function
101 >>> Sum(k, (k, 1, m))
102 Sum(k, (k, 1, m))
103 >>> Sum(k, (k, 1, m)).doit()
104 m**2/2 + m/2
105 >>> Sum(k**2, (k, 1, m))
106 Sum(k**2, (k, 1, m))
107 >>> Sum(k**2, (k, 1, m)).doit()
108 m**3/3 + m**2/2 + m/6
109 >>> Sum(x**k, (k, 0, oo))
110 Sum(x**k, (k, 0, oo))
111 >>> Sum(x**k, (k, 0, oo)).doit()
112 Piecewise((1/(1 - x), Abs(x) < 1), (Sum(x**k, (k, 0, oo)), True))
113 >>> Sum(x**k/factorial(k), (k, 0, oo)).doit()
114 exp(x)
116 Here are examples to do summation with symbolic indices. You
117 can use either Function of IndexedBase classes:
119 >>> f = Function('f')
120 >>> Sum(f(n), (n, 0, 3)).doit()
121 f(0) + f(1) + f(2) + f(3)
122 >>> Sum(f(n), (n, 0, oo)).doit()
123 Sum(f(n), (n, 0, oo))
124 >>> f = IndexedBase('f')
125 >>> Sum(f[n]**2, (n, 0, 3)).doit()
126 f[0]**2 + f[1]**2 + f[2]**2 + f[3]**2
128 An example showing that the symbolic result of a summation is still
129 valid for seemingly nonsensical values of the limits. Then the Karr
130 convention allows us to give a perfectly valid interpretation to
131 those sums by interchanging the limits according to the above rules:
133 >>> S = Sum(i, (i, 1, n)).doit()
134 >>> S
135 n**2/2 + n/2
136 >>> S.subs(n, -4)
137 6
138 >>> Sum(i, (i, 1, -4)).doit()
139 6
140 >>> Sum(-i, (i, -3, 0)).doit()
141 6
143 An explicit example of the Karr summation convention:
145 >>> S1 = Sum(i**2, (i, m, m+n-1)).doit()
146 >>> S1
147 m**2*n + m*n**2 - m*n + n**3/3 - n**2/2 + n/6
148 >>> S2 = Sum(i**2, (i, m+n, m-1)).doit()
149 >>> S2
150 -m**2*n - m*n**2 + m*n - n**3/3 + n**2/2 - n/6
151 >>> S1 + S2
152 0
153 >>> S3 = Sum(i, (i, m, m-1)).doit()
154 >>> S3
155 0
157 See Also
158 ========
160 summation
161 Product, sympy.concrete.products.product
163 References
164 ==========
166 .. [1] Michael Karr, "Summation in Finite Terms", Journal of the ACM,
167 Volume 28 Issue 2, April 1981, Pages 305-350
168 https://dl.acm.org/doi/10.1145/322248.322255
169 .. [2] https://en.wikipedia.org/wiki/Summation#Capital-sigma_notation
170 .. [3] https://en.wikipedia.org/wiki/Empty_sum
171 """
173 __slots__ = ()
175 limits: tTuple[tTuple[Symbol, Expr, Expr]]
177 def __new__(cls, function, *symbols, **assumptions):
178 obj = AddWithLimits.__new__(cls, function, *symbols, **assumptions)
179 if not hasattr(obj, 'limits'):
180 return obj
181 if any(len(l) != 3 or None in l for l in obj.limits):
182 raise ValueError('Sum requires values for lower and upper bounds.')
184 return obj
186 def _eval_is_zero(self):
187 # a Sum is only zero if its function is zero or if all terms
188 # cancel out. This only answers whether the summand is zero; if
189 # not then None is returned since we don't analyze whether all
190 # terms cancel out.
191 if self.function.is_zero or self.has_empty_sequence:
192 return True
194 def _eval_is_extended_real(self):
195 if self.has_empty_sequence:
196 return True
197 return self.function.is_extended_real
199 def _eval_is_positive(self):
200 if self.has_finite_limits and self.has_reversed_limits is False:
201 return self.function.is_positive
203 def _eval_is_negative(self):
204 if self.has_finite_limits and self.has_reversed_limits is False:
205 return self.function.is_negative
207 def _eval_is_finite(self):
208 if self.has_finite_limits and self.function.is_finite:
209 return True
211 def doit(self, **hints):
212 if hints.get('deep', True):
213 f = self.function.doit(**hints)
214 else:
215 f = self.function
217 # first make sure any definite limits have summation
218 # variables with matching assumptions
219 reps = {}
220 for xab in self.limits:
221 d = _dummy_with_inherited_properties_concrete(xab)
222 if d:
223 reps[xab[0]] = d
224 if reps:
225 undo = {v: k for k, v in reps.items()}
226 did = self.xreplace(reps).doit(**hints)
227 if isinstance(did, tuple): # when separate=True
228 did = tuple([i.xreplace(undo) for i in did])
229 elif did is not None:
230 did = did.xreplace(undo)
231 else:
232 did = self
233 return did
236 if self.function.is_Matrix:
237 expanded = self.expand()
238 if self != expanded:
239 return expanded.doit()
240 return _eval_matrix_sum(self)
242 for n, limit in enumerate(self.limits):
243 i, a, b = limit
244 dif = b - a
245 if dif == -1:
246 # Any summation over an empty set is zero
247 return S.Zero
248 if dif.is_integer and dif.is_negative:
249 a, b = b + 1, a - 1
250 f = -f
252 newf = eval_sum(f, (i, a, b))
253 if newf is None:
254 if f == self.function:
255 zeta_function = self.eval_zeta_function(f, (i, a, b))
256 if zeta_function is not None:
257 return zeta_function
258 return self
259 else:
260 return self.func(f, *self.limits[n:])
261 f = newf
263 if hints.get('deep', True):
264 # eval_sum could return partially unevaluated
265 # result with Piecewise. In this case we won't
266 # doit() recursively.
267 if not isinstance(f, Piecewise):
268 return f.doit(**hints)
270 return f
272 def eval_zeta_function(self, f, limits):
273 """
274 Check whether the function matches with the zeta function.
276 If it matches, then return a `Piecewise` expression because
277 zeta function does not converge unless `s > 1` and `q > 0`
278 """
279 i, a, b = limits
280 w, y, z = Wild('w', exclude=[i]), Wild('y', exclude=[i]), Wild('z', exclude=[i])
281 result = f.match((w * i + y) ** (-z))
282 if result is not None and b is S.Infinity:
283 coeff = 1 / result[w] ** result[z]
284 s = result[z]
285 q = result[y] / result[w] + a
286 return Piecewise((coeff * zeta(s, q), And(q > 0, s > 1)), (self, True))
288 def _eval_derivative(self, x):
289 """
290 Differentiate wrt x as long as x is not in the free symbols of any of
291 the upper or lower limits.
293 Explanation
294 ===========
296 Sum(a*b*x, (x, 1, a)) can be differentiated wrt x or b but not `a`
297 since the value of the sum is discontinuous in `a`. In a case
298 involving a limit variable, the unevaluated derivative is returned.
299 """
301 # diff already confirmed that x is in the free symbols of self, but we
302 # don't want to differentiate wrt any free symbol in the upper or lower
303 # limits
304 # XXX remove this test for free_symbols when the default _eval_derivative is in
305 if isinstance(x, Symbol) and x not in self.free_symbols:
306 return S.Zero
308 # get limits and the function
309 f, limits = self.function, list(self.limits)
311 limit = limits.pop(-1)
313 if limits: # f is the argument to a Sum
314 f = self.func(f, *limits)
316 _, a, b = limit
317 if x in a.free_symbols or x in b.free_symbols:
318 return None
319 df = Derivative(f, x, evaluate=True)
320 rv = self.func(df, limit)
321 return rv
323 def _eval_difference_delta(self, n, step):
324 k, _, upper = self.args[-1]
325 new_upper = upper.subs(n, n + step)
327 if len(self.args) == 2:
328 f = self.args[0]
329 else:
330 f = self.func(*self.args[:-1])
332 return Sum(f, (k, upper + 1, new_upper)).doit()
334 def _eval_simplify(self, **kwargs):
336 function = self.function
338 if kwargs.get('deep', True):
339 function = function.simplify(**kwargs)
341 # split the function into adds
342 terms = Add.make_args(expand(function))
343 s_t = [] # Sum Terms
344 o_t = [] # Other Terms
346 for term in terms:
347 if term.has(Sum):
348 # if there is an embedded sum here
349 # it is of the form x * (Sum(whatever))
350 # hence we make a Mul out of it, and simplify all interior sum terms
351 subterms = Mul.make_args(expand(term))
352 out_terms = []
353 for subterm in subterms:
354 # go through each term
355 if isinstance(subterm, Sum):
356 # if it's a sum, simplify it
357 out_terms.append(subterm._eval_simplify(**kwargs))
358 else:
359 # otherwise, add it as is
360 out_terms.append(subterm)
362 # turn it back into a Mul
363 s_t.append(Mul(*out_terms))
364 else:
365 o_t.append(term)
367 # next try to combine any interior sums for further simplification
368 from sympy.simplify.simplify import factor_sum, sum_combine
369 result = Add(sum_combine(s_t), *o_t)
371 return factor_sum(result, limits=self.limits)
373 def is_convergent(self):
374 r"""
375 Checks for the convergence of a Sum.
377 Explanation
378 ===========
380 We divide the study of convergence of infinite sums and products in
381 two parts.
383 First Part:
384 One part is the question whether all the terms are well defined, i.e.,
385 they are finite in a sum and also non-zero in a product. Zero
386 is the analogy of (minus) infinity in products as
387 :math:`e^{-\infty} = 0`.
389 Second Part:
390 The second part is the question of convergence after infinities,
391 and zeros in products, have been omitted assuming that their number
392 is finite. This means that we only consider the tail of the sum or
393 product, starting from some point after which all terms are well
394 defined.
396 For example, in a sum of the form:
398 .. math::
400 \sum_{1 \leq i < \infty} \frac{1}{n^2 + an + b}
402 where a and b are numbers. The routine will return true, even if there
403 are infinities in the term sequence (at most two). An analogous
404 product would be:
406 .. math::
408 \prod_{1 \leq i < \infty} e^{\frac{1}{n^2 + an + b}}
410 This is how convergence is interpreted. It is concerned with what
411 happens at the limit. Finding the bad terms is another independent
412 matter.
414 Note: It is responsibility of user to see that the sum or product
415 is well defined.
417 There are various tests employed to check the convergence like
418 divergence test, root test, integral test, alternating series test,
419 comparison tests, Dirichlet tests. It returns true if Sum is convergent
420 and false if divergent and NotImplementedError if it cannot be checked.
422 References
423 ==========
425 .. [1] https://en.wikipedia.org/wiki/Convergence_tests
427 Examples
428 ========
430 >>> from sympy import factorial, S, Sum, Symbol, oo
431 >>> n = Symbol('n', integer=True)
432 >>> Sum(n/(n - 1), (n, 4, 7)).is_convergent()
433 True
434 >>> Sum(n/(2*n + 1), (n, 1, oo)).is_convergent()
435 False
436 >>> Sum(factorial(n)/5**n, (n, 1, oo)).is_convergent()
437 False
438 >>> Sum(1/n**(S(6)/5), (n, 1, oo)).is_convergent()
439 True
441 See Also
442 ========
444 Sum.is_absolutely_convergent
445 sympy.concrete.products.Product.is_convergent
446 """
447 p, q, r = symbols('p q r', cls=Wild)
449 sym = self.limits[0][0]
450 lower_limit = self.limits[0][1]
451 upper_limit = self.limits[0][2]
452 sequence_term = self.function.simplify()
454 if len(sequence_term.free_symbols) > 1:
455 raise NotImplementedError("convergence checking for more than one symbol "
456 "containing series is not handled")
458 if lower_limit.is_finite and upper_limit.is_finite:
459 return S.true
461 # transform sym -> -sym and swap the upper_limit = S.Infinity
462 # and lower_limit = - upper_limit
463 if lower_limit is S.NegativeInfinity:
464 if upper_limit is S.Infinity:
465 return Sum(sequence_term, (sym, 0, S.Infinity)).is_convergent() and \
466 Sum(sequence_term, (sym, S.NegativeInfinity, 0)).is_convergent()
467 from sympy.simplify.simplify import simplify
468 sequence_term = simplify(sequence_term.xreplace({sym: -sym}))
469 lower_limit = -upper_limit
470 upper_limit = S.Infinity
472 sym_ = Dummy(sym.name, integer=True, positive=True)
473 sequence_term = sequence_term.xreplace({sym: sym_})
474 sym = sym_
476 interval = Interval(lower_limit, upper_limit)
478 # Piecewise function handle
479 if sequence_term.is_Piecewise:
480 for func, cond in sequence_term.args:
481 # see if it represents something going to oo
482 if cond == True or cond.as_set().sup is S.Infinity:
483 s = Sum(func, (sym, lower_limit, upper_limit))
484 return s.is_convergent()
485 return S.true
487 ### -------- Divergence test ----------- ###
488 try:
489 lim_val = limit_seq(sequence_term, sym)
490 if lim_val is not None and lim_val.is_zero is False:
491 return S.false
492 except NotImplementedError:
493 pass
495 try:
496 lim_val_abs = limit_seq(abs(sequence_term), sym)
497 if lim_val_abs is not None and lim_val_abs.is_zero is False:
498 return S.false
499 except NotImplementedError:
500 pass
502 order = O(sequence_term, (sym, S.Infinity))
504 ### --------- p-series test (1/n**p) ---------- ###
505 p_series_test = order.expr.match(sym**p)
506 if p_series_test is not None:
507 if p_series_test[p] < -1:
508 return S.true
509 if p_series_test[p] >= -1:
510 return S.false
512 ### ------------- comparison test ------------- ###
513 # 1/(n**p*log(n)**q*log(log(n))**r) comparison
514 n_log_test = (order.expr.match(1/(sym**p*log(1/sym)**q*log(-log(1/sym))**r)) or
515 order.expr.match(1/(sym**p*(-log(1/sym))**q*log(-log(1/sym))**r)))
516 if n_log_test is not None:
517 if (n_log_test[p] > 1 or
518 (n_log_test[p] == 1 and n_log_test[q] > 1) or
519 (n_log_test[p] == n_log_test[q] == 1 and n_log_test[r] > 1)):
520 return S.true
521 return S.false
523 ### ------------- Limit comparison test -----------###
524 # (1/n) comparison
525 try:
526 lim_comp = limit_seq(sym*sequence_term, sym)
527 if lim_comp is not None and lim_comp.is_number and lim_comp > 0:
528 return S.false
529 except NotImplementedError:
530 pass
532 ### ----------- ratio test ---------------- ###
533 next_sequence_term = sequence_term.xreplace({sym: sym + 1})
534 from sympy.simplify.combsimp import combsimp
535 from sympy.simplify.powsimp import powsimp
536 ratio = combsimp(powsimp(next_sequence_term/sequence_term))
537 try:
538 lim_ratio = limit_seq(ratio, sym)
539 if lim_ratio is not None and lim_ratio.is_number:
540 if abs(lim_ratio) > 1:
541 return S.false
542 if abs(lim_ratio) < 1:
543 return S.true
544 except NotImplementedError:
545 lim_ratio = None
547 ### ---------- Raabe's test -------------- ###
548 if lim_ratio == 1: # ratio test inconclusive
549 test_val = sym*(sequence_term/
550 sequence_term.subs(sym, sym + 1) - 1)
551 test_val = test_val.gammasimp()
552 try:
553 lim_val = limit_seq(test_val, sym)
554 if lim_val is not None and lim_val.is_number:
555 if lim_val > 1:
556 return S.true
557 if lim_val < 1:
558 return S.false
559 except NotImplementedError:
560 pass
562 ### ----------- root test ---------------- ###
563 # lim = Limit(abs(sequence_term)**(1/sym), sym, S.Infinity)
564 try:
565 lim_evaluated = limit_seq(abs(sequence_term)**(1/sym), sym)
566 if lim_evaluated is not None and lim_evaluated.is_number:
567 if lim_evaluated < 1:
568 return S.true
569 if lim_evaluated > 1:
570 return S.false
571 except NotImplementedError:
572 pass
574 ### ------------- alternating series test ----------- ###
575 dict_val = sequence_term.match(S.NegativeOne**(sym + p)*q)
576 if not dict_val[p].has(sym) and is_decreasing(dict_val[q], interval):
577 return S.true
579 ### ------------- integral test -------------- ###
580 check_interval = None
581 from sympy.solvers.solveset import solveset
582 maxima = solveset(sequence_term.diff(sym), sym, interval)
583 if not maxima:
584 check_interval = interval
585 elif isinstance(maxima, FiniteSet) and maxima.sup.is_number:
586 check_interval = Interval(maxima.sup, interval.sup)
587 if (check_interval is not None and
588 (is_decreasing(sequence_term, check_interval) or
589 is_decreasing(-sequence_term, check_interval))):
590 integral_val = Integral(
591 sequence_term, (sym, lower_limit, upper_limit))
592 try:
593 integral_val_evaluated = integral_val.doit()
594 if integral_val_evaluated.is_number:
595 return S(integral_val_evaluated.is_finite)
596 except NotImplementedError:
597 pass
599 ### ----- Dirichlet and bounded times convergent tests ----- ###
600 # TODO
601 #
602 # Dirichlet_test
603 # https://en.wikipedia.org/wiki/Dirichlet%27s_test
604 #
605 # Bounded times convergent test
606 # It is based on comparison theorems for series.
607 # In particular, if the general term of a series can
608 # be written as a product of two terms a_n and b_n
609 # and if a_n is bounded and if Sum(b_n) is absolutely
610 # convergent, then the original series Sum(a_n * b_n)
611 # is absolutely convergent and so convergent.
612 #
613 # The following code can grows like 2**n where n is the
614 # number of args in order.expr
615 # Possibly combined with the potentially slow checks
616 # inside the loop, could make this test extremely slow
617 # for larger summation expressions.
619 if order.expr.is_Mul:
620 args = order.expr.args
621 argset = set(args)
623 ### -------------- Dirichlet tests -------------- ###
624 m = Dummy('m', integer=True)
625 def _dirichlet_test(g_n):
626 try:
627 ing_val = limit_seq(Sum(g_n, (sym, interval.inf, m)).doit(), m)
628 if ing_val is not None and ing_val.is_finite:
629 return S.true
630 except NotImplementedError:
631 pass
633 ### -------- bounded times convergent test ---------###
634 def _bounded_convergent_test(g1_n, g2_n):
635 try:
636 lim_val = limit_seq(g1_n, sym)
637 if lim_val is not None and (lim_val.is_finite or (
638 isinstance(lim_val, AccumulationBounds)
639 and (lim_val.max - lim_val.min).is_finite)):
640 if Sum(g2_n, (sym, lower_limit, upper_limit)).is_absolutely_convergent():
641 return S.true
642 except NotImplementedError:
643 pass
645 for n in range(1, len(argset)):
646 for a_tuple in itertools.combinations(args, n):
647 b_set = argset - set(a_tuple)
648 a_n = Mul(*a_tuple)
649 b_n = Mul(*b_set)
651 if is_decreasing(a_n, interval):
652 dirich = _dirichlet_test(b_n)
653 if dirich is not None:
654 return dirich
656 bc_test = _bounded_convergent_test(a_n, b_n)
657 if bc_test is not None:
658 return bc_test
660 _sym = self.limits[0][0]
661 sequence_term = sequence_term.xreplace({sym: _sym})
662 raise NotImplementedError("The algorithm to find the Sum convergence of %s "
663 "is not yet implemented" % (sequence_term))
665 def is_absolutely_convergent(self):
666 """
667 Checks for the absolute convergence of an infinite series.
669 Same as checking convergence of absolute value of sequence_term of
670 an infinite series.
672 References
673 ==========
675 .. [1] https://en.wikipedia.org/wiki/Absolute_convergence
677 Examples
678 ========
680 >>> from sympy import Sum, Symbol, oo
681 >>> n = Symbol('n', integer=True)
682 >>> Sum((-1)**n, (n, 1, oo)).is_absolutely_convergent()
683 False
684 >>> Sum((-1)**n/n**2, (n, 1, oo)).is_absolutely_convergent()
685 True
687 See Also
688 ========
690 Sum.is_convergent
691 """
692 return Sum(abs(self.function), self.limits).is_convergent()
694 def euler_maclaurin(self, m=0, n=0, eps=0, eval_integral=True):
695 """
696 Return an Euler-Maclaurin approximation of self, where m is the
697 number of leading terms to sum directly and n is the number of
698 terms in the tail.
700 With m = n = 0, this is simply the corresponding integral
701 plus a first-order endpoint correction.
703 Returns (s, e) where s is the Euler-Maclaurin approximation
704 and e is the estimated error (taken to be the magnitude of
705 the first omitted term in the tail):
707 >>> from sympy.abc import k, a, b
708 >>> from sympy import Sum
709 >>> Sum(1/k, (k, 2, 5)).doit().evalf()
710 1.28333333333333
711 >>> s, e = Sum(1/k, (k, 2, 5)).euler_maclaurin()
712 >>> s
713 -log(2) + 7/20 + log(5)
714 >>> from sympy import sstr
715 >>> print(sstr((s.evalf(), e.evalf()), full_prec=True))
716 (1.26629073187415, 0.0175000000000000)
718 The endpoints may be symbolic:
720 >>> s, e = Sum(1/k, (k, a, b)).euler_maclaurin()
721 >>> s
722 -log(a) + log(b) + 1/(2*b) + 1/(2*a)
723 >>> e
724 Abs(1/(12*b**2) - 1/(12*a**2))
726 If the function is a polynomial of degree at most 2n+1, the
727 Euler-Maclaurin formula becomes exact (and e = 0 is returned):
729 >>> Sum(k, (k, 2, b)).euler_maclaurin()
730 (b**2/2 + b/2 - 1, 0)
731 >>> Sum(k, (k, 2, b)).doit()
732 b**2/2 + b/2 - 1
734 With a nonzero eps specified, the summation is ended
735 as soon as the remainder term is less than the epsilon.
736 """
737 m = int(m)
738 n = int(n)
739 f = self.function
740 if len(self.limits) != 1:
741 raise ValueError("More than 1 limit")
742 i, a, b = self.limits[0]
743 if (a > b) == True:
744 if a - b == 1:
745 return S.Zero, S.Zero
746 a, b = b + 1, a - 1
747 f = -f
748 s = S.Zero
749 if m:
750 if b.is_Integer and a.is_Integer:
751 m = min(m, b - a + 1)
752 if not eps or f.is_polynomial(i):
753 s = Add(*[f.subs(i, a + k) for k in range(m)])
754 else:
755 term = f.subs(i, a)
756 if term:
757 test = abs(term.evalf(3)) < eps
758 if test == True:
759 return s, abs(term)
760 elif not (test == False):
761 # a symbolic Relational class, can't go further
762 return term, S.Zero
763 s = term
764 for k in range(1, m):
765 term = f.subs(i, a + k)
766 if abs(term.evalf(3)) < eps and term != 0:
767 return s, abs(term)
768 s += term
769 if b - a + 1 == m:
770 return s, S.Zero
771 a += m
772 x = Dummy('x')
773 I = Integral(f.subs(i, x), (x, a, b))
774 if eval_integral:
775 I = I.doit()
776 s += I
778 def fpoint(expr):
779 if b is S.Infinity:
780 return expr.subs(i, a), 0
781 return expr.subs(i, a), expr.subs(i, b)
782 fa, fb = fpoint(f)
783 iterm = (fa + fb)/2
784 g = f.diff(i)
785 for k in range(1, n + 2):
786 ga, gb = fpoint(g)
787 term = bernoulli(2*k)/factorial(2*k)*(gb - ga)
788 if k > n:
789 break
790 if eps and term:
791 term_evalf = term.evalf(3)
792 if term_evalf is S.NaN:
793 return S.NaN, S.NaN
794 if abs(term_evalf) < eps:
795 break
796 s += term
797 g = g.diff(i, 2, simplify=False)
798 return s + iterm, abs(term)
801 def reverse_order(self, *indices):
802 """
803 Reverse the order of a limit in a Sum.
805 Explanation
806 ===========
808 ``reverse_order(self, *indices)`` reverses some limits in the expression
809 ``self`` which can be either a ``Sum`` or a ``Product``. The selectors in
810 the argument ``indices`` specify some indices whose limits get reversed.
811 These selectors are either variable names or numerical indices counted
812 starting from the inner-most limit tuple.
814 Examples
815 ========
817 >>> from sympy import Sum
818 >>> from sympy.abc import x, y, a, b, c, d
820 >>> Sum(x, (x, 0, 3)).reverse_order(x)
821 Sum(-x, (x, 4, -1))
822 >>> Sum(x*y, (x, 1, 5), (y, 0, 6)).reverse_order(x, y)
823 Sum(x*y, (x, 6, 0), (y, 7, -1))
824 >>> Sum(x, (x, a, b)).reverse_order(x)
825 Sum(-x, (x, b + 1, a - 1))
826 >>> Sum(x, (x, a, b)).reverse_order(0)
827 Sum(-x, (x, b + 1, a - 1))
829 While one should prefer variable names when specifying which limits
830 to reverse, the index counting notation comes in handy in case there
831 are several symbols with the same name.
833 >>> S = Sum(x**2, (x, a, b), (x, c, d))
834 >>> S
835 Sum(x**2, (x, a, b), (x, c, d))
836 >>> S0 = S.reverse_order(0)
837 >>> S0
838 Sum(-x**2, (x, b + 1, a - 1), (x, c, d))
839 >>> S1 = S0.reverse_order(1)
840 >>> S1
841 Sum(x**2, (x, b + 1, a - 1), (x, d + 1, c - 1))
843 Of course we can mix both notations:
845 >>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(x, 1)
846 Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))
847 >>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(y, x)
848 Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))
850 See Also
851 ========
853 sympy.concrete.expr_with_intlimits.ExprWithIntLimits.index, reorder_limit,
854 sympy.concrete.expr_with_intlimits.ExprWithIntLimits.reorder
856 References
857 ==========
859 .. [1] Michael Karr, "Summation in Finite Terms", Journal of the ACM,
860 Volume 28 Issue 2, April 1981, Pages 305-350
861 https://dl.acm.org/doi/10.1145/322248.322255
862 """
863 l_indices = list(indices)
865 for i, indx in enumerate(l_indices):
866 if not isinstance(indx, int):
867 l_indices[i] = self.index(indx)
869 e = 1
870 limits = []
871 for i, limit in enumerate(self.limits):
872 l = limit
873 if i in l_indices:
874 e = -e
875 l = (limit[0], limit[2] + 1, limit[1] - 1)
876 limits.append(l)
878 return Sum(e * self.function, *limits)
880 def _eval_rewrite_as_Product(self, *args, **kwargs):
881 from sympy.concrete.products import Product
882 if self.function.is_extended_real:
883 return log(Product(exp(self.function), *self.limits))
886def summation(f, *symbols, **kwargs):
887 r"""
888 Compute the summation of f with respect to symbols.
890 Explanation
891 ===========
893 The notation for symbols is similar to the notation used in Integral.
894 summation(f, (i, a, b)) computes the sum of f with respect to i from a to b,
895 i.e.,
897 ::
899 b
900 ____
901 \ `
902 summation(f, (i, a, b)) = ) f
903 /___,
904 i = a
906 If it cannot compute the sum, it returns an unevaluated Sum object.
907 Repeated sums can be computed by introducing additional symbols tuples::
909 Examples
910 ========
912 >>> from sympy import summation, oo, symbols, log
913 >>> i, n, m = symbols('i n m', integer=True)
915 >>> summation(2*i - 1, (i, 1, n))
916 n**2
917 >>> summation(1/2**i, (i, 0, oo))
918 2
919 >>> summation(1/log(n)**n, (n, 2, oo))
920 Sum(log(n)**(-n), (n, 2, oo))
921 >>> summation(i, (i, 0, n), (n, 0, m))
922 m**3/6 + m**2/2 + m/3
924 >>> from sympy.abc import x
925 >>> from sympy import factorial
926 >>> summation(x**n/factorial(n), (n, 0, oo))
927 exp(x)
929 See Also
930 ========
932 Sum
933 Product, sympy.concrete.products.product
935 """
936 return Sum(f, *symbols, **kwargs).doit(deep=False)
939def telescopic_direct(L, R, n, limits):
940 """
941 Returns the direct summation of the terms of a telescopic sum
943 Explanation
944 ===========
946 L is the term with lower index
947 R is the term with higher index
948 n difference between the indexes of L and R
950 Examples
951 ========
953 >>> from sympy.concrete.summations import telescopic_direct
954 >>> from sympy.abc import k, a, b
955 >>> telescopic_direct(1/k, -1/(k+2), 2, (k, a, b))
956 -1/(b + 2) - 1/(b + 1) + 1/(a + 1) + 1/a
958 """
959 (i, a, b) = limits
960 return Add(*[L.subs(i, a + m) + R.subs(i, b - m) for m in range(n)])
963def telescopic(L, R, limits):
964 '''
965 Tries to perform the summation using the telescopic property.
967 Return None if not possible.
968 '''
969 (i, a, b) = limits
970 if L.is_Add or R.is_Add:
971 return None
973 # We want to solve(L.subs(i, i + m) + R, m)
974 # First we try a simple match since this does things that
975 # solve doesn't do, e.g. solve(cos(k+m)-cos(k), m) gives
976 # a more complicated solution than m == 0.
978 k = Wild("k")
979 sol = (-R).match(L.subs(i, i + k))
980 s = None
981 if sol and k in sol:
982 s = sol[k]
983 if not (s.is_Integer and L.subs(i, i + s) + R == 0):
984 # invalid match or match didn't work
985 s = None
987 # But there are things that match doesn't do that solve
988 # can do, e.g. determine that 1/(x + m) = 1/(1 - x) when m = 1
990 if s is None:
991 m = Dummy('m')
992 try:
993 from sympy.solvers.solvers import solve
994 sol = solve(L.subs(i, i + m) + R, m) or []
995 except NotImplementedError:
996 return None
997 sol = [si for si in sol if si.is_Integer and
998 (L.subs(i, i + si) + R).expand().is_zero]
999 if len(sol) != 1:
1000 return None
1001 s = sol[0]
1003 if s < 0:
1004 return telescopic_direct(R, L, abs(s), (i, a, b))
1005 elif s > 0:
1006 return telescopic_direct(L, R, s, (i, a, b))
1009def eval_sum(f, limits):
1010 (i, a, b) = limits
1011 if f.is_zero:
1012 return S.Zero
1013 if i not in f.free_symbols:
1014 return f*(b - a + 1)
1015 if a == b:
1016 return f.subs(i, a)
1017 if isinstance(f, Piecewise):
1018 if not any(i in arg.args[1].free_symbols for arg in f.args):
1019 # Piecewise conditions do not depend on the dummy summation variable,
1020 # therefore we can fold: Sum(Piecewise((e, c), ...), limits)
1021 # --> Piecewise((Sum(e, limits), c), ...)
1022 newargs = []
1023 for arg in f.args:
1024 newexpr = eval_sum(arg.expr, limits)
1025 if newexpr is None:
1026 return None
1027 newargs.append((newexpr, arg.cond))
1028 return f.func(*newargs)
1030 if f.has(KroneckerDelta):
1031 from .delta import deltasummation, _has_simple_delta
1032 f = f.replace(
1033 lambda x: isinstance(x, Sum),
1034 lambda x: x.factor()
1035 )
1036 if _has_simple_delta(f, limits[0]):
1037 return deltasummation(f, limits)
1039 dif = b - a
1040 definite = dif.is_Integer
1041 # Doing it directly may be faster if there are very few terms.
1042 if definite and (dif < 100):
1043 return eval_sum_direct(f, (i, a, b))
1044 if isinstance(f, Piecewise):
1045 return None
1046 # Try to do it symbolically. Even when the number of terms is
1047 # known, this can save time when b-a is big.
1048 value = eval_sum_symbolic(f.expand(), (i, a, b))
1049 if value is not None:
1050 return value
1051 # Do it directly
1052 if definite:
1053 return eval_sum_direct(f, (i, a, b))
1056def eval_sum_direct(expr, limits):
1057 """
1058 Evaluate expression directly, but perform some simple checks first
1059 to possibly result in a smaller expression and faster execution.
1060 """
1061 (i, a, b) = limits
1063 dif = b - a
1064 # Linearity
1065 if expr.is_Mul:
1066 # Try factor out everything not including i
1067 without_i, with_i = expr.as_independent(i)
1068 if without_i != 1:
1069 s = eval_sum_direct(with_i, (i, a, b))
1070 if s:
1071 r = without_i*s
1072 if r is not S.NaN:
1073 return r
1074 else:
1075 # Try term by term
1076 L, R = expr.as_two_terms()
1078 if not L.has(i):
1079 sR = eval_sum_direct(R, (i, a, b))
1080 if sR:
1081 return L*sR
1083 if not R.has(i):
1084 sL = eval_sum_direct(L, (i, a, b))
1085 if sL:
1086 return sL*R
1088 # do this whether its an Add or Mul
1089 # e.g. apart(1/(25*i**2 + 45*i + 14)) and
1090 # apart(1/((5*i + 2)*(5*i + 7))) ->
1091 # -1/(5*(5*i + 7)) + 1/(5*(5*i + 2))
1092 try:
1093 expr = apart(expr, i) # see if it becomes an Add
1094 except PolynomialError:
1095 pass
1097 if expr.is_Add:
1098 # Try factor out everything not including i
1099 without_i, with_i = expr.as_independent(i)
1100 if without_i != 0:
1101 s = eval_sum_direct(with_i, (i, a, b))
1102 if s:
1103 r = without_i*(dif + 1) + s
1104 if r is not S.NaN:
1105 return r
1106 else:
1107 # Try term by term
1108 L, R = expr.as_two_terms()
1109 lsum = eval_sum_direct(L, (i, a, b))
1110 rsum = eval_sum_direct(R, (i, a, b))
1112 if None not in (lsum, rsum):
1113 r = lsum + rsum
1114 if r is not S.NaN:
1115 return r
1117 return Add(*[expr.subs(i, a + j) for j in range(dif + 1)])
1120def eval_sum_symbolic(f, limits):
1121 f_orig = f
1122 (i, a, b) = limits
1123 if not f.has(i):
1124 return f*(b - a + 1)
1126 # Linearity
1127 if f.is_Mul:
1128 # Try factor out everything not including i
1129 without_i, with_i = f.as_independent(i)
1130 if without_i != 1:
1131 s = eval_sum_symbolic(with_i, (i, a, b))
1132 if s:
1133 r = without_i*s
1134 if r is not S.NaN:
1135 return r
1136 else:
1137 # Try term by term
1138 L, R = f.as_two_terms()
1140 if not L.has(i):
1141 sR = eval_sum_symbolic(R, (i, a, b))
1142 if sR:
1143 return L*sR
1145 if not R.has(i):
1146 sL = eval_sum_symbolic(L, (i, a, b))
1147 if sL:
1148 return sL*R
1150 # do this whether its an Add or Mul
1151 # e.g. apart(1/(25*i**2 + 45*i + 14)) and
1152 # apart(1/((5*i + 2)*(5*i + 7))) ->
1153 # -1/(5*(5*i + 7)) + 1/(5*(5*i + 2))
1154 try:
1155 f = apart(f, i)
1156 except PolynomialError:
1157 pass
1159 if f.is_Add:
1160 L, R = f.as_two_terms()
1161 lrsum = telescopic(L, R, (i, a, b))
1163 if lrsum:
1164 return lrsum
1166 # Try factor out everything not including i
1167 without_i, with_i = f.as_independent(i)
1168 if without_i != 0:
1169 s = eval_sum_symbolic(with_i, (i, a, b))
1170 if s:
1171 r = without_i*(b - a + 1) + s
1172 if r is not S.NaN:
1173 return r
1174 else:
1175 # Try term by term
1176 lsum = eval_sum_symbolic(L, (i, a, b))
1177 rsum = eval_sum_symbolic(R, (i, a, b))
1179 if None not in (lsum, rsum):
1180 r = lsum + rsum
1181 if r is not S.NaN:
1182 return r
1185 # Polynomial terms with Faulhaber's formula
1186 n = Wild('n')
1187 result = f.match(i**n)
1189 if result is not None:
1190 n = result[n]
1192 if n.is_Integer:
1193 if n >= 0:
1194 if (b is S.Infinity and a is not S.NegativeInfinity) or \
1195 (a is S.NegativeInfinity and b is not S.Infinity):
1196 return S.Infinity
1197 return ((bernoulli(n + 1, b + 1) - bernoulli(n + 1, a))/(n + 1)).expand()
1198 elif a.is_Integer and a >= 1:
1199 if n == -1:
1200 return harmonic(b) - harmonic(a - 1)
1201 else:
1202 return harmonic(b, abs(n)) - harmonic(a - 1, abs(n))
1204 if not (a.has(S.Infinity, S.NegativeInfinity) or
1205 b.has(S.Infinity, S.NegativeInfinity)):
1206 # Geometric terms
1207 c1 = Wild('c1', exclude=[i])
1208 c2 = Wild('c2', exclude=[i])
1209 c3 = Wild('c3', exclude=[i])
1210 wexp = Wild('wexp')
1212 # Here we first attempt powsimp on f for easier matching with the
1213 # exponential pattern, and attempt expansion on the exponent for easier
1214 # matching with the linear pattern.
1215 e = f.powsimp().match(c1 ** wexp)
1216 if e is not None:
1217 e_exp = e.pop(wexp).expand().match(c2*i + c3)
1218 if e_exp is not None:
1219 e.update(e_exp)
1221 p = (c1**c3).subs(e)
1222 q = (c1**c2).subs(e)
1223 r = p*(q**a - q**(b + 1))/(1 - q)
1224 l = p*(b - a + 1)
1225 return Piecewise((l, Eq(q, S.One)), (r, True))
1227 r = gosper_sum(f, (i, a, b))
1229 if isinstance(r, (Mul,Add)):
1230 from sympy.simplify.radsimp import denom
1231 from sympy.solvers.solvers import solve
1232 non_limit = r.free_symbols - Tuple(*limits[1:]).free_symbols
1233 den = denom(together(r))
1234 den_sym = non_limit & den.free_symbols
1235 args = []
1236 for v in ordered(den_sym):
1237 try:
1238 s = solve(den, v)
1239 m = Eq(v, s[0]) if s else S.false
1240 if m != False:
1241 args.append((Sum(f_orig.subs(*m.args), limits).doit(), m))
1242 break
1243 except NotImplementedError:
1244 continue
1246 args.append((r, True))
1247 return Piecewise(*args)
1249 if r not in (None, S.NaN):
1250 return r
1252 h = eval_sum_hyper(f_orig, (i, a, b))
1253 if h is not None:
1254 return h
1256 r = eval_sum_residue(f_orig, (i, a, b))
1257 if r is not None:
1258 return r
1260 factored = f_orig.factor()
1261 if factored != f_orig:
1262 return eval_sum_symbolic(factored, (i, a, b))
1265def _eval_sum_hyper(f, i, a):
1266 """ Returns (res, cond). Sums from a to oo. """
1267 if a != 0:
1268 return _eval_sum_hyper(f.subs(i, i + a), i, 0)
1270 if f.subs(i, 0) == 0:
1271 from sympy.simplify.simplify import simplify
1272 if simplify(f.subs(i, Dummy('i', integer=True, positive=True))) == 0:
1273 return S.Zero, True
1274 return _eval_sum_hyper(f.subs(i, i + 1), i, 0)
1276 from sympy.simplify.simplify import hypersimp
1277 hs = hypersimp(f, i)
1278 if hs is None:
1279 return None
1281 if isinstance(hs, Float):
1282 from sympy.simplify.simplify import nsimplify
1283 hs = nsimplify(hs)
1285 from sympy.simplify.combsimp import combsimp
1286 from sympy.simplify.hyperexpand import hyperexpand
1287 from sympy.simplify.radsimp import fraction
1288 numer, denom = fraction(factor(hs))
1289 top, topl = numer.as_coeff_mul(i)
1290 bot, botl = denom.as_coeff_mul(i)
1291 ab = [top, bot]
1292 factors = [topl, botl]
1293 params = [[], []]
1294 for k in range(2):
1295 for fac in factors[k]:
1296 mul = 1
1297 if fac.is_Pow:
1298 mul = fac.exp
1299 fac = fac.base
1300 if not mul.is_Integer:
1301 return None
1302 p = Poly(fac, i)
1303 if p.degree() != 1:
1304 return None
1305 m, n = p.all_coeffs()
1306 ab[k] *= m**mul
1307 params[k] += [n/m]*mul
1309 # Add "1" to numerator parameters, to account for implicit n! in
1310 # hypergeometric series.
1311 ap = params[0] + [1]
1312 bq = params[1]
1313 x = ab[0]/ab[1]
1314 h = hyper(ap, bq, x)
1315 f = combsimp(f)
1316 return f.subs(i, 0)*hyperexpand(h), h.convergence_statement
1319def eval_sum_hyper(f, i_a_b):
1320 i, a, b = i_a_b
1322 if f.is_hypergeometric(i) is False:
1323 return
1325 if (b - a).is_Integer:
1326 # We are never going to do better than doing the sum in the obvious way
1327 return None
1329 old_sum = Sum(f, (i, a, b))
1331 if b != S.Infinity:
1332 if a is S.NegativeInfinity:
1333 res = _eval_sum_hyper(f.subs(i, -i), i, -b)
1334 if res is not None:
1335 return Piecewise(res, (old_sum, True))
1336 else:
1337 n_illegal = lambda x: sum(x.count(_) for _ in _illegal)
1338 had = n_illegal(f)
1339 # check that no extra illegals are introduced
1340 res1 = _eval_sum_hyper(f, i, a)
1341 if res1 is None or n_illegal(res1) > had:
1342 return
1343 res2 = _eval_sum_hyper(f, i, b + 1)
1344 if res2 is None or n_illegal(res2) > had:
1345 return
1346 (res1, cond1), (res2, cond2) = res1, res2
1347 cond = And(cond1, cond2)
1348 if cond == False:
1349 return None
1350 return Piecewise((res1 - res2, cond), (old_sum, True))
1352 if a is S.NegativeInfinity:
1353 res1 = _eval_sum_hyper(f.subs(i, -i), i, 1)
1354 res2 = _eval_sum_hyper(f, i, 0)
1355 if res1 is None or res2 is None:
1356 return None
1357 res1, cond1 = res1
1358 res2, cond2 = res2
1359 cond = And(cond1, cond2)
1360 if cond == False or cond.as_set() == S.EmptySet:
1361 return None
1362 return Piecewise((res1 + res2, cond), (old_sum, True))
1364 # Now b == oo, a != -oo
1365 res = _eval_sum_hyper(f, i, a)
1366 if res is not None:
1367 r, c = res
1368 if c == False:
1369 if r.is_number:
1370 f = f.subs(i, Dummy('i', integer=True, positive=True) + a)
1371 if f.is_positive or f.is_zero:
1372 return S.Infinity
1373 elif f.is_negative:
1374 return S.NegativeInfinity
1375 return None
1376 return Piecewise(res, (old_sum, True))
1379def eval_sum_residue(f, i_a_b):
1380 r"""Compute the infinite summation with residues
1382 Notes
1383 =====
1385 If $f(n), g(n)$ are polynomials with $\deg(g(n)) - \deg(f(n)) \ge 2$,
1386 some infinite summations can be computed by the following residue
1387 evaluations.
1389 .. math::
1390 \sum_{n=-\infty, g(n) \ne 0}^{\infty} \frac{f(n)}{g(n)} =
1391 -\pi \sum_{\alpha|g(\alpha)=0}
1392 \text{Res}(\cot(\pi x) \frac{f(x)}{g(x)}, \alpha)
1394 .. math::
1395 \sum_{n=-\infty, g(n) \ne 0}^{\infty} (-1)^n \frac{f(n)}{g(n)} =
1396 -\pi \sum_{\alpha|g(\alpha)=0}
1397 \text{Res}(\csc(\pi x) \frac{f(x)}{g(x)}, \alpha)
1399 Examples
1400 ========
1402 >>> from sympy import Sum, oo, Symbol
1403 >>> x = Symbol('x')
1405 Doubly infinite series of rational functions.
1407 >>> Sum(1 / (x**2 + 1), (x, -oo, oo)).doit()
1408 pi/tanh(pi)
1410 Doubly infinite alternating series of rational functions.
1412 >>> Sum((-1)**x / (x**2 + 1), (x, -oo, oo)).doit()
1413 pi/sinh(pi)
1415 Infinite series of even rational functions.
1417 >>> Sum(1 / (x**2 + 1), (x, 0, oo)).doit()
1418 1/2 + pi/(2*tanh(pi))
1420 Infinite series of alternating even rational functions.
1422 >>> Sum((-1)**x / (x**2 + 1), (x, 0, oo)).doit()
1423 pi/(2*sinh(pi)) + 1/2
1425 This also have heuristics to transform arbitrarily shifted summand or
1426 arbitrarily shifted summation range to the canonical problem the
1427 formula can handle.
1429 >>> Sum(1 / (x**2 + 2*x + 2), (x, -1, oo)).doit()
1430 1/2 + pi/(2*tanh(pi))
1431 >>> Sum(1 / (x**2 + 4*x + 5), (x, -2, oo)).doit()
1432 1/2 + pi/(2*tanh(pi))
1433 >>> Sum(1 / (x**2 + 1), (x, 1, oo)).doit()
1434 -1/2 + pi/(2*tanh(pi))
1435 >>> Sum(1 / (x**2 + 1), (x, 2, oo)).doit()
1436 -1 + pi/(2*tanh(pi))
1438 References
1439 ==========
1441 .. [#] http://www.supermath.info/InfiniteSeriesandtheResidueTheorem.pdf
1443 .. [#] Asmar N.H., Grafakos L. (2018) Residue Theory.
1444 In: Complex Analysis with Applications.
1445 Undergraduate Texts in Mathematics. Springer, Cham.
1446 https://doi.org/10.1007/978-3-319-94063-2_5
1447 """
1448 i, a, b = i_a_b
1450 def is_even_function(numer, denom):
1451 """Test if the rational function is an even function"""
1452 numer_even = all(i % 2 == 0 for (i,) in numer.monoms())
1453 denom_even = all(i % 2 == 0 for (i,) in denom.monoms())
1454 numer_odd = all(i % 2 == 1 for (i,) in numer.monoms())
1455 denom_odd = all(i % 2 == 1 for (i,) in denom.monoms())
1456 return (numer_even and denom_even) or (numer_odd and denom_odd)
1458 def match_rational(f, i):
1459 numer, denom = f.as_numer_denom()
1460 try:
1461 (numer, denom), opt = parallel_poly_from_expr((numer, denom), i)
1462 except (PolificationFailed, PolynomialError):
1463 return None
1464 return numer, denom
1466 def get_poles(denom):
1467 roots = denom.sqf_part().all_roots()
1468 roots = sift(roots, lambda x: x.is_integer)
1469 if None in roots:
1470 return None
1471 int_roots, nonint_roots = roots[True], roots[False]
1472 return int_roots, nonint_roots
1474 def get_shift(denom):
1475 n = denom.degree(i)
1476 a = denom.coeff_monomial(i**n)
1477 b = denom.coeff_monomial(i**(n-1))
1478 shift = - b / a / n
1479 return shift
1481 #Need a dummy symbol with no assumptions set for get_residue_factor
1482 z = Dummy('z')
1484 def get_residue_factor(numer, denom, alternating):
1485 residue_factor = (numer.as_expr() / denom.as_expr()).subs(i, z)
1486 if not alternating:
1487 residue_factor *= cot(S.Pi * z)
1488 else:
1489 residue_factor *= csc(S.Pi * z)
1490 return residue_factor
1492 # We don't know how to deal with symbolic constants in summand
1493 if f.free_symbols - {i}:
1494 return None
1496 if not (a.is_Integer or a in (S.Infinity, S.NegativeInfinity)):
1497 return None
1498 if not (b.is_Integer or b in (S.Infinity, S.NegativeInfinity)):
1499 return None
1501 # Quick exit heuristic for the sums which doesn't have infinite range
1502 if a != S.NegativeInfinity and b != S.Infinity:
1503 return None
1505 match = match_rational(f, i)
1506 if match:
1507 alternating = False
1508 numer, denom = match
1509 else:
1510 match = match_rational(f / S.NegativeOne**i, i)
1511 if match:
1512 alternating = True
1513 numer, denom = match
1514 else:
1515 return None
1517 if denom.degree(i) - numer.degree(i) < 2:
1518 return None
1520 if (a, b) == (S.NegativeInfinity, S.Infinity):
1521 poles = get_poles(denom)
1522 if poles is None:
1523 return None
1524 int_roots, nonint_roots = poles
1526 if int_roots:
1527 return None
1529 residue_factor = get_residue_factor(numer, denom, alternating)
1530 residues = [residue(residue_factor, z, root) for root in nonint_roots]
1531 return -S.Pi * sum(residues)
1533 if not (a.is_finite and b is S.Infinity):
1534 return None
1536 if not is_even_function(numer, denom):
1537 # Try shifting summation and check if the summand can be made
1538 # and even function from the origin.
1539 # Sum(f(n), (n, a, b)) => Sum(f(n + s), (n, a - s, b - s))
1540 shift = get_shift(denom)
1542 if not shift.is_Integer:
1543 return None
1544 if shift == 0:
1545 return None
1547 numer = numer.shift(shift)
1548 denom = denom.shift(shift)
1550 if not is_even_function(numer, denom):
1551 return None
1553 if alternating:
1554 f = S.NegativeOne**i * (S.NegativeOne**shift * numer.as_expr() / denom.as_expr())
1555 else:
1556 f = numer.as_expr() / denom.as_expr()
1557 return eval_sum_residue(f, (i, a-shift, b-shift))
1559 poles = get_poles(denom)
1560 if poles is None:
1561 return None
1562 int_roots, nonint_roots = poles
1564 if int_roots:
1565 int_roots = [int(root) for root in int_roots]
1566 int_roots_max = max(int_roots)
1567 int_roots_min = min(int_roots)
1568 # Integer valued poles must be next to each other
1569 # and also symmetric from origin (Because the function is even)
1570 if not len(int_roots) == int_roots_max - int_roots_min + 1:
1571 return None
1573 # Check whether the summation indices contain poles
1574 if a <= max(int_roots):
1575 return None
1577 residue_factor = get_residue_factor(numer, denom, alternating)
1578 residues = [residue(residue_factor, z, root) for root in int_roots + nonint_roots]
1579 full_sum = -S.Pi * sum(residues)
1581 if not int_roots:
1582 # Compute Sum(f, (i, 0, oo)) by adding a extraneous evaluation
1583 # at the origin.
1584 half_sum = (full_sum + f.xreplace({i: 0})) / 2
1586 # Add and subtract extraneous evaluations
1587 extraneous_neg = [f.xreplace({i: i0}) for i0 in range(int(a), 0)]
1588 extraneous_pos = [f.xreplace({i: i0}) for i0 in range(0, int(a))]
1589 result = half_sum + sum(extraneous_neg) - sum(extraneous_pos)
1591 return result
1593 # Compute Sum(f, (i, min(poles) + 1, oo))
1594 half_sum = full_sum / 2
1596 # Subtract extraneous evaluations
1597 extraneous = [f.xreplace({i: i0}) for i0 in range(max(int_roots) + 1, int(a))]
1598 result = half_sum - sum(extraneous)
1600 return result
1603def _eval_matrix_sum(expression):
1604 f = expression.function
1605 for n, limit in enumerate(expression.limits):
1606 i, a, b = limit
1607 dif = b - a
1608 if dif.is_Integer:
1609 if (dif < 0) == True:
1610 a, b = b + 1, a - 1
1611 f = -f
1613 newf = eval_sum_direct(f, (i, a, b))
1614 if newf is not None:
1615 return newf.doit()
1618def _dummy_with_inherited_properties_concrete(limits):
1619 """
1620 Return a Dummy symbol that inherits as many assumptions as possible
1621 from the provided symbol and limits.
1623 If the symbol already has all True assumption shared by the limits
1624 then return None.
1625 """
1626 x, a, b = limits
1627 l = [a, b]
1629 assumptions_to_consider = ['extended_nonnegative', 'nonnegative',
1630 'extended_nonpositive', 'nonpositive',
1631 'extended_positive', 'positive',
1632 'extended_negative', 'negative',
1633 'integer', 'rational', 'finite',
1634 'zero', 'real', 'extended_real']
1636 assumptions_to_keep = {}
1637 assumptions_to_add = {}
1638 for assum in assumptions_to_consider:
1639 assum_true = x._assumptions.get(assum, None)
1640 if assum_true:
1641 assumptions_to_keep[assum] = True
1642 elif all(getattr(i, 'is_' + assum) for i in l):
1643 assumptions_to_add[assum] = True
1644 if assumptions_to_add:
1645 assumptions_to_keep.update(assumptions_to_add)
1646 return Dummy('d', **assumptions_to_keep)