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
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
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
19 L=lim log|f(x)| / log|g(x)| (for x -> oo)
21 We define >, < ~ according to::
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
34 - f is lower than any power of g
36 3. f ~ g .... L!=0,+-oo
39 - both f and g are bounded from above and below by suitable integral
46 2 < x < exp(x) < exp(x**2) < exp(exp(x))
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))
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.
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::
66 f = c0*w^e0 + c1*w^e1 + ... + O(w^en), where e0<e1<...<en, c0!=0
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::
72 for e0<0, lim f = +-oo (the sign depends on the sign of c0)
73 for e0=0, lim f = lim c0
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.
78 Important functions from the implementation:
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
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.
92 This code is almost exact rewrite of the Maple code inside the Gruntz thesis.
97 from sympy
import Basic
, Add
, Mul
, Pow
, Function
, log
, oo
, Rational
, exp
, \
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)
109 return lambda *args
, **kwargs
: func(*args
, **kwargs
)
115 def indent(s
,type=1):
125 if len(subtrees
)==0: return ""
127 for a
in subtrees
[:-1]:
129 f
+= indent(subtrees
[-1],2)
134 def maketree(f
,*args
,**kw
):
144 s
= "%s%s = %s\n" % (f
.func_name
,args
,r
)
145 if tmp
!=[]: s
+= tree(tmp
)
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
)
167 "Returns a python set of most rapidly varying (mrv) subexpressions of 'e'"
168 assert isinstance(e
, Basic
)
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
):
181 return mrv(exp(e
.exp
* log(e
.base
)), x
)
183 return mrv(e
.base
, x
)
184 elif isinstance(e
, log
):
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
)
191 elif isinstance(e
, Function
):
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
)
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.
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
212 c
=compare(list(f
)[0], list(g
)[0], x
)
213 if c
== ">": return f
214 elif c
== "<": return g
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
225 Returns the rewritten e in terms of w and log(w). See test_rewrite1()
226 for examples and correct results.
228 assert isinstance(Omega
, set)
230 #all items in Omega must be exponentials
231 for t
in Omega
: assert isinstance(t
, exp
)
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"
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"
245 c
=mrv_leadterm(f
[0]/g
[0], x
)
246 #the c is a constant, because both f and g are from Omega:
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
252 for a
,b
in zip(Omega
,O2
):
255 #tmp.append("Omega=%s; O2=%s; w=%s; wsym=%s\n"%(Omega,O2,g,wsym))
257 #finally compute the logarithm of w (logw).
259 if sig
: logw
=-logw
#log(w)->log(1/w)=-log(w)
264 """Returns a sign of an expression e(x) for x->oo.
270 assert isinstance(e
, Basic
)
271 if isinstance(e
, (Rational
, Real
)):
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
):
291 elif isinstance(e
, Pow
):
292 if sign(e
.base
, x
) == 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
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
)
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)
309 #the leading term shouldn't be 0:
312 elif sig
==0: return limitinf(c0
, x
) #e0=0: lim f = lim c0
315 return [e
.subs(x
,exp(x
)) for e
in l
]
318 return [e
.subs(x
,log(x
)) for e
in l
]
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
325 return e
.subs(sub
, Symbol("x", dummy
=True)) != e
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
)]
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
))
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
)
359 def __init__(self
,e
,x
,x0
):
362 self
._args
.append(self
.sympify(e
))
363 self
._args
.append(self
.sympify(x
))
364 self
._args
.append(self
.sympify(x0
))
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
))
375 return r
"\lim_{%s \to %s}%s" % (self
.x
.__latex
__(), \
392 return limit(self
.e
,self
.x
,self
.x0
)
394 def __mathml__(self
):
397 import xml
.dom
.minidom
398 dom
= xml
.dom
.minidom
.Document()
399 x
= dom
.createElement("apply")
400 x
.appendChild(dom
.createElement(self
.mathml_tag
))
402 x_1
= dom
.createElement('bvar')
404 x_2
= dom
.createElement('lowlimit')
408 x
.appendChild( self
.e
.__mathml
__() )
409 x
.childNodes
[1].appendChild( self
.x
.__mathml
__() )
410 x
.childNodes
[2].appendChild( self
.x0
.__mathml__() )
415 def limit(e
, z
, z0
, dir="+"):
417 Compute the limit of e(z) at the point z0.
419 z0 can be any expression, including oo and -oo.
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.
425 if not isinstance(z
, Symbol
):
426 raise NotImplementedError("Second argument must be a Symbol")
428 #convert all limits to the limit z->oo
430 return limitinf(e
, z
)
432 return limitinf(e
.subs(z
,-z
), z
)
434 x
= Symbol("x", dummy
=True)
436 e0
= e
.subs(z
,z0
-1/x
)
438 e0
= e
.subs(z
,z0
+1/x
)
440 raise NotImplementedError("dir must be '+' or '-'")
441 return limitinf(e0
, x
)