sympy/series/limits_series.py

   1 # This is a very simplified (and not that reliable) version of the limit
   2 # algorithm that is only used for series expansion. Use this file only in the
   3 # series expansion code. Everywere else, use the general and robust limit
   4 # algorithm in limits.py.
   5
   6 # The idea is to use the general and robust algorithm in limits.py even in the
   7 # series expansion, but currently the limits.py is slower and there were some
   8 # other bugs (mainly recursion), if it were used in the series expansion. So
   9 # currently we use limits_series, until we move to limits.py completely.
  10
  11 from sympy.core.basic import Basic, S
  12 from sympy.core.methods import RelMeths, ArithMeths
  13 from sympy.core.cache import cache_it, cache_it_immutable
  14
  15 class Limit(Basic, RelMeths, ArithMeths):
  16     """ Find the limit of the expression under process x->xlim.
  17
  18     Limit(expr, x, xlim)
  19     """
  20     @cache_it_immutable
  21     def __new__(cls, expr, x, xlim, direction='<', **assumptions):
  22         expr = Basic.sympify(expr)
  23         x = Basic.sympify(x)
  24         xlim = Basic.sympify(xlim)
  25         if not isinstance(x, Basic.Symbol):
  26             raise ValueError("Limit 2nd argument must be Symbol instance (got %s)" % (x))
  27         assert isinstance(x, Basic.Symbol),`x`
  28
  29         if not expr.has(x):
  30             return expr
  31
  32         if isinstance(xlim, Basic.NegativeInfinity):
  33             xoo = InfLimit.limit_process_symbol()
  34             if expr.has(xoo):
  35                 xoo = Basic.Symbol(x.name + '_oo',dummy=True,positive=True,unbounded=True)
  36             return InfLimit(expr.subs(x,-xoo), xoo)
  37         if isinstance(xlim, Basic.Infinity):
  38             return InfLimit(expr, x)
  39         else:
  40             xoo = InfLimit.limit_process_symbol()
  41             if expr.has(xoo):
  42                 xoo = Basic.Symbol(x.name + '_oo',dummy=True,positive=True,unbounded=True)
  43             if direction=='<':
  44                 return InfLimit(expr.subs(x, xlim+1/xoo), xoo)
  45             elif direction=='>':
  46                 return InfLimit(expr.subs(x, xlim-1/xoo), xoo)
  47             else:
  48                 raise ValueError("Limit direction must be < or > (got %s)" % (direction))
  49
  50         # XXX This code is currently unreachable
  51         obj = Basic.__new__(cls, expr, x, xlim, **assumptions)
  52         obj.direction = direction
  53         return obj
  54
  55     def _hashable_content(self):
  56         return self._args + (self.direction,)
  57
  58     @property
  59     def expr(self):
  60         return self._args[0]
  61
  62     @property
  63     def var(self):
  64         return self._args[1]
  65
  66     @property
  67     def varlim(self):
  68         return self._args[2]
  69
  70     def tostr(self, level=0):
  71         if isinstance(self.varlim,(Basic.Infinity, Basic.NegativeInfinity)): s = ''
  72         elif self.direction=='<': s = '+'
  73         else: s = '-'
  74         r = 'lim[%s->%s%s](%s)' % (self.var.tostr(), self.varlim.tostr(),s,self.expr.tostr())
  75         if self.precedence <= level: r = '(%s)' % (r)
  76         return r
  77
  78 class InfLimit(Basic):
  79     @staticmethod
  80     @cache_it_immutable
  81     def limit_process_symbol():
  82         return Basic.Symbol('xoo', dummy=True, unbounded=True, positive=True)
  83
  84     @cache_it_immutable
  85     def __new__(cls, expr, x):
  86         expr = orig_expr = Basic.sympify(expr)
  87         orig_x = Basic.sympify(x)
  88         assert isinstance(orig_x,Basic.Symbol),`orig_x`
  89
  90         # handle trivial results
  91         if orig_expr==orig_x:
  92             return S.Infinity
  93         elif not orig_expr.has(orig_x):
  94             return orig_expr
  95
  96         x = InfLimit.limit_process_symbol()
  97         if not orig_expr.has(x):
  98             expr = orig_expr.subs(orig_x, x)
  99         elif orig_x==x:
 100             expr = orig_expr
 101         else:
 102             x = Basic.Symbol(orig_x.name + '_oo', dummy=True, unbounded=True, positive=True)
 103             expr = orig_expr.subs(orig_x, x)
 104
 105         result = None
 106         if hasattr(expr,'_eval_inflimit'):
 107             # support for callbacks
 108             result = getattr(expr,'_eval_inflimit')(x)
 109         elif isinstance(expr, Basic.Add):
 110             result, factors = expr.as_coeff_factors(x)
 111             for f in factors:
 112                 result += f.inflimit(x)
 113                 if isinstance(result, Basic.NaN):
 114                     result = None
 115                     break
 116         elif isinstance(expr, Basic.Mul):
 117             result, terms = expr.as_coeff_terms(x)
 118             for t in terms:
 119                 result *= t.inflimit(x)
 120                 if isinstance(result, Basic.NaN):
 121                     result = None
 122                     break
 123         elif isinstance(expr, Basic.Pow):
 124             if not expr.exp.has(x):
 125                 result = expr.base.inflimit(x) ** expr.exp
 126             elif not expr.base.has(x):
 127                 result = expr.base ** expr.exp.inflimit(x)
 128             else:
 129                 result = Basic.exp(expr.exp * Basic.log(expr.base)).inflimit(x)
 130         elif isinstance(expr, Basic.Function):
 131             # warning: assume that
 132             #  lim_x f(g1(x),g2(x),..) = f(lim_x g1(x), lim_x g2(x))
 133             # if this is incorrect, one must define f._eval_inflimit(x) method
 134             result = expr.func(*[a.inflimit(x) for a in expr.args])
 135
 136         if result is None:
 137             result = mrv_inflimit(expr, x)
 138
 139         return result
 140
 141 @cache_it_immutable
 142 def mrv_inflimit(expr, x, _cache = {}):
 143     if _cache.has_key((expr, x)):
 144         raise RuntimeError('Detected recursion while computing mrv_inflimit(%s, %s)' % (expr, x))
 145     _cache[(expr, x)] = 1
 146     expr_map = {}
 147     mrv_map = {}
 148     newexpr = mrv2(expr, x, expr_map, mrv_map)
 149     if mrv_map.has_key(x):
 150         t = Basic.Temporary(unbounded=True, positive=True)
 151         r = mrv_inflimit(expr.subs(Basic.log(x), t).subs(x, Basic.exp(t)).subs(t, x), x)
 152         del _cache[(expr, x)]
 153         return r
 154     w = Basic.Symbol('w_0',dummy=True, positive=True, infinitesimal=True)
 155     germ, new_mrv_map = rewrite_mrv_map(mrv_map, x, w)
 156     new_expr = rewrite_expr(newexpr, germ, new_mrv_map, w)
 157     lt = new_expr.as_leading_term(w)
 158     if germ is not None:
 159         lt = lt.subs(Basic.log(w), -germ.args[0])
 160     c,e = lt.as_coeff_exponent(w)
 161     assert not c.has(w),`c`
 162     if e==0:
 163         r = c.inflimit(x)
 164         del _cache[(expr, x)]
 165         return r
 166     if e.is_positive:
 167         del _cache[(expr, x)]
 168         return S.Zero
 169     if e.is_negative:
 170         del _cache[(expr, x)]
 171         return Basic.sign(c) * S.Infinity
 172     raise RuntimeError('Failed to compute mrv_inflimit(%s, %s), got lt=%s' % (self, x, lt))
 173
 174 @cache_it_immutable
 175 def cmp_ops_count(e1,e2):
 176     return cmp(e1.count_ops(symbolic=False), e2.count_ops(symbolic=False))
 177
 178 @cache_it_immutable
 179 def mrv_compare(f, g, x):
 180     log = Basic.log
 181     if isinstance(f, Basic.exp): f = f.args[0]
 182     else: f = log(f)
 183     if isinstance(g, Basic.exp): g = g.args[0]
 184     else: g = log(g)
 185     c = (f/g).inflimit(x)
 186     if c==0:
 187         return '<'
 188     if isinstance(abs(c), Basic.Infinity):
 189         return '>'
 190     if not c.is_comparable:
 191         raise ValueError("non-comparable result: %s" % (c))
 192     return '='
 193
 194 def mrv2(expr, x, d, md):
 195     """
 196     Compute a set of most rapidly varying subexpressions of expr with respect to x.
 197
 198     d = {}
 199     md = {}
 200     mrv2(x + exp(x),x,d) -> x+se, d={x:x, exp(x):se}, md={exp(x):se}
 201     """
 202     if d.has_key(expr): return d[expr]
 203     if not expr.has(x):
 204         return expr
 205     if expr==x:
 206         if not md: md[x] = x
 207         return x
 208     if isinstance(expr, (Basic.Add, Basic.Mul)):
 209         r = expr.__class__(*[mrv2(t, x, d, md) for t in expr.args])
 210         d[expr] = r
 211         return r
 212     log = Basic.log
 213     exp = Basic.exp
 214     if isinstance(expr, Basic.Pow):
 215         if not expr.exp.has(x):
 216             r = mrv2(expr.base, x, d, md)**expr.exp
 217         else:
 218             r = mrv2(exp(expr.exp * log(expr.base)), x, d, md)
 219         d[expr] = r
 220         return r
 221     if isinstance(expr, Basic.exp):
 222         e = expr.args[0]
 223         l = e.inflimit(x)
 224         r = exp(mrv2(e, x, d, md))
 225         if isinstance(l, Basic.Infinity):
 226             # e -> oo as x->oo
 227             en = e
 228         elif isinstance(l, Basic.NegativeInfinity):
 229             # e -> -oo as x->oo
 230             # rewrite to ensure that exp(e) -> oo
 231             en = -e
 232         else:
 233             # |e| < oo as x->oo
 234             d[expr] = r
 235             return r
 236         # normalize exp(2*e) -> exp(e)
 237         coeff, terms = en.as_coeff_terms()
 238         new_terms = []
 239         for t in terms:
 240             if t.has(x):
 241                 pass
 242             elif t.is_positive:
 243                 continue
 244             elif t.is_negative:
 245                 coeff *= -1
 246                 continue
 247             new_terms.append(t)
 248         terms = new_terms
 249         coeff = Basic.sign(coeff)
 250         if not isinstance(coeff, Basic.One):
 251             terms.insert(0,coeff)
 252         en = Basic.Mul(*terms)
 253         nexpr = exp(en)
 254         #print coeff,terms,nexpr
 255         if md.has_key(x):
 256             c = '>'
 257         else:
 258             lst = md.keys()
 259             lst.sort(cmp_ops_count)
 260             c = mrv_compare(nexpr, lst[0], x)
 261         if c !='<':
 262             if c=='>':
 263                 md.clear()
 264             if md.has_key(nexpr):
 265                 tmp = md[nexpr]
 266             else:
 267                 tmp = Basic.Temporary()
 268                 md[nexpr] = tmp
 269             r = expr.subs(nexpr, tmp)
 270         d[expr] = r
 271         return r
 272     if isinstance(expr, Basic.Function):
 273         r = expr.func(*[mrv2(a, x, d, md) for a in expr.args])
 274         d[expr] = r
 275         return r
 276     raise NotImplementedError("don't know how to find mrv2(%s,%s)" % (expr,x))
 277
 278 def rewrite_mrv_map(mrv_map, x, w):
 279     germs = mrv_map.keys()
 280     germs.sort(cmp_ops_count)
 281     if germs:
 282         g = germs[0]
 283         gname = mrv_map[g]
 284         garg = g.args[0]
 285     else:
 286         g = None
 287     d = {}
 288     for germ in germs:
 289         name = mrv_map[germ]
 290         if name==gname:
 291             d[name] = 1/w
 292             continue
 293         arg = germ.args[0]
 294         c = (arg/garg).inflimit(x)
 295         Aarg = arg-c*garg
 296         Aarg = Aarg.subs(g, 1/w)
 297         A = Basic.exp(Aarg)
 298         new_germ = A * w ** -c
 299         d[name] = new_germ
 300     return g, d
 301
 302 def rewrite_expr(expr, germ, mrv_map, w):
 303     tmps = expr.atoms(Basic.Temporary)
 304     e = expr
 305     for t in tmps:
 306         try:
 307             g = mrv_map[t]
 308         except KeyError:
 309             continue
 310         e = e.subs(t, g)
 311     if germ is not None:
 312         mrvlog = Basic.MrvLog
 313         log = Basic.log
 314         e = e.subs(log, mrvlog).subs(germ.args[0], -log(w)).subs(mrvlog, log)
 315     return e
 316
 317 Basic.singleton['limit'] = lambda : Limit
 318
 319 def limit(e, x, x0):
 320     return e.limit(x, x0)