sympy/series/limits2.py

   1 # This file is the original (unmodified) limits.py from the oldcore
   2 # please don't modify it, it's here for the reference
   3 # see limits.py for more info
   4
   5 """
   6 Limits
   7 ======
   8
   9 Implemented according to the PhD thesis
  10 http://www.cybertester.com/data/gruntz.pdf, which contains very thorough
  11 descriptions of the algorithm including many examples.  We summarize here the
  12 gist of it.
  13
  14
  15 All functions are sorted according to how rapidly varying they are at infinity
  16 using the following rules. Any two functions f and g can be compared using the
  17 properties of L:
  18
  19 L=lim  log|f(x)| / log|g(x)|           (for x -> oo)
  20
  21 We define >, < ~ according to::
  22
  23     1. f > g .... L=+-oo
  24
  25         we say that:
  26         - f is greater than any power of g
  27         - f is more rapidly varying than g
  28         - f goes to infinity/zero faster than g
  29
  30
  31     2. f < g .... L=0
  32
  33         we say that:
  34         - f is lower than any power of g
  35
  36     3. f ~ g .... L!=0,+-oo
  37
  38         we say that:
  39         - both f and g are bounded from above and below by suitable integral
  40           powers of the other
  41
  42
  43 Examples
  44 ========
  45 ::
  46     2 < x < exp(x) < exp(x**2) < exp(exp(x))
  47     2 ~ 3 ~ -5
  48     x ~ x**2 ~ x**3 ~ 1/x ~ x**m ~ -x
  49     exp(x) ~ exp(-x) ~ exp(2x) ~ exp(x)**2 ~ exp(x+exp(-x))
  50     f ~ 1/f
  51
  52 So we can divide all the functions into comparability classes (x and x^2 belong
  53 to one class, exp(x) and exp(-x) belong to some other class). In principle, we
  54 could compare any two functions, but in our algorithm, we don't compare
  55 anything below the class 2~3~-5 (for example log(x) is below this), so we set
  56 2~3~-5 as the lowest comparability class.
  57
  58 Given the function f, we find the list of most rapidly varying (mrv set)
  59 subexpressions of it. This list belongs to the same comparability class. Let's
  60 say it is {exp(x), exp(2x)}. Using the rule f ~ 1/f we find an element "w"
  61 (either from the list or a new one) from the same comparability class which
  62 goes to zero at infinity. In our example we set w=exp(-x) (but we could also
  63 set w=exp(-2x) or w=exp(-3x) ...). We rewrite the mrv set using w, in our case
  64 {1/w,1/w^2}, and substitute it into f. Then we expand f into a series in w::
  65
  66     f = c0*w^e0 + c1*w^e1 + ... + O(w^en),        where e0<e1<...<en, c0!=0
  67
  68 but for x->oo, lim f = lim c0*w^e0, because all the other terms go to zero,
  69 because w goes to zero faster than the ci and ei. So::
  70
  71     for e0>0, lim f = 0
  72     for e0<0, lim f = +-oo   (the sign depends on the sign of c0)
  73     for e0=0, lim f = lim c0
  74
  75 We need to recursively compute limits at several places of the algorithm, but
  76 as is shown in the PhD thesis, it always finishes.
  77
  78 Important functions from the implementation:
  79
  80 compare(a,b,x) compares "a" and "b" by computing the limit L.
  81 mrv(e,x) returns the list of most rapidly varying (mrv) subexpressions of "e"
  82 rewrite(e,Omega,x,wsym) rewrites "e" in terms of w
  83 leadterm(f,x) returns the lowest power term in the series of f
  84 mrv_leadterm(e,x) returns the lead term (c0,e0) for e
  85 limitinf(e,x) computes lim e  (for x->oo)
  86 limit(e,z,z0) computes any limit by converting it to the case x->oo
  87
  88 All the functions are really simple and straightforward except rewrite(), which
  89 is the most difficult/complex part of the algorithm. When the algorithm fails,
  90 the bugs are usually in the series expansion (i.e. in SymPy) or in rewrite.
  91
  92 This code is almost exact rewrite of the Maple code inside the Gruntz thesis.
  93
  94 """
  95
  96 import sympy
  97 from sympy import Basic, Add, Mul, Pow, Function, log, oo, Rational, exp, \
  98         Real, Order, Symbol
  99 O = Order
 100
 101 def debug(func):
 102     def decorated(*args, **kwargs):
 103         #r = func(*args, **kwargs)
 104         r = maketree(func, *args, **kwargs)
 105         #print "%s = %s(%s, %s)" % (r, func.__name__, args, kwargs)
 106         return r
 107     if 1:
 108         #normal mode
 109         return lambda *args, **kwargs: func(*args, **kwargs)
 110     else:
 111         #debug mode
 112         return decorated
 113
 114 def tree(subtrees):
 115     def indent(s,type=1):
 116         x = s.split("\n")
 117         r = "+-%s\n"%x[0]
 118         for a in x[1:]:
 119             if a=="": continue
 120             if type==1:
 121                 r += "| %s\n"%a
 122             else:
 123                 r += "  %s\n"%a
 124         return r
 125     if len(subtrees)==0: return ""
 126     f="";
 127     for a in subtrees[:-1]:
 128         f += indent(a)
 129     f += indent(subtrees[-1],2)
 130     return f
 131
 132 tmp=[]
 133 iter=0
 134 def maketree(f,*args,**kw):
 135     global tmp
 136     global iter
 137     oldtmp=tmp
 138     tmp=[]
 139     iter+=1
 140
 141     r = f(*args,**kw)
 142
 143     iter-=1
 144     s = "%s%s = %s\n" % (f.func_name,args,r)
 145     if tmp!=[]: s += tree(tmp)
 146     tmp=oldtmp
 147     tmp.append(s)
 148     if iter == 0:
 149         print tmp[0]
 150         tmp=[]
 151     return r
 152
 153 def compare(a,b,x):
 154     """Returns "<" if a<b, "=" for a==b, ">" for a>b"""
 155     #use the sympy's broken limit as the starting point (bootstrapping) :)
 156     #c = limitinf(log(a)/log(b), x)
 157     c = (log(a)/log(b)).inflimit(x)
 158     if c == 0:
 159         return "<"
 160     elif c in [oo,-oo]:
 161         return ">"
 162     else:
 163         return "="
 164
 165 @debug
 166 def mrv(e, x):
 167     "Returns a python set of  most rapidly varying (mrv) subexpressions of 'e'"
 168     assert isinstance(e, Basic)
 169     if not e.has(x):
 170         return set([])
 171     elif e == x:
 172         return set([x])
 173     elif isinstance(e, Mul):
 174         a, b = e.as_two_terms()
 175         return mrv_max(mrv(a,x), mrv(b,x), x)
 176     elif isinstance(e, Add):
 177         a, b = e.as_two_terms()
 178         return mrv_max(mrv(a,x), mrv(b,x), x)
 179     elif isinstance(e, Pow):
 180         if e.exp.has(x):
 181             return mrv(exp(e.exp * log(e.base)), x)
 182         else:
 183             return mrv(e.base, x)
 184     elif isinstance(e, log):
 185         return mrv(e[0], x)
 186     elif isinstance(e, exp):
 187         if e[0].inflimit(x) in [oo,-oo]:
 188             return mrv_max(set([e]), mrv(e[0], x), x)
 189         else:
 190             return mrv(e[0], x)
 191     elif isinstance(e, Function):
 192         if len(e) == 1:
 193             return mrv(e[0], x)
 194         #only functions of 1 argument currently implemented
 195         raise NotImplementedError("Functions with more arguments: '%s'" % e)
 196     raise NotImplementedError("Don't know how to calculate the mrv of '%s'" % e)
 197
 198 def mrv_max(f, g, x):
 199     """Computes the maximum of two sets of expressions f and g, which
 200     are in the same comparability class, i.e. max() compares (two elements of)
 201     f and g and returns the set, which is in the higher comparability class
 202     of the union of both, if they have the same order of variation.
 203     """
 204     assert isinstance(f, set)
 205     assert isinstance(g, set)
 206     if f==set([]): return g
 207     elif g==set([]): return f
 208     elif f.intersection(g) != set([]): return f.union(g)
 209     elif x in f: return g
 210     elif x in g: return f
 211
 212     c=compare(list(f)[0], list(g)[0], x)
 213     if c == ">": return f
 214     elif c == "<": return g
 215     else:
 216         assert c == "="
 217         return f.union(g)
 218
 219 @debug
 220 def rewrite(e,Omega,x,wsym):
 221     """e(x) ... the function
 222     Omega ... the mrv set
 223     wsym ... the symbol which is going to be used for w
 224
 225     Returns the rewritten e in terms of w and log(w). See test_rewrite1()
 226     for examples and correct results.
 227     """
 228     assert isinstance(Omega, set)
 229     assert len(Omega)!=0
 230     #all items in Omega must be exponentials
 231     for t in Omega: assert isinstance(t, exp)
 232     def cmpfunc(a,b):
 233         #FIXME: this is really, really slow...
 234         return -cmp(len(mrv(a,x)), len(mrv(b,x)))
 235     #sort Omega (mrv set) from the most complicated to the simplest ones
 236     #the complexity of "a" from Omega: the length of the mrv set of "a"
 237     Omega = list(Omega)
 238     Omega.sort(cmp=cmpfunc)
 239     g=Omega[-1] #g is going to be the "w" - the simplest one in the mrv set
 240     sig = (sign(g[0], x) == 1)
 241     if sig: wsym=1/wsym #if g goes to oo, substitute 1/w
 242     #O2 is a list, which results by rewriting each item in Omega using "w"
 243     O2=[]
 244     for f in Omega:
 245         c=mrv_leadterm(f[0]/g[0], x)
 246         #the c is a constant, because both f and g are from Omega:
 247         assert c[1] == 0
 248         O2.append(exp((f[0]-c[0]*g[0]).expand())*wsym**c[0])
 249     #Remember that Omega contains subexpressions of "e". So now we find
 250     #them in "e" and substitute them for our rewriting, stored in O2
 251     f=e
 252     for a,b in zip(Omega,O2):
 253         f=f.subs(a,b)
 254
 255     #tmp.append("Omega=%s; O2=%s; w=%s; wsym=%s\n"%(Omega,O2,g,wsym))
 256
 257     #finally compute the logarithm of w (logw).
 258     logw=g[0]
 259     if sig: logw=-logw     #log(w)->log(1/w)=-log(w)
 260     return f, logw
 261
 262 @debug
 263 def sign(e, x):
 264     """Returns a sign of an expression e(x) for x->oo.
 265
 266         e >  0 ...  1
 267         e == 0 ...  0
 268         e <  0 ... -1
 269     """
 270     assert isinstance(e, Basic)
 271     if isinstance(e, (Rational, Real)):
 272         if e == 0:
 273             return 0
 274         elif e.evalf() > 0:
 275             return 1
 276         else:
 277             return -1
 278     elif not e.has(x):
 279         f= e.evalf()
 280         if f > 0:
 281             return 1
 282         else:
 283             return -1
 284     elif e == x:
 285         return 1
 286     elif isinstance(e, Mul):
 287         a,b = e.as_two_terms()
 288         return sign(a, x) * sign(b, x)
 289     elif isinstance(e, exp):
 290         return 1
 291     elif isinstance(e, Pow):
 292         if sign(e.base, x) == 1:
 293             return 1
 294     elif isinstance(e, log):
 295         return sign(e[0] -1, x)
 296     elif isinstance(e, Add):
 297         return sign(e.inflimit(x), x)
 298     raise "cannot determine the sign of %s"%e
 299
 300 @debug
 301 def limitinf(e, x):
 302     """Limit e(x) for x-> oo"""
 303     if not e.has(x): return e #e is a constant
 304     c0, e0 = mrv_leadterm(e,x)
 305     sig=sign(e0,x)
 306     if sig==1: return Rational(0) # e0>0: lim f = 0
 307     elif sig==-1: #e0<0: lim f = +-oo   (the sign depends on the sign of c0)
 308         s = sign(c0, x)
 309         #the leading term shouldn't be 0:
 310         assert s != 0
 311         return s * oo
 312     elif sig==0: return limitinf(c0, x) #e0=0: lim f = lim c0
 313
 314 def moveup(l, x):
 315     return [e.subs(x,exp(x)) for e in l]
 316
 317 def movedown(l, x):
 318     return [e.subs(x,log(x)) for e in l]
 319
 320 def subexp(e,sub):
 321     """Is "sub" a subexpression of "e"? """
 322     #we substitute some symbol for the "sub", and if the
 323     #expression changes, the substitution was successful, thus the answer
 324     #is yes.
 325     return e.subs(sub, Symbol("x", dummy=True)) != e
 326
 327 @debug
 328 def mrv_leadterm(e, x, Omega=[]):
 329     """Returns (c0, e0) for e."""
 330     if not e.has(x): return (e, 0)
 331     Omega = [t for t in Omega if subexp(e,t)]
 332     if Omega == []:
 333         Omega = mrv(e,x)
 334     if x in set(Omega):
 335         #move the whole omega up (exponentiate each term):
 336         Omega_up = set(moveup(Omega,x))
 337         e_up = moveup([e],x)[0]
 338         #calculate the lead term
 339         mrv_leadterm_up = mrv_leadterm(e_up, x, Omega_up)
 340         #move the result (c0, e0) down
 341         return tuple(movedown(mrv_leadterm_up, x))
 342     wsym = Symbol("w", dummy=True)
 343     f, logw=rewrite(e, set(Omega), x, wsym)
 344     series=f.expand().oseries(O(wsym**2, wsym))
 345     assert series!=0
 346     assert not isinstance(series,O)
 347     #print "sss1",series,type(series),f,n
 348     #series = series.removeO()
 349     #print "sss2",series,type(series)
 350     series=series.subs(log(wsym), logw)
 351     #print "sss3",series,type(series)
 352     return series.leadterm(wsym)
 353
 354
 355 class Limit2(Basic):
 356
 357     mathml_tag = 'limit'
 358
 359     def __init__(self,e,x,x0):
 360         Basic.__init__(self)
 361         self._args = list()
 362         self._args.append(self.sympify(e))
 363         self._args.append(self.sympify(x))
 364         self._args.append(self.sympify(x0))
 365
 366
 367     def __pretty__(self):
 368          e, x, t = [a.__pretty__() for a in (self.e,self.x,self.x0)]
 369          a = prettyForm('lim')
 370          a = prettyForm(*a.below('%s->%s' % (x, t)))
 371          a = prettyForm(*stringPict.next(a, e))
 372          return a
 373
 374     def __latex__(self):
 375          return r"\lim_{%s \to %s}%s" % (self.x.__latex__(), \
 376                                                  self.x0.__latex__(),
 377                                                  self.e.__latex__() )
 378
 379     @property
 380     def e(self):
 381         return self._args[0]
 382
 383     @property
 384     def x(self):
 385         return self._args[1]
 386
 387     @property
 388     def x0(self):
 389         return self._args[2]
 390
 391     def doit(self):
 392         return limit(self.e,self.x,self.x0)
 393
 394     def __mathml__(self):
 395         if self._mathml:
 396             return self._mathml
 397         import xml.dom.minidom
 398         dom = xml.dom.minidom.Document()
 399         x = dom.createElement("apply")
 400         x.appendChild(dom.createElement(self.mathml_tag))
 401
 402         x_1 = dom.createElement('bvar')
 403
 404         x_2 = dom.createElement('lowlimit')
 405
 406         x.appendChild(x_1)
 407         x.appendChild(x_2)
 408         x.appendChild( self.e.__mathml__() )
 409         x.childNodes[1].appendChild( self.x.__mathml__() )
 410         x.childNodes[2].appendChild( self.x0.__mathml__() )
 411         self._mathml = x
 412
 413         return self._mathml
 414
 415 def limit(e, z, z0, dir="+"):
 416     """
 417     Compute the limit of e(z) at the point z0.
 418
 419     z0 can be any expression, including oo and -oo.
 420
 421     For dir="+" (default) it calculates the limit from the right
 422     (z->z0+) and for dir="-" the limit from the left (z->z0-). For infinite z0
 423     (oo or -oo), the dir argument doesn't matter.
 424     """
 425     if not isinstance(z, Symbol):
 426         raise NotImplementedError("Second argument must be a Symbol")
 427
 428     #convert all limits to the limit z->oo
 429     elif z0 == oo:
 430         return limitinf(e, z)
 431     elif z0 == -oo:
 432         return limitinf(e.subs(z,-z), z)
 433     else:
 434         x = Symbol("x", dummy=True)
 435         if dir == "-":
 436             e0 = e.subs(z,z0-1/x)
 437         elif dir == "+":
 438             e0 = e.subs(z,z0+1/x)
 439         else:
 440             raise NotImplementedError("dir must be '+' or '-'")
 441         return limitinf(e0, x)