1""" 2Limits 3====== 4 5Implemented according to the PhD thesis 6http://www.cybertester.com/data/gruntz.pdf, which contains very thorough 7descriptions of the algorithm including many examples. We summarize here 8the gist of it. 9 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: 13 14L=lim log|f(x)| / log|g(x)| (for x -> oo) 15 16We define >, < ~ according to:: 17 18 1. f > g .... L=+-oo 19 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 24 25 2. f < g .... L=0 26 27 we say that: 28 - f is lower than any power of g 29 30 3. f ~ g .... L!=0, +-oo 31 32 we say that: 33 - both f and g are bounded from above and below by suitable integral 34 powers of the other 35 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 44 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 48don't 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. 50 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:: 59 60 f = c0*w^e0 + c1*w^e1 + ... + O(w^en), where e0<e1<...<en, c0!=0 61 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:: 64 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 68 69We need to recursively compute limits at several places of the algorithm, but 70as is shown in the PhD thesis, it always finishes. 71 72Important functions from the implementation: 73 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 81 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. 86 87This code is almost exact rewrite of the Maple code inside the Gruntz 88thesis. 89 90Debugging 91--------- 92 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: 97 98[user@localhost]: SYMPY_DEBUG=True ./bin/isympy 99 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 114 115And check manually which line is wrong. Then go to the source code and 116debug this function to figure out the exact problem. 117 118""" 119from functools import reduce 120 121from sympy import cacheit 122from sympy.core import Basic, S, oo, I, Dummy, Wild, Mul, PoleError 123 124from sympy.functions import log, exp 125from sympy.series.order import Order 126from sympy.simplify import logcombine 127from sympy.simplify.powsimp import powsimp, powdenest 128 129from sympy.utilities.misc import debug_decorator as debug 130from sympy.utilities.timeutils import timethis 131timeit = timethis('gruntz') 132 133 134 135def compare(a, b, x): 136 """Returns "<" if a<b, "=" for a == b, ">" for a>b""" 137 # log(exp(...)) must always be simplified here for termination 138 la, lb = log(a), log(b) 139 if isinstance(a, Basic) and (isinstance(a, exp) or (a.is_Pow and a.base == S.Exp1)): 140 la = a.exp 141 if isinstance(b, Basic) and (isinstance(b, exp) or (b.is_Pow and b.base == S.Exp1)): 142 lb = b.exp 143 144 c = limitinf(la/lb, x) 145 if c == 0: 146 return "<" 147 elif c.is_infinite: 148 return ">" 149 else: 150 return "=" 151 152 153class SubsSet(dict): 154 """ 155 Stores (expr, dummy) pairs, and how to rewrite expr-s. 156 157 Explanation 158 =========== 159 160 The gruntz algorithm needs to rewrite certain expressions in term of a new 161 variable w. We cannot use subs, because it is just too smart for us. For 162 example:: 163 164 > Omega=[exp(exp(_p - exp(-_p))/(1 - 1/_p)), exp(exp(_p))] 165 > O2=[exp(-exp(_p) + exp(-exp(-_p))*exp(_p)/(1 - 1/_p))/_w, 1/_w] 166 > e = exp(exp(_p - exp(-_p))/(1 - 1/_p)) - exp(exp(_p)) 167 > e.subs(Omega[0],O2[0]).subs(Omega[1],O2[1]) 168 -1/w + exp(exp(p)*exp(-exp(-p))/(1 - 1/p)) 169 170 is really not what we want! 171 172 So we do it the hard way and keep track of all the things we potentially 173 want to substitute by dummy variables. Consider the expression:: 174 175 exp(x - exp(-x)) + exp(x) + x. 176 177 The mrv set is {exp(x), exp(-x), exp(x - exp(-x))}. 178 We introduce corresponding dummy variables d1, d2, d3 and rewrite:: 179 180 d3 + d1 + x. 181 182 This class first of all keeps track of the mapping expr->variable, i.e. 183 will at this stage be a dictionary:: 184 185 {exp(x): d1, exp(-x): d2, exp(x - exp(-x)): d3}. 186 187 [It turns out to be more convenient this way round.] 188 But sometimes expressions in the mrv set have other expressions from the 189 mrv set as subexpressions, and we need to keep track of that as well. In 190 this case, d3 is really exp(x - d2), so rewrites at this stage is:: 191 192 {d3: exp(x-d2)}. 193 194 The function rewrite uses all this information to correctly rewrite our 195 expression in terms of w. In this case w can be chosen to be exp(-x), 196 i.e. d2. The correct rewriting then is:: 197 198 exp(-w)/w + 1/w + x. 199 """ 200 def __init__(self): 201 self.rewrites = {} 202 203 def __repr__(self): 204 return super().__repr__() + ', ' + self.rewrites.__repr__() 205 206 def __getitem__(self, key): 207 if not key in self: 208 self[key] = Dummy() 209 return dict.__getitem__(self, key) 210 211 def do_subs(self, e): 212 """Substitute the variables with expressions""" 213 for expr, var in self.items(): 214 e = e.xreplace({var: expr}) 215 return e 216 217 def meets(self, s2): 218 """Tell whether or not self and s2 have non-empty intersection""" 219 return set(self.keys()).intersection(list(s2.keys())) != set() 220 221 def union(self, s2, exps=None): 222 """Compute the union of self and s2, adjusting exps""" 223 res = self.copy() 224 tr = {} 225 for expr, var in s2.items(): 226 if expr in self: 227 if exps: 228 exps = exps.xreplace({var: res[expr]}) 229 tr[var] = res[expr] 230 else: 231 res[expr] = var 232 for var, rewr in s2.rewrites.items(): 233 res.rewrites[var] = rewr.xreplace(tr) 234 return res, exps 235 236 def copy(self): 237 """Create a shallow copy of SubsSet""" 238 r = SubsSet() 239 r.rewrites = self.rewrites.copy() 240 for expr, var in self.items(): 241 r[expr] = var 242 return r 243 244 245@debug 246def mrv(e, x): 247 """Returns a SubsSet of most rapidly varying (mrv) subexpressions of 'e', 248 and e rewritten in terms of these""" 249 e = powsimp(e, deep=True, combine='exp') 250 if not isinstance(e, Basic): 251 raise TypeError("e should be an instance of Basic") 252 if not e.has(x): 253 return SubsSet(), e 254 elif e == x: 255 s = SubsSet() 256 return s, s[x] 257 elif e.is_Mul or e.is_Add: 258 i, d = e.as_independent(x) # throw away x-independent terms 259 if d.func != e.func: 260 s, expr = mrv(d, x) 261 return s, e.func(i, expr) 262 a, b = d.as_two_terms() 263 s1, e1 = mrv(a, x) 264 s2, e2 = mrv(b, x) 265 return mrv_max1(s1, s2, e.func(i, e1, e2), x) 266 elif e.is_Pow and e.base != S.Exp1: 267 e1 = S.One 268 while e.is_Pow: 269 b1 = e.base 270 e1 *= e.exp 271 e = b1 272 if b1 == 1: 273 return SubsSet(), b1 274 if e1.has(x): 275 base_lim = limitinf(b1, x) 276 if base_lim is S.One: 277 return mrv(exp(e1 * (b1 - 1)), x) 278 return mrv(exp(e1 * log(b1)), x) 279 else: 280 s, expr = mrv(b1, x) 281 return s, expr**e1 282 elif isinstance(e, log): 283 s, expr = mrv(e.args[0], x) 284 return s, log(expr) 285 elif isinstance(e, exp) or (e.is_Pow and e.base == S.Exp1): 286 # We know from the theory of this algorithm that exp(log(...)) may always 287 # be simplified here, and doing so is vital for termination. 288 if isinstance(e.exp, log): 289 return mrv(e.exp.args[0], x) 290 # if a product has an infinite factor the result will be 291 # infinite if there is no zero, otherwise NaN; here, we 292 # consider the result infinite if any factor is infinite 293 li = limitinf(e.exp, x) 294 if any(_.is_infinite for _ in Mul.make_args(li)): 295 s1 = SubsSet() 296 e1 = s1[e] 297 s2, e2 = mrv(e.exp, x) 298 su = s1.union(s2)[0] 299 su.rewrites[e1] = exp(e2) 300 return mrv_max3(s1, e1, s2, exp(e2), su, e1, x) 301 else: 302 s, expr = mrv(e.exp, x) 303 return s, exp(expr) 304 elif e.is_Function: 305 l = [mrv(a, x) for a in e.args] 306 l2 = [s for (s, _) in l if s != SubsSet()] 307 if len(l2) != 1: 308 # e.g. something like BesselJ(x, x) 309 raise NotImplementedError("MRV set computation for functions in" 310 " several variables not implemented.") 311 s, ss = l2[0], SubsSet() 312 args = [ss.do_subs(x[1]) for x in l] 313 return s, e.func(*args) 314 elif e.is_Derivative: 315 raise NotImplementedError("MRV set computation for derviatives" 316 " not implemented yet.") 317 raise NotImplementedError( 318 "Don't know how to calculate the mrv of '%s'" % e) 319 320 321def mrv_max3(f, expsf, g, expsg, union, expsboth, x): 322 """ 323 Computes the maximum of two sets of expressions f and g, which 324 are in the same comparability class, i.e. max() compares (two elements of) 325 f and g and returns either (f, expsf) [if f is larger], (g, expsg) 326 [if g is larger] or (union, expsboth) [if f, g are of the same class]. 327 """ 328 if not isinstance(f, SubsSet): 329 raise TypeError("f should be an instance of SubsSet") 330 if not isinstance(g, SubsSet): 331 raise TypeError("g should be an instance of SubsSet") 332 if f == SubsSet(): 333 return g, expsg 334 elif g == SubsSet(): 335 return f, expsf 336 elif f.meets(g): 337 return union, expsboth 338 339 c = compare(list(f.keys())[0], list(g.keys())[0], x) 340 if c == ">": 341 return f, expsf 342 elif c == "<": 343 return g, expsg 344 else: 345 if c != "=": 346 raise ValueError("c should be =") 347 return union, expsboth 348 349 350def mrv_max1(f, g, exps, x): 351 """Computes the maximum of two sets of expressions f and g, which 352 are in the same comparability class, i.e. mrv_max1() compares (two elements of) 353 f and g and returns the set, which is in the higher comparability class 354 of the union of both, if they have the same order of variation. 355 Also returns exps, with the appropriate substitutions made. 356 """ 357 u, b = f.union(g, exps) 358 return mrv_max3(f, g.do_subs(exps), g, f.do_subs(exps), 359 u, b, x) 360 361 362@debug 363@cacheit 364@timeit 365def sign(e, x): 366 """ 367 Returns a sign of an expression e(x) for x->oo. 368 369 :: 370 371 e > 0 for x sufficiently large ... 1 372 e == 0 for x sufficiently large ... 0 373 e < 0 for x sufficiently large ... -1 374 375 The result of this function is currently undefined if e changes sign 376 arbitrarily often for arbitrarily large x (e.g. sin(x)). 377 378 Note that this returns zero only if e is *constantly* zero 379 for x sufficiently large. [If e is constant, of course, this is just 380 the same thing as the sign of e.] 381 """ 382 from sympy import sign as _sign 383 if not isinstance(e, Basic): 384 raise TypeError("e should be an instance of Basic") 385 386 if e.is_positive: 387 return 1 388 elif e.is_negative: 389 return -1 390 elif e.is_zero: 391 return 0 392 393 elif not e.has(x): 394 e = logcombine(e) 395 return _sign(e) 396 elif e == x: 397 return 1 398 elif e.is_Mul: 399 a, b = e.as_two_terms() 400 sa = sign(a, x) 401 if not sa: 402 return 0 403 return sa * sign(b, x) 404 elif isinstance(e, exp): 405 return 1 406 elif e.is_Pow: 407 if e.base == S.Exp1: 408 return 1 409 s = sign(e.base, x) 410 if s == 1: 411 return 1 412 if e.exp.is_Integer: 413 return s**e.exp 414 elif isinstance(e, log): 415 return sign(e.args[0] - 1, x) 416 417 # if all else fails, do it the hard way 418 c0, e0 = mrv_leadterm(e, x) 419 return sign(c0, x) 420 421 422@debug 423@timeit 424@cacheit 425def limitinf(e, x, leadsimp=False): 426 """Limit e(x) for x-> oo. 427 428 Explanation 429 =========== 430 431 If ``leadsimp`` is True, an attempt is made to simplify the leading 432 term of the series expansion of ``e``. That may succeed even if 433 ``e`` cannot be simplified. 434 """ 435 # rewrite e in terms of tractable functions only 436 437 if not e.has(x): 438 return e # e is a constant 439 if e.has(Order): 440 e = e.expand().removeO() 441 if not x.is_positive or x.is_integer: 442 # We make sure that x.is_positive is True and x.is_integer is None 443 # so we get all the correct mathematical behavior from the expression. 444 # We need a fresh variable. 445 p = Dummy('p', positive=True) 446 e = e.subs(x, p) 447 x = p 448 e = e.rewrite('tractable', deep=True, limitvar=x) 449 e = powdenest(e) 450 c0, e0 = mrv_leadterm(e, x) 451 sig = sign(e0, x) 452 if sig == 1: 453 return S.Zero # e0>0: lim f = 0 454 elif sig == -1: # e0<0: lim f = +-oo (the sign depends on the sign of c0) 455 if c0.match(I*Wild("a", exclude=[I])): 456 return c0*oo 457 s = sign(c0, x) 458 # the leading term shouldn't be 0: 459 if s == 0: 460 raise ValueError("Leading term should not be 0") 461 return s*oo 462 elif sig == 0: 463 if leadsimp: 464 c0 = c0.simplify() 465 return limitinf(c0, x, leadsimp) # e0=0: lim f = lim c0 466 else: 467 raise ValueError("{} could not be evaluated".format(sig)) 468 469 470def moveup2(s, x): 471 r = SubsSet() 472 for expr, var in s.items(): 473 r[expr.xreplace({x: exp(x)})] = var 474 for var, expr in s.rewrites.items(): 475 r.rewrites[var] = s.rewrites[var].xreplace({x: exp(x)}) 476 return r 477 478 479def moveup(l, x): 480 return [e.xreplace({x: exp(x)}) for e in l] 481 482 483@debug 484@timeit 485def calculate_series(e, x, logx=None): 486 """ Calculates at least one term of the series of ``e`` in ``x``. 487 488 This is a place that fails most often, so it is in its own function. 489 """ 490 from sympy.polys import cancel 491 from sympy.simplify import bottom_up 492 493 for t in e.lseries(x, logx=logx): 494 # bottom_up function is required for a specific case - when e is 495 # -exp(p/(p + 1)) + exp(-p**2/(p + 1) + p). No current simplification 496 # methods reduce this to 0 while not expanding polynomials. 497 t = bottom_up(t, lambda w: getattr(w, 'normal', lambda: w)()) 498 t = cancel(t, expand=False).factor() 499 500 if t.has(exp) and t.has(log): 501 t = powdenest(t) 502 503 if not t.is_zero: 504 break 505 506 return t 507 508 509@debug 510@timeit 511@cacheit 512def mrv_leadterm(e, x): 513 """Returns (c0, e0) for e.""" 514 Omega = SubsSet() 515 if not e.has(x): 516 return (e, S.Zero) 517 if Omega == SubsSet(): 518 Omega, exps = mrv(e, x) 519 if not Omega: 520 # e really does not depend on x after simplification 521 return exps, S.Zero 522 if x in Omega: 523 # move the whole omega up (exponentiate each term): 524 Omega_up = moveup2(Omega, x) 525 exps_up = moveup([exps], x)[0] 526 # NOTE: there is no need to move this down! 527 Omega = Omega_up 528 exps = exps_up 529 # 530 # The positive dummy, w, is used here so log(w*2) etc. will expand; 531 # a unique dummy is needed in this algorithm 532 # 533 # For limits of complex functions, the algorithm would have to be 534 # improved, or just find limits of Re and Im components separately. 535 # 536 w = Dummy("w", real=True, positive=True) 537 f, logw = rewrite(exps, Omega, x, w) 538 series = calculate_series(f, w, logx=logw) 539 try: 540 lt = series.leadterm(w, logx=logw) 541 except (ValueError, PoleError): 542 lt = f.as_coeff_exponent(w) 543 # as_coeff_exponent won't always split in required form. It may simply 544 # return (f, 0) when a better form may be obtained. Example (-x)**(-pi) 545 # can be written as (-1**(-pi), -pi) which as_coeff_exponent does not return 546 if lt[0].has(w): 547 base = f.as_base_exp()[0].as_coeff_exponent(w) 548 ex = f.as_base_exp()[1] 549 lt = (base[0]**ex, base[1]*ex) 550 return (lt[0].subs(log(w), logw), lt[1]) 551 552 553def build_expression_tree(Omega, rewrites): 554 r""" Helper function for rewrite. 555 556 We need to sort Omega (mrv set) so that we replace an expression before 557 we replace any expression in terms of which it has to be rewritten:: 558 559 e1 ---> e2 ---> e3 560 \ 561 -> e4 562 563 Here we can do e1, e2, e3, e4 or e1, e2, e4, e3. 564 To do this we assemble the nodes into a tree, and sort them by height. 565 566 This function builds the tree, rewrites then sorts the nodes. 567 """ 568 class Node: 569 def ht(self): 570 return reduce(lambda x, y: x + y, 571 [x.ht() for x in self.before], 1) 572 nodes = {} 573 for expr, v in Omega: 574 n = Node() 575 n.before = [] 576 n.var = v 577 n.expr = expr 578 nodes[v] = n 579 for _, v in Omega: 580 if v in rewrites: 581 n = nodes[v] 582 r = rewrites[v] 583 for _, v2 in Omega: 584 if r.has(v2): 585 n.before.append(nodes[v2]) 586 587 return nodes 588 589 590@debug 591@timeit 592def rewrite(e, Omega, x, wsym): 593 """e(x) ... the function 594 Omega ... the mrv set 595 wsym ... the symbol which is going to be used for w 596 597 Returns the rewritten e in terms of w and log(w). See test_rewrite1() 598 for examples and correct results. 599 """ 600 from sympy import ilcm 601 if not isinstance(Omega, SubsSet): 602 raise TypeError("Omega should be an instance of SubsSet") 603 if len(Omega) == 0: 604 raise ValueError("Length can not be 0") 605 # all items in Omega must be exponentials 606 for t in Omega.keys(): 607 if not isinstance(t, exp): 608 raise ValueError("Value should be exp") 609 rewrites = Omega.rewrites 610 Omega = list(Omega.items()) 611 612 nodes = build_expression_tree(Omega, rewrites) 613 Omega.sort(key=lambda x: nodes[x[1]].ht(), reverse=True) 614 615 # make sure we know the sign of each exp() term; after the loop, 616 # g is going to be the "w" - the simplest one in the mrv set 617 for g, _ in Omega: 618 sig = sign(g.exp, x) 619 if sig != 1 and sig != -1: 620 raise NotImplementedError('Result depends on the sign of %s' % sig) 621 if sig == 1: 622 wsym = 1/wsym # if g goes to oo, substitute 1/w 623 # O2 is a list, which results by rewriting each item in Omega using "w" 624 O2 = [] 625 denominators = [] 626 for f, var in Omega: 627 c = limitinf(f.exp/g.exp, x) 628 if c.is_Rational: 629 denominators.append(c.q) 630 arg = f.exp 631 if var in rewrites: 632 if not isinstance(rewrites[var], exp): 633 raise ValueError("Value should be exp") 634 arg = rewrites[var].args[0] 635 O2.append((var, exp((arg - c*g.exp).expand())*wsym**c)) 636 637 # Remember that Omega contains subexpressions of "e". So now we find 638 # them in "e" and substitute them for our rewriting, stored in O2 639 640 # the following powsimp is necessary to automatically combine exponentials, 641 # so that the .xreplace() below succeeds: 642 # TODO this should not be necessary 643 f = powsimp(e, deep=True, combine='exp') 644 for a, b in O2: 645 f = f.xreplace({a: b}) 646 647 for _, var in Omega: 648 assert not f.has(var) 649 650 # finally compute the logarithm of w (logw). 651 logw = g.exp 652 if sig == 1: 653 logw = -logw # log(w)->log(1/w)=-log(w) 654 655 # Some parts of sympy have difficulty computing series expansions with 656 # non-integral exponents. The following heuristic improves the situation: 657 exponent = reduce(ilcm, denominators, 1) 658 f = f.subs({wsym: wsym**exponent}) 659 logw /= exponent 660 661 return f, logw 662 663 664def gruntz(e, z, z0, dir="+"): 665 """ 666 Compute the limit of e(z) at the point z0 using the Gruntz algorithm. 667 668 Explanation 669 =========== 670 671 ``z0`` can be any expression, including oo and -oo. 672 673 For ``dir="+"`` (default) it calculates the limit from the right 674 (z->z0+) and for ``dir="-"`` the limit from the left (z->z0-). For infinite z0 675 (oo or -oo), the dir argument doesn't matter. 676 677 This algorithm is fully described in the module docstring in the gruntz.py 678 file. It relies heavily on the series expansion. Most frequently, gruntz() 679 is only used if the faster limit() function (which uses heuristics) fails. 680 """ 681 if not z.is_symbol: 682 raise NotImplementedError("Second argument must be a Symbol") 683 684 # convert all limits to the limit z->oo; sign of z is handled in limitinf 685 r = None 686 if z0 == oo: 687 e0 = e 688 elif z0 == -oo: 689 e0 = e.subs(z, -z) 690 else: 691 if str(dir) == "-": 692 e0 = e.subs(z, z0 - 1/z) 693 elif str(dir) == "+": 694 e0 = e.subs(z, z0 + 1/z) 695 else: 696 raise NotImplementedError("dir must be '+' or '-'") 697 698 try: 699 r = limitinf(e0, z) 700 except ValueError: 701 r = limitinf(e0, z, leadsimp=True) 702 703 # This is a bit of a heuristic for nice results... we always rewrite 704 # tractable functions in terms of familiar intractable ones. 705 # It might be nicer to rewrite the exactly to what they were initially, 706 # but that would take some work to implement. 707 return r.rewrite('intractable', deep=True) 708