Coverage for /usr/lib/python3/dist-packages/sympy/series/gruntz.py: 13%
360 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Limits
3======
5Implemented according to the PhD thesis
6https://www.cybertester.com/data/gruntz.pdf, which contains very thorough
7descriptions of the algorithm including many examples. We summarize here
8the gist of it.
10All functions are sorted according to how rapidly varying they are at
11infinity using the following rules. Any two functions f and g can be
12compared using the properties of L:
14L=lim log|f(x)| / log|g(x)| (for x -> oo)
16We define >, < ~ according to::
18 1. f > g .... L=+-oo
20 we say that:
21 - f is greater than any power of g
22 - f is more rapidly varying than g
23 - f goes to infinity/zero faster than g
25 2. f < g .... L=0
27 we say that:
28 - f is lower than any power of g
30 3. f ~ g .... L!=0, +-oo
32 we say that:
33 - both f and g are bounded from above and below by suitable integral
34 powers of the other
36Examples
37========
38::
39 2 < x < exp(x) < exp(x**2) < exp(exp(x))
40 2 ~ 3 ~ -5
41 x ~ x**2 ~ x**3 ~ 1/x ~ x**m ~ -x
42 exp(x) ~ exp(-x) ~ exp(2x) ~ exp(x)**2 ~ exp(x+exp(-x))
43 f ~ 1/f
45So we can divide all the functions into comparability classes (x and x^2
46belong to one class, exp(x) and exp(-x) belong to some other class). In
47principle, we could compare any two functions, but in our algorithm, we
48do not compare anything below the class 2~3~-5 (for example log(x) is
49below this), so we set 2~3~-5 as the lowest comparability class.
51Given the function f, we find the list of most rapidly varying (mrv set)
52subexpressions of it. This list belongs to the same comparability class.
53Let's say it is {exp(x), exp(2x)}. Using the rule f ~ 1/f we find an
54element "w" (either from the list or a new one) from the same
55comparability class which goes to zero at infinity. In our example we
56set w=exp(-x) (but we could also set w=exp(-2x) or w=exp(-3x) ...). We
57rewrite the mrv set using w, in our case {1/w, 1/w^2}, and substitute it
58into f. Then we expand f into a series in w::
60 f = c0*w^e0 + c1*w^e1 + ... + O(w^en), where e0<e1<...<en, c0!=0
62but for x->oo, lim f = lim c0*w^e0, because all the other terms go to zero,
63because w goes to zero faster than the ci and ei. So::
65 for e0>0, lim f = 0
66 for e0<0, lim f = +-oo (the sign depends on the sign of c0)
67 for e0=0, lim f = lim c0
69We need to recursively compute limits at several places of the algorithm, but
70as is shown in the PhD thesis, it always finishes.
72Important functions from the implementation:
74compare(a, b, x) compares "a" and "b" by computing the limit L.
75mrv(e, x) returns list of most rapidly varying (mrv) subexpressions of "e"
76rewrite(e, Omega, x, wsym) rewrites "e" in terms of w
77leadterm(f, x) returns the lowest power term in the series of f
78mrv_leadterm(e, x) returns the lead term (c0, e0) for e
79limitinf(e, x) computes lim e (for x->oo)
80limit(e, z, z0) computes any limit by converting it to the case x->oo
82All the functions are really simple and straightforward except
83rewrite(), which is the most difficult/complex part of the algorithm.
84When the algorithm fails, the bugs are usually in the series expansion
85(i.e. in SymPy) or in rewrite.
87This code is almost exact rewrite of the Maple code inside the Gruntz
88thesis.
90Debugging
91---------
93Because the gruntz algorithm is highly recursive, it's difficult to
94figure out what went wrong inside a debugger. Instead, turn on nice
95debug prints by defining the environment variable SYMPY_DEBUG. For
96example:
98[user@localhost]: SYMPY_DEBUG=True ./bin/isympy
100In [1]: limit(sin(x)/x, x, 0)
101limitinf(_x*sin(1/_x), _x) = 1
102+-mrv_leadterm(_x*sin(1/_x), _x) = (1, 0)
103| +-mrv(_x*sin(1/_x), _x) = set([_x])
104| | +-mrv(_x, _x) = set([_x])
105| | +-mrv(sin(1/_x), _x) = set([_x])
106| | +-mrv(1/_x, _x) = set([_x])
107| | +-mrv(_x, _x) = set([_x])
108| +-mrv_leadterm(exp(_x)*sin(exp(-_x)), _x, set([exp(_x)])) = (1, 0)
109| +-rewrite(exp(_x)*sin(exp(-_x)), set([exp(_x)]), _x, _w) = (1/_w*sin(_w), -_x)
110| +-sign(_x, _x) = 1
111| +-mrv_leadterm(1, _x) = (1, 0)
112+-sign(0, _x) = 0
113+-limitinf(1, _x) = 1
115And check manually which line is wrong. Then go to the source code and
116debug this function to figure out the exact problem.
118"""
119from functools import reduce
121from sympy.core import Basic, S, Mul, PoleError, expand_mul
122from sympy.core.cache import cacheit
123from sympy.core.numbers import ilcm, I, oo
124from sympy.core.symbol import Dummy, Wild
125from sympy.core.traversal import bottom_up
127from sympy.functions import log, exp, sign as _sign
128from sympy.series.order import Order
129from sympy.utilities.exceptions import SymPyDeprecationWarning
130from sympy.utilities.misc import debug_decorator as debug
131from sympy.utilities.timeutils import timethis
133timeit = timethis('gruntz')
136def compare(a, b, x):
137 """Returns "<" if a<b, "=" for a == b, ">" for a>b"""
138 # log(exp(...)) must always be simplified here for termination
139 la, lb = log(a), log(b)
140 if isinstance(a, Basic) and (isinstance(a, exp) or (a.is_Pow and a.base == S.Exp1)):
141 la = a.exp
142 if isinstance(b, Basic) and (isinstance(b, exp) or (b.is_Pow and b.base == S.Exp1)):
143 lb = b.exp
145 c = limitinf(la/lb, x)
146 if c == 0:
147 return "<"
148 elif c.is_infinite:
149 return ">"
150 else:
151 return "="
154class SubsSet(dict):
155 """
156 Stores (expr, dummy) pairs, and how to rewrite expr-s.
158 Explanation
159 ===========
161 The gruntz algorithm needs to rewrite certain expressions in term of a new
162 variable w. We cannot use subs, because it is just too smart for us. For
163 example::
165 > Omega=[exp(exp(_p - exp(-_p))/(1 - 1/_p)), exp(exp(_p))]
166 > O2=[exp(-exp(_p) + exp(-exp(-_p))*exp(_p)/(1 - 1/_p))/_w, 1/_w]
167 > e = exp(exp(_p - exp(-_p))/(1 - 1/_p)) - exp(exp(_p))
168 > e.subs(Omega[0],O2[0]).subs(Omega[1],O2[1])
169 -1/w + exp(exp(p)*exp(-exp(-p))/(1 - 1/p))
171 is really not what we want!
173 So we do it the hard way and keep track of all the things we potentially
174 want to substitute by dummy variables. Consider the expression::
176 exp(x - exp(-x)) + exp(x) + x.
178 The mrv set is {exp(x), exp(-x), exp(x - exp(-x))}.
179 We introduce corresponding dummy variables d1, d2, d3 and rewrite::
181 d3 + d1 + x.
183 This class first of all keeps track of the mapping expr->variable, i.e.
184 will at this stage be a dictionary::
186 {exp(x): d1, exp(-x): d2, exp(x - exp(-x)): d3}.
188 [It turns out to be more convenient this way round.]
189 But sometimes expressions in the mrv set have other expressions from the
190 mrv set as subexpressions, and we need to keep track of that as well. In
191 this case, d3 is really exp(x - d2), so rewrites at this stage is::
193 {d3: exp(x-d2)}.
195 The function rewrite uses all this information to correctly rewrite our
196 expression in terms of w. In this case w can be chosen to be exp(-x),
197 i.e. d2. The correct rewriting then is::
199 exp(-w)/w + 1/w + x.
200 """
201 def __init__(self):
202 self.rewrites = {}
204 def __repr__(self):
205 return super().__repr__() + ', ' + self.rewrites.__repr__()
207 def __getitem__(self, key):
208 if key not in self:
209 self[key] = Dummy()
210 return dict.__getitem__(self, key)
212 def do_subs(self, e):
213 """Substitute the variables with expressions"""
214 for expr, var in self.items():
215 e = e.xreplace({var: expr})
216 return e
218 def meets(self, s2):
219 """Tell whether or not self and s2 have non-empty intersection"""
220 return set(self.keys()).intersection(list(s2.keys())) != set()
222 def union(self, s2, exps=None):
223 """Compute the union of self and s2, adjusting exps"""
224 res = self.copy()
225 tr = {}
226 for expr, var in s2.items():
227 if expr in self:
228 if exps:
229 exps = exps.xreplace({var: res[expr]})
230 tr[var] = res[expr]
231 else:
232 res[expr] = var
233 for var, rewr in s2.rewrites.items():
234 res.rewrites[var] = rewr.xreplace(tr)
235 return res, exps
237 def copy(self):
238 """Create a shallow copy of SubsSet"""
239 r = SubsSet()
240 r.rewrites = self.rewrites.copy()
241 for expr, var in self.items():
242 r[expr] = var
243 return r
246@debug
247def mrv(e, x):
248 """Returns a SubsSet of most rapidly varying (mrv) subexpressions of 'e',
249 and e rewritten in terms of these"""
250 from sympy.simplify.powsimp import powsimp
251 e = powsimp(e, deep=True, combine='exp')
252 if not isinstance(e, Basic):
253 raise TypeError("e should be an instance of Basic")
254 if not e.has(x):
255 return SubsSet(), e
256 elif e == x:
257 s = SubsSet()
258 return s, s[x]
259 elif e.is_Mul or e.is_Add:
260 i, d = e.as_independent(x) # throw away x-independent terms
261 if d.func != e.func:
262 s, expr = mrv(d, x)
263 return s, e.func(i, expr)
264 a, b = d.as_two_terms()
265 s1, e1 = mrv(a, x)
266 s2, e2 = mrv(b, x)
267 return mrv_max1(s1, s2, e.func(i, e1, e2), x)
268 elif e.is_Pow and e.base != S.Exp1:
269 e1 = S.One
270 while e.is_Pow:
271 b1 = e.base
272 e1 *= e.exp
273 e = b1
274 if b1 == 1:
275 return SubsSet(), b1
276 if e1.has(x):
277 base_lim = limitinf(b1, x)
278 if base_lim is S.One:
279 return mrv(exp(e1 * (b1 - 1)), x)
280 return mrv(exp(e1 * log(b1)), x)
281 else:
282 s, expr = mrv(b1, x)
283 return s, expr**e1
284 elif isinstance(e, log):
285 s, expr = mrv(e.args[0], x)
286 return s, log(expr)
287 elif isinstance(e, exp) or (e.is_Pow and e.base == S.Exp1):
288 # We know from the theory of this algorithm that exp(log(...)) may always
289 # be simplified here, and doing so is vital for termination.
290 if isinstance(e.exp, log):
291 return mrv(e.exp.args[0], x)
292 # if a product has an infinite factor the result will be
293 # infinite if there is no zero, otherwise NaN; here, we
294 # consider the result infinite if any factor is infinite
295 li = limitinf(e.exp, x)
296 if any(_.is_infinite for _ in Mul.make_args(li)):
297 s1 = SubsSet()
298 e1 = s1[e]
299 s2, e2 = mrv(e.exp, x)
300 su = s1.union(s2)[0]
301 su.rewrites[e1] = exp(e2)
302 return mrv_max3(s1, e1, s2, exp(e2), su, e1, x)
303 else:
304 s, expr = mrv(e.exp, x)
305 return s, exp(expr)
306 elif e.is_Function:
307 l = [mrv(a, x) for a in e.args]
308 l2 = [s for (s, _) in l if s != SubsSet()]
309 if len(l2) != 1:
310 # e.g. something like BesselJ(x, x)
311 raise NotImplementedError("MRV set computation for functions in"
312 " several variables not implemented.")
313 s, ss = l2[0], SubsSet()
314 args = [ss.do_subs(x[1]) for x in l]
315 return s, e.func(*args)
316 elif e.is_Derivative:
317 raise NotImplementedError("MRV set computation for derivatives"
318 " not implemented yet.")
319 raise NotImplementedError(
320 "Don't know how to calculate the mrv of '%s'" % e)
323def mrv_max3(f, expsf, g, expsg, union, expsboth, x):
324 """
325 Computes the maximum of two sets of expressions f and g, which
326 are in the same comparability class, i.e. max() compares (two elements of)
327 f and g and returns either (f, expsf) [if f is larger], (g, expsg)
328 [if g is larger] or (union, expsboth) [if f, g are of the same class].
329 """
330 if not isinstance(f, SubsSet):
331 raise TypeError("f should be an instance of SubsSet")
332 if not isinstance(g, SubsSet):
333 raise TypeError("g should be an instance of SubsSet")
334 if f == SubsSet():
335 return g, expsg
336 elif g == SubsSet():
337 return f, expsf
338 elif f.meets(g):
339 return union, expsboth
341 c = compare(list(f.keys())[0], list(g.keys())[0], x)
342 if c == ">":
343 return f, expsf
344 elif c == "<":
345 return g, expsg
346 else:
347 if c != "=":
348 raise ValueError("c should be =")
349 return union, expsboth
352def mrv_max1(f, g, exps, x):
353 """Computes the maximum of two sets of expressions f and g, which
354 are in the same comparability class, i.e. mrv_max1() compares (two elements of)
355 f and g and returns the set, which is in the higher comparability class
356 of the union of both, if they have the same order of variation.
357 Also returns exps, with the appropriate substitutions made.
358 """
359 u, b = f.union(g, exps)
360 return mrv_max3(f, g.do_subs(exps), g, f.do_subs(exps),
361 u, b, x)
364@debug
365@cacheit
366@timeit
367def sign(e, x):
368 """
369 Returns a sign of an expression e(x) for x->oo.
371 ::
373 e > 0 for x sufficiently large ... 1
374 e == 0 for x sufficiently large ... 0
375 e < 0 for x sufficiently large ... -1
377 The result of this function is currently undefined if e changes sign
378 arbitrarily often for arbitrarily large x (e.g. sin(x)).
380 Note that this returns zero only if e is *constantly* zero
381 for x sufficiently large. [If e is constant, of course, this is just
382 the same thing as the sign of e.]
383 """
384 if not isinstance(e, Basic):
385 raise TypeError("e should be an instance of Basic")
387 if e.is_positive:
388 return 1
389 elif e.is_negative:
390 return -1
391 elif e.is_zero:
392 return 0
394 elif not e.has(x):
395 from sympy.simplify import logcombine
396 e = logcombine(e)
397 return _sign(e)
398 elif e == x:
399 return 1
400 elif e.is_Mul:
401 a, b = e.as_two_terms()
402 sa = sign(a, x)
403 if not sa:
404 return 0
405 return sa * sign(b, x)
406 elif isinstance(e, exp):
407 return 1
408 elif e.is_Pow:
409 if e.base == S.Exp1:
410 return 1
411 s = sign(e.base, x)
412 if s == 1:
413 return 1
414 if e.exp.is_Integer:
415 return s**e.exp
416 elif isinstance(e, log):
417 return sign(e.args[0] - 1, x)
419 # if all else fails, do it the hard way
420 c0, e0 = mrv_leadterm(e, x)
421 return sign(c0, x)
424@debug
425@timeit
426@cacheit
427def limitinf(e, x):
428 """Limit e(x) for x-> oo."""
429 # rewrite e in terms of tractable functions only
431 old = e
432 if not e.has(x):
433 return e # e is a constant
434 from sympy.simplify.powsimp import powdenest
435 from sympy.calculus.util import AccumBounds
436 if e.has(Order):
437 e = e.expand().removeO()
438 if not x.is_positive or x.is_integer:
439 # We make sure that x.is_positive is True and x.is_integer is None
440 # so we get all the correct mathematical behavior from the expression.
441 # We need a fresh variable.
442 p = Dummy('p', positive=True)
443 e = e.subs(x, p)
444 x = p
445 e = e.rewrite('tractable', deep=True, limitvar=x)
446 e = powdenest(e)
447 if isinstance(e, AccumBounds):
448 if mrv_leadterm(e.min, x) != mrv_leadterm(e.max, x):
449 raise NotImplementedError
450 c0, e0 = mrv_leadterm(e.min, x)
451 else:
452 c0, e0 = mrv_leadterm(e, x)
453 sig = sign(e0, x)
454 if sig == 1:
455 return S.Zero # e0>0: lim f = 0
456 elif sig == -1: # e0<0: lim f = +-oo (the sign depends on the sign of c0)
457 if c0.match(I*Wild("a", exclude=[I])):
458 return c0*oo
459 s = sign(c0, x)
460 # the leading term shouldn't be 0:
461 if s == 0:
462 raise ValueError("Leading term should not be 0")
463 return s*oo
464 elif sig == 0:
465 if c0 == old:
466 c0 = c0.cancel()
467 return limitinf(c0, x) # e0=0: lim f = lim c0
468 else:
469 raise ValueError("{} could not be evaluated".format(sig))
472def moveup2(s, x):
473 r = SubsSet()
474 for expr, var in s.items():
475 r[expr.xreplace({x: exp(x)})] = var
476 for var, expr in s.rewrites.items():
477 r.rewrites[var] = s.rewrites[var].xreplace({x: exp(x)})
478 return r
481def moveup(l, x):
482 return [e.xreplace({x: exp(x)}) for e in l]
485@debug
486@timeit
487def calculate_series(e, x, logx=None):
488 """ Calculates at least one term of the series of ``e`` in ``x``.
490 This is a place that fails most often, so it is in its own function.
491 """
493 SymPyDeprecationWarning(
494 feature="calculate_series",
495 useinstead="series() with suitable n, or as_leading_term",
496 issue=21838,
497 deprecated_since_version="1.12"
498 ).warn()
500 from sympy.simplify.powsimp import powdenest
502 for t in e.lseries(x, logx=logx):
503 # bottom_up function is required for a specific case - when e is
504 # -exp(p/(p + 1)) + exp(-p**2/(p + 1) + p)
505 t = bottom_up(t, lambda w:
506 getattr(w, 'normal', lambda: w)())
507 # And the expression
508 # `(-sin(1/x) + sin((x + exp(x))*exp(-x)/x))*exp(x)`
509 # from the first test of test_gruntz_eval_special needs to
510 # be expanded. But other forms need to be have at least
511 # factor_terms applied. `factor` accomplishes both and is
512 # faster than using `factor_terms` for the gruntz suite. It
513 # does not appear that use of `cancel` is necessary.
514 # t = cancel(t, expand=False)
515 t = t.factor()
517 if t.has(exp) and t.has(log):
518 t = powdenest(t)
520 if not t.is_zero:
521 break
523 return t
526@debug
527@timeit
528@cacheit
529def mrv_leadterm(e, x):
530 """Returns (c0, e0) for e."""
531 Omega = SubsSet()
532 if not e.has(x):
533 return (e, S.Zero)
534 if Omega == SubsSet():
535 Omega, exps = mrv(e, x)
536 if not Omega:
537 # e really does not depend on x after simplification
538 return exps, S.Zero
539 if x in Omega:
540 # move the whole omega up (exponentiate each term):
541 Omega_up = moveup2(Omega, x)
542 exps_up = moveup([exps], x)[0]
543 # NOTE: there is no need to move this down!
544 Omega = Omega_up
545 exps = exps_up
546 #
547 # The positive dummy, w, is used here so log(w*2) etc. will expand;
548 # a unique dummy is needed in this algorithm
549 #
550 # For limits of complex functions, the algorithm would have to be
551 # improved, or just find limits of Re and Im components separately.
552 #
553 w = Dummy("w", positive=True)
554 f, logw = rewrite(exps, Omega, x, w)
555 try:
556 lt = f.leadterm(w, logx=logw)
557 except (NotImplementedError, PoleError, ValueError):
558 n0 = 1
559 _series = Order(1)
560 incr = S.One
561 while _series.is_Order:
562 _series = f._eval_nseries(w, n=n0+incr, logx=logw)
563 incr *= 2
564 series = _series.expand().removeO()
565 try:
566 lt = series.leadterm(w, logx=logw)
567 except (NotImplementedError, PoleError, ValueError):
568 lt = f.as_coeff_exponent(w)
569 if lt[0].has(w):
570 base = f.as_base_exp()[0].as_coeff_exponent(w)
571 ex = f.as_base_exp()[1]
572 lt = (base[0]**ex, base[1]*ex)
573 return (lt[0].subs(log(w), logw), lt[1])
576def build_expression_tree(Omega, rewrites):
577 r""" Helper function for rewrite.
579 We need to sort Omega (mrv set) so that we replace an expression before
580 we replace any expression in terms of which it has to be rewritten::
582 e1 ---> e2 ---> e3
583 \
584 -> e4
586 Here we can do e1, e2, e3, e4 or e1, e2, e4, e3.
587 To do this we assemble the nodes into a tree, and sort them by height.
589 This function builds the tree, rewrites then sorts the nodes.
590 """
591 class Node:
592 def __init__(self):
593 self.before = []
594 self.expr = None
595 self.var = None
596 def ht(self):
597 return reduce(lambda x, y: x + y,
598 [x.ht() for x in self.before], 1)
599 nodes = {}
600 for expr, v in Omega:
601 n = Node()
602 n.var = v
603 n.expr = expr
604 nodes[v] = n
605 for _, v in Omega:
606 if v in rewrites:
607 n = nodes[v]
608 r = rewrites[v]
609 for _, v2 in Omega:
610 if r.has(v2):
611 n.before.append(nodes[v2])
613 return nodes
616@debug
617@timeit
618def rewrite(e, Omega, x, wsym):
619 """e(x) ... the function
620 Omega ... the mrv set
621 wsym ... the symbol which is going to be used for w
623 Returns the rewritten e in terms of w and log(w). See test_rewrite1()
624 for examples and correct results.
625 """
627 from sympy import AccumBounds
628 if not isinstance(Omega, SubsSet):
629 raise TypeError("Omega should be an instance of SubsSet")
630 if len(Omega) == 0:
631 raise ValueError("Length cannot be 0")
632 # all items in Omega must be exponentials
633 for t in Omega.keys():
634 if not isinstance(t, exp):
635 raise ValueError("Value should be exp")
636 rewrites = Omega.rewrites
637 Omega = list(Omega.items())
639 nodes = build_expression_tree(Omega, rewrites)
640 Omega.sort(key=lambda x: nodes[x[1]].ht(), reverse=True)
642 # make sure we know the sign of each exp() term; after the loop,
643 # g is going to be the "w" - the simplest one in the mrv set
644 for g, _ in Omega:
645 sig = sign(g.exp, x)
646 if sig != 1 and sig != -1 and not sig.has(AccumBounds):
647 raise NotImplementedError('Result depends on the sign of %s' % sig)
648 if sig == 1:
649 wsym = 1/wsym # if g goes to oo, substitute 1/w
650 # O2 is a list, which results by rewriting each item in Omega using "w"
651 O2 = []
652 denominators = []
653 for f, var in Omega:
654 c = limitinf(f.exp/g.exp, x)
655 if c.is_Rational:
656 denominators.append(c.q)
657 arg = f.exp
658 if var in rewrites:
659 if not isinstance(rewrites[var], exp):
660 raise ValueError("Value should be exp")
661 arg = rewrites[var].args[0]
662 O2.append((var, exp((arg - c*g.exp).expand())*wsym**c))
664 # Remember that Omega contains subexpressions of "e". So now we find
665 # them in "e" and substitute them for our rewriting, stored in O2
667 # the following powsimp is necessary to automatically combine exponentials,
668 # so that the .xreplace() below succeeds:
669 # TODO this should not be necessary
670 from sympy.simplify.powsimp import powsimp
671 f = powsimp(e, deep=True, combine='exp')
672 for a, b in O2:
673 f = f.xreplace({a: b})
675 for _, var in Omega:
676 assert not f.has(var)
678 # finally compute the logarithm of w (logw).
679 logw = g.exp
680 if sig == 1:
681 logw = -logw # log(w)->log(1/w)=-log(w)
683 # Some parts of SymPy have difficulty computing series expansions with
684 # non-integral exponents. The following heuristic improves the situation:
685 exponent = reduce(ilcm, denominators, 1)
686 f = f.subs({wsym: wsym**exponent})
687 logw /= exponent
689 # bottom_up function is required for a specific case - when f is
690 # -exp(p/(p + 1)) + exp(-p**2/(p + 1) + p). No current simplification
691 # methods reduce this to 0 while not expanding polynomials.
692 f = bottom_up(f, lambda w: getattr(w, 'normal', lambda: w)())
693 f = expand_mul(f)
695 return f, logw
698def gruntz(e, z, z0, dir="+"):
699 """
700 Compute the limit of e(z) at the point z0 using the Gruntz algorithm.
702 Explanation
703 ===========
705 ``z0`` can be any expression, including oo and -oo.
707 For ``dir="+"`` (default) it calculates the limit from the right
708 (z->z0+) and for ``dir="-"`` the limit from the left (z->z0-). For infinite z0
709 (oo or -oo), the dir argument does not matter.
711 This algorithm is fully described in the module docstring in the gruntz.py
712 file. It relies heavily on the series expansion. Most frequently, gruntz()
713 is only used if the faster limit() function (which uses heuristics) fails.
714 """
715 if not z.is_symbol:
716 raise NotImplementedError("Second argument must be a Symbol")
718 # convert all limits to the limit z->oo; sign of z is handled in limitinf
719 r = None
720 if z0 in (oo, I*oo):
721 e0 = e
722 elif z0 in (-oo, -I*oo):
723 e0 = e.subs(z, -z)
724 else:
725 if str(dir) == "-":
726 e0 = e.subs(z, z0 - 1/z)
727 elif str(dir) == "+":
728 e0 = e.subs(z, z0 + 1/z)
729 else:
730 raise NotImplementedError("dir must be '+' or '-'")
732 r = limitinf(e0, z)
734 # This is a bit of a heuristic for nice results... we always rewrite
735 # tractable functions in terms of familiar intractable ones.
736 # It might be nicer to rewrite the exactly to what they were initially,
737 # but that would take some work to implement.
738 return r.rewrite('intractable', deep=True)