1r""" 2The Gruntz Algorithm 3==================== 4 5This section explains the basics of the algorithm :cite:`Gruntz1996limits` used for computing 6limits. Most of the time the :py:func:`~diofant.series.limits.limit` function 7should just work. However it is still useful to keep in mind how it is 8implemented in case something does not work as expected. 9 10First we define an ordering on functions of single variable `x` according 11to how rapidly varying they at infinity. Any two functions `f(x)` and 12`g(x)` can be compared using the properties of: 13 14 .. math:: 15 L = \lim\limits_{x\to\infty}\frac{\log|f(x)|}{\log|g(x)|} 16 17We shall say that `f(x)` *dominates* `g(x)`, written `f(x) \succ 18g(x)`, iff `L=\pm\infty`. We also say that `f(x)` and `g(x)` are *of the 19same comparability class* if neither `f(x) \succ g(x)` nor `g(x) \succ 20f(x)` and shall denote it as `f(x) \asymp g(x)`. 21 22It is easy to show the following examples: 23 24* `e^{e^x} \succ e^{x^2} \succ e^x \succ x \succ 42` 25* `2 \asymp 3 \asymp -5` 26* `x \asymp x^2 \asymp x^3 \asymp -x` 27* `e^x \asymp e^{-x} \asymp e^{2x} \asymp e^{x + e^{-x}}` 28* `f(x) \asymp 1/f(x)` 29 30Using these definitions yields the following strategy for 31computing `\lim_{x \to \infty} f(x)`: 32 331. Given the function `f(x)`, we find the set of *most rapidly varying 34 subexpressions* (MRV set) of it. All items of this set belongs to the 35 same comparability class. Let's say it is `\{e^x, e^{2x}\}`. 36 372. Choose an expression `\omega` which is positive and tends to zero and 38 which is in the same comparability class as any element of the MRV set. 39 Such element always exists. Then we rewrite the MRV set using `\omega`, 40 in our case `\{\omega^{-1}, \omega^{-2}\}`, and substitute it into `f(x)`. 41 423. Let `f(\omega)` be the function which is obtained from `f(x)` after the 43 rewrite step above. Consider all expressions independent of `\omega` as 44 constants and compute the leading term of the power series of `f(\omega)` 45 around `\omega = 0^+`: 46 47 .. math:: f(\omega) = c_0 \omega^{e_0} + c_1 \omega^{e_1} + \ldots 48 49 where `e_0 < e_1 < e_2 \ldots` 50 514. If the leading exponent `e_0 > 0` then the limit is `0`. If `e_0 < 0`, 52 then the answer is `\pm\infty` (depends on sign of `c_0`). Finally, 53 if `e_0 = 0`, the limit is the limit of the leading coefficient `c_0`. 54 55Notes 56----- 57This exposition glossed over several details. For example, limits could be 58computed recursively (steps 1 and 4). Please address to the Gruntz thesis :cite:`Gruntz1996limits` 59for proof of the termination (pp. 52-60). 60 61""" 62 63import functools 64 65from ..core import Add, Dummy, E, Float, Integer, Mul, cacheit, oo 66from ..core.evaluate import evaluate 67from ..functions import Abs, exp, log, sign 68from ..utilities import ordered 69 70 71def compare(a, b, x): 72 r""" 73 Determine order relation between two functons. 74 75 Returns 76 ======= 77 78 {1, 0, -1} 79 Respectively, if `a(x) \succ b(x)`, `a(x) \asymp b(x)` 80 or `b(x) \succ a(x)`. 81 82 Examples 83 ======== 84 85 >>> x = Symbol('x', real=True, positive=True) 86 87 >>> compare(exp(x), x**5, x) 88 1 89 90 """ 91 # The log(exp(...)) must always be simplified here for termination. 92 la = a.exp if a.is_Pow and a.base is E else log(a) 93 lb = b.exp if b.is_Pow and b.base is E else log(b) 94 95 c = limitinf(la/lb, x) 96 if c.is_zero: 97 return -1 98 elif c.is_infinite: 99 return 1 100 else: 101 return 0 102 103 104def mrv(e, x): 105 """ 106 Calculate the MRV set of expression. 107 108 Examples 109 ======== 110 111 >>> x = Symbol('x', real=True, positive=True) 112 113 >>> mrv(log(x - log(x))/log(x), x) 114 {x} 115 116 """ 117 if not e.has(x): 118 return set() 119 elif e == x: 120 return {x} 121 elif e.is_Mul or e.is_Add: 122 a, b = e.as_two_terms() 123 return mrv_max(mrv(a, x), mrv(b, x), x) 124 elif e.is_Pow and e.base is E: 125 if e.exp == x: 126 return {e} 127 elif any(a.is_infinite for a in Mul.make_args(limitinf(e.exp, x))): 128 return mrv_max({e}, mrv(e.exp, x), x) 129 else: 130 return mrv(e.exp, x) 131 elif e.is_Pow: 132 assert not e.exp.has(x) 133 return mrv(e.base, x) 134 elif isinstance(e, log): 135 return mrv(e.args[0], x) 136 elif e.is_Function: 137 return functools.reduce(lambda a, b: mrv_max(a, b, x), 138 [mrv(a, x) for a in e.args]) 139 else: 140 raise NotImplementedError(f"Don't know how to calculate the mrv of '{e}'") 141 142 143def mrv_max(f, g, x): 144 """Computes the maximum of two MRV sets.""" 145 if not f or not g or f & g: 146 return f | g 147 148 c = compare(list(f)[0], list(g)[0], x) 149 if c: 150 return f if c > 0 else g 151 else: 152 return f | g 153 154 155@cacheit 156def signinf(e, x): 157 r""" 158 Determine a sign of an expression at infinity. 159 160 Returns 161 ======= 162 163 {1, 0, -1} 164 One or minus one, if `e > 0` or `e < 0` for `x` sufficiently 165 large and zero if `e` is *constantly* zero for `x\to\infty`. 166 167 The result of this function is currently undefined if `e` changes 168 sign arbitrarily often at infinity (e.g. `\sin(x)`). 169 170 """ 171 if not e.has(x): 172 return sign(e).simplify() 173 elif e == x: 174 return 1 175 elif e.is_Mul: 176 a, b = e.as_two_terms() 177 return signinf(a, x)*signinf(b, x) 178 elif e.is_Pow: 179 s = signinf(e.base, x) 180 if s == 1: 181 return 1 182 183 c0, _ = mrv_leadterm(e, x) 184 return signinf(c0, x) 185 186 187@cacheit 188def limitinf(e, x): 189 """ 190 Compute limit of the expression at the infinity. 191 192 Examples 193 ======== 194 195 >>> x = Symbol('x', real=True, positive=True) 196 197 >>> limitinf(exp(x)*(exp(1/x - exp(-x)) - exp(1/x)), x) 198 -1 199 200 """ 201 assert x.is_real and x.is_positive 202 assert not e.has(Float) 203 204 # Rewrite e in terms of tractable functions only: 205 e = e.rewrite('tractable', deep=True) 206 207 def transform_abs(f): 208 s = sign(limitinf(f.args[0], x)) 209 return s*f.args[0] if s in (1, -1) else f 210 211 e = e.replace(lambda f: isinstance(f, Abs) and f.has(x), 212 transform_abs) 213 214 if not e.has(x): 215 # This is a bit of a heuristic for nice results. We always rewrite 216 # tractable functions in terms of familiar intractable ones. 217 # TODO: It might be nicer to rewrite the exactly to what they were 218 # initially, but that would take some work to implement. 219 return e.rewrite('intractable', deep=True) 220 221 c0, e0 = mrv_leadterm(e, x) 222 sig = signinf(e0, x) 223 if sig == 1: 224 return Integer(0) 225 elif sig == -1: 226 s = signinf(c0, x) 227 assert s != 0 228 return s*oo 229 elif sig == 0: 230 return limitinf(c0, x) 231 else: 232 raise NotImplementedError(f'Result depends on the sign of {sig}') 233 234 235@cacheit 236def mrv_leadterm(e, x): 237 """ 238 Compute the leading term of the series. 239 240 Returns 241 ======= 242 243 tuple 244 The leading term `c_0 w^{e_0}` of the series of `e` in terms 245 of the most rapidly varying subexpression `w` in form of 246 the pair ``(c0, e0)`` of Expr. 247 248 Examples 249 ======== 250 251 >>> x = Symbol('x', real=True, positive=True) 252 253 >>> mrv_leadterm(1/exp(-x + exp(-x)) - exp(x), x) 254 (-1, 0) 255 256 """ 257 if not e.has(x): 258 return e, Integer(0) 259 260 e = e.replace(lambda f: f.is_Pow and f.exp.has(x), 261 lambda f: exp(log(f.base)*f.exp)) 262 e = e.replace(lambda f: f.is_Mul and sum(a.is_Pow for a in f.args) > 1, 263 lambda f: Mul(exp(Add(*[a.exp for a in f.args if a.is_Pow and a.base is E])), 264 *[a for a in f.args if not a.is_Pow or a.base is not E])) 265 266 # The positive dummy, w, is used here so log(w*2) etc. will expand. 267 # TODO: For limits of complex functions, the algorithm would have to 268 # be improved, or just find limits of Re and Im components separately. 269 w = Dummy('w', real=True, positive=True) 270 e, logw = rewrite(e, x, w) 271 272 lt = e.compute_leading_term(w, logx=logw) 273 c0, e0 = lt.as_coeff_exponent(w) 274 if c0.has(w): 275 raise NotImplementedError(f'Cannot compute mrv_leadterm({e}, {x}). ' 276 'The coefficient should have been free of ' 277 f'{w}, but got {c0}.') 278 return c0, e0 279 280 281def rewrite(e, x, w): 282 r""" 283 Rewrites expression in terms of the most rapidly varying subexpression. 284 285 Parameters 286 ========== 287 288 e : Expr 289 an expression 290 x : Symbol 291 variable of the `e` 292 w : Symbol 293 The symbol which is going to be used for substitution in place 294 of the most rapidly varying in `x` subexpression. 295 296 Returns 297 ======= 298 299 tuple 300 A pair: rewritten (in `w`) expression and `\log(w)`. 301 302 Examples 303 ======== 304 305 >>> x = Symbol('x', real=True, positive=True) 306 >>> m = Symbol('m', real=True, positive=True) 307 308 >>> rewrite(exp(x)*log(log(exp(x))), x, m) 309 (log(x)/m, -x) 310 311 """ 312 Omega = mrv(e, x) 313 if not Omega: 314 return e, None # e really does not depend on x 315 316 assert all(e.has(t) for t in Omega) 317 318 if x in Omega: 319 # Moving up in the asymptotical scale: 320 with evaluate(False): 321 e = e.xreplace({x: exp(x)}) 322 Omega = {s.xreplace({x: exp(x)}) for s in Omega} 323 324 Omega = list(ordered(Omega, keys=lambda a: -len(mrv(a, x)))) 325 326 for g in Omega: 327 sig = signinf(g.exp, x) 328 if sig not in (1, -1): 329 raise NotImplementedError(f'Result depends on the sign of {sig}') 330 331 if sig == 1: 332 w = 1/w # if g goes to oo, substitute 1/w 333 334 # Rewrite and substitute subexpressions in the Omega. 335 for a in Omega: 336 c = limitinf(a.exp/g.exp, x) 337 b = exp(a.exp - c*g.exp)*w**c # exponential must never be expanded here 338 with evaluate(False): 339 e = e.xreplace({a: b}) 340 341 return e, -sig*g.exp 342