sympy/concrete/summations.py

   1
   2 from sympy.core import (Basic, S, C, Add, Mul, Symbol, Equality, Interval,
   3     sympify, symbols, Wild)
   4 from sympy.functions import factorial
   5 from sympy.solvers import solve
   6
   7 class Sum(Basic):
   8     """Represents unevaluated summation."""
   9
  10     def __new__(cls, f, *symbols, **assumptions):
  11         f = sympify(f)
  12
  13         if f.is_Number:
  14             if f is S.NaN:
  15                 return S.NaN
  16             elif f is S.Zero:
  17                 return S.Zero
  18
  19         if not symbols:
  20             limits = f.atoms(Symbol)
  21
  22             if not limits:
  23                 return f
  24         else:
  25             limits = []
  26
  27             for V in symbols:
  28                 if isinstance(V, Symbol):
  29                     limits.append(V)
  30                     continue
  31                 elif isinstance(V, Equality):
  32                     if isinstance(V.lhs, Symbol):
  33                         if isinstance(V.rhs, Interval):
  34                             limits.append((V.lhs, V.rhs.start, V.rhs.end))
  35                         else:
  36                             limits.append((V.lhs, V.rhs))
  37
  38                         continue
  39                 elif isinstance(V, (tuple, list)):
  40                     if len(V) == 1:
  41                         if isinstance(V[0], Symbol):
  42                             limits.append(V[0])
  43                             continue
  44                     elif len(V) in (2, 3):
  45                         if isinstance(V[0], Symbol):
  46                             limits.append(tuple(map(sympify, V)))
  47                             continue
  48
  49                 raise ValueError("Invalid summation variable or limits")
  50
  51         obj = Basic.__new__(cls, **assumptions)
  52         obj._args = (f, tuple(limits))
  53
  54         return obj
  55
  56     @property
  57     def function(self):
  58         return self._args[0]
  59
  60     @property
  61     def limits(self):
  62         return self._args[1]
  63
  64     def doit(self, **hints):
  65         #if not hints.get('sums', True):
  66         #    return self
  67         f = self.function
  68         for i, a, b in self.limits:
  69             f = eval_sum(f, (i, a, b))
  70             if f is None:
  71                 return self
  72         return f
  73
  74     def _eval_summation(self, f, x):
  75         return
  76
  77     def euler_maclaurin(self, m=0, n=0, eps=0, eval_integral=True):
  78         """
  79         Return an Euler-Maclaurin approximation of self, where m is the
  80         number of leading terms to sum directly and n is the number of
  81         terms in the tail.
  82
  83         With m = n = 0, this is simply the corresponding integral
  84         plus a first-order endpoint correction.
  85
  86         Returns (s, e) where s is the Euler-Maclaurin approximation
  87         and e is the estimated error (taken to be the magnitude of
  88         the first omitted term in the tail):
  89
  90             >>> k = Symbol('k')
  91             >>> Sum(1/k, (k, 2, 5)).doit().evalf()
  92             1.28333333333333
  93             >>> s, e = Sum(1/k, (k, 2, 5)).euler_maclaurin()
  94             >>> s
  95             7/20 - log(2) + log(5)
  96             >>> s.evalf(), e.evalf()
  97             (1.26629073187416, 0.0175000000000000)
  98
  99         The endpoints may be symbolic:
 100
 101             >>> k, a, b = symbols('kab')
 102             >>> s, e = Sum(1/k, (k, a, b)).euler_maclaurin()
 103             >>> s
 104             -log(a) + log(b) + 1/(2*a) + 1/(2*b)
 105             >>> e
 106             abs(-1/(12*b**2) + 1/(12*a**2))
 107
 108         If the function is a polynomial of degree at most 2n+1, the
 109         Euler-Maclaurin formula becomes exact (and e = 0 is returned):
 110
 111             >>> Sum(k, (k, 2, b)).euler_maclaurin()
 112             (-1 + b/2 + 1/2*b**2, 0)
 113             >>> Sum(k, (k, 2, b)).doit()
 114             -1 + b/2 + 1/2*b**2
 115
 116         With a nonzero eps specified, the summation is ended
 117         as soon as the remainder term is less than the epsilon.
 118         """
 119         m = int(m)
 120         n = int(n)
 121         f = self.function
 122         assert len(self.limits) == 1
 123         i, a, b = self.limits[0]
 124         s = S.Zero
 125         if m:
 126             for k in range(m):
 127                 term = f.subs(i, a+k)
 128                 if (eps and term and abs(term.evalf(3)) < eps):
 129                     return s, abs(term)
 130                 s += term
 131             a += m
 132         x = Symbol('x', dummy=True)
 133         I = C.Integral(f.subs(i, x), (x, a, b))
 134         if eval_integral:
 135             I = I.doit()
 136         s += I
 137         def fpoint(expr):
 138             if b is S.Infinity:
 139                 return expr.subs(i, a), 0
 140             return expr.subs(i, a), expr.subs(i, b)
 141         fa, fb = fpoint(f)
 142         iterm = (fa + fb)/2
 143         g = f.diff(i)
 144         for k in xrange(1, n+2):
 145             ga, gb = fpoint(g)
 146             term = C.bernoulli(2*k)/C.Factorial(2*k)*(gb-ga)
 147             if (eps and term and abs(term.evalf(3)) < eps) or (k > n):
 148                 break
 149             s += term
 150             g = g.diff(i, 2)
 151         return s + iterm, abs(term)
 152
 153
 154 def sum(*args, **kwargs):
 155     summation = Sum(*args, **kwargs)
 156
 157     if isinstance(summation, Sum):
 158         return summation.doit()
 159     else:
 160         return summation
 161
 162
 163 def getab(expr):
 164     cls = expr.__class__
 165     return cls(expr.args[0]), cls(*expr.args[1:])
 166
 167 def telescopic_direct(L, R, n, (i, a, b)):
 168     '''Returns the direct summation of the terms of a telescopic sum
 169
 170     L is the term with lower index
 171     R is the term with higher index
 172     n difference between the indexes of L and R
 173
 174     For example:
 175
 176     >>> k,a,b = symbols('kab')
 177     >>> telescopic_direct(1/k, -1/(k+2), 2, (k, a, b))
 178     1/a + 1/(1 + a) - 1/(1 + b) - 1/(2 + b)
 179     '''
 180     s = 0
 181     for m in xrange(n):
 182         s += L.subs(i,a+m) + R.subs(i,b-m)
 183     return s
 184
 185 def telescopic(L, R, (i, a, b)):
 186     '''Tries to perform the summation using the telescopic property
 187
 188     return None if not possible
 189     '''
 190     if L.is_Add or R.is_Add:
 191         return None
 192     s = None
 193     #First we try to solve using match
 194     #Maybe this should go inside solve
 195     k = Wild("k")
 196     sol = (-R).match(L.subs(i, i + k))
 197     if sol and sol.has_key(k):
 198         if L.subs(i,i + sol[k]) == -R:
 199             #sometimes match fail(f(x+2).match(-f(x+k))->{k: -2 - 2x}))
 200             s = sol[k]
 201     #Then we try to solve using solve
 202     if not s or not s.is_Integer:
 203         m = Symbol("m")
 204         try:
 205             s = solve(L.subs(i, i + m) + R, m)[0]
 206         except ValueError:
 207             pass
 208     if s and s.is_Integer:
 209         if s < 0:
 210             return telescopic_direct(R, L, abs(s), (i, a, b))
 211         elif s > 0:
 212             return telescopic_direct(L, R, s, (i, a, b))
 213     return None
 214
 215 def eval_sum(f, (i, a, b)):
 216     if not f.has(i):
 217         return f*(b-a+1)
 218     definite = a.is_Integer and b.is_Integer
 219     # Doing it directly may be faster if there are very few terms.
 220     if definite and (b-a < 100):
 221         return eval_sum_direct(f, (i, a, b))
 222     # Try to do it symbolically. Even when the number of terms is known,
 223     # this can save time when b-a is big.
 224     # We should try to transform to partial fractions
 225     value = eval_sum_symbolic(f.expand(), (i, a, b))
 226     if value is not None:
 227         return value
 228     # Do it directly
 229     if definite:
 230         return eval_sum_direct(f, (i, a, b))
 231
 232 def eval_sum_symbolic(f, (i, a, b)):
 233     if not f.has(i):
 234         return f*(b-a+1)
 235     # Linearity
 236     if f.is_Mul:
 237         L, R = getab(f)
 238         if not L.has(i):
 239             sR = eval_sum_symbolic(R, (i, a, b))
 240             if sR: return L*sR
 241         if not R.has(i):
 242             sL = eval_sum_symbolic(L, (i, a, b))
 243             if sL: return R*sL
 244     if f.is_Add:
 245         L, R = getab(f)
 246         lrsum = telescopic(L, R, (i, a, b))
 247         if lrsum: return lrsum
 248         lsum = eval_sum_symbolic(L, (i, a, b))
 249         rsum = eval_sum_symbolic(R, (i, a, b))
 250         if None not in (lsum, rsum):
 251             return lsum + rsum
 252     # Polynomial terms with Faulhaber's formula
 253     p = C.Wild('p')
 254     e = f.match(i**p)
 255     if e != None:
 256         c = p.subs(e)
 257         B = C.bernoulli
 258         if c.is_integer and c >= 0:
 259             s = (B(c+1, b+1) - B(c+1, a))/(c+1)
 260             return s.expand()
 261     # Geometric terms
 262     c1 = C.Wild('c1', exclude=[i])
 263     c2 = C.Wild('c2', exclude=[i])
 264     c3 = C.Wild('c3', exclude=[i])
 265     e = f.match(c1**(c2*i+c3))
 266     if e is not None:
 267         c1 = c1.subs(e)
 268         c2 = c2.subs(e)
 269         c3 = c3.subs(e)
 270         # TODO: more general limit handling
 271         return c1**c3 * (c1**(a*c2) - c1**(c2+b*c2)) / (1 - c1**c2)
 272     return None
 273
 274 def eval_sum_direct(expr, (i, a, b)):
 275     s = S.Zero
 276     if expr.has(i):
 277         for j in xrange(a, b+1):
 278             s += expr.subs(i, j)
 279     else:
 280         for j in xrange(a, b+1):
 281             s += expr
 282     return s