limits.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394
  1. from sympy.calculus.accumulationbounds import AccumBounds
  2. from sympy.core import S, Symbol, Add, sympify, Expr, PoleError, Mul
  3. from sympy.core.exprtools import factor_terms
  4. from sympy.core.numbers import Float, _illegal
  5. from sympy.core.function import AppliedUndef
  6. from sympy.core.symbol import Dummy
  7. from sympy.functions.combinatorial.factorials import factorial
  8. from sympy.functions.elementary.complexes import (Abs, sign, arg, re)
  9. from sympy.functions.elementary.exponential import (exp, log)
  10. from sympy.functions.special.gamma_functions import gamma
  11. from sympy.polys import PolynomialError, factor
  12. from sympy.series.order import Order
  13. from .gruntz import gruntz
  14. def limit(e, z, z0, dir="+"):
  15. """Computes the limit of ``e(z)`` at the point ``z0``.
  16. Parameters
  17. ==========
  18. e : expression, the limit of which is to be taken
  19. z : symbol representing the variable in the limit.
  20. Other symbols are treated as constants. Multivariate limits
  21. are not supported.
  22. z0 : the value toward which ``z`` tends. Can be any expression,
  23. including ``oo`` and ``-oo``.
  24. dir : string, optional (default: "+")
  25. The limit is bi-directional if ``dir="+-"``, from the right
  26. (z->z0+) if ``dir="+"``, and from the left (z->z0-) if
  27. ``dir="-"``. For infinite ``z0`` (``oo`` or ``-oo``), the ``dir``
  28. argument is determined from the direction of the infinity
  29. (i.e., ``dir="-"`` for ``oo``).
  30. Examples
  31. ========
  32. >>> from sympy import limit, sin, oo
  33. >>> from sympy.abc import x
  34. >>> limit(sin(x)/x, x, 0)
  35. 1
  36. >>> limit(1/x, x, 0) # default dir='+'
  37. oo
  38. >>> limit(1/x, x, 0, dir="-")
  39. -oo
  40. >>> limit(1/x, x, 0, dir='+-')
  41. zoo
  42. >>> limit(1/x, x, oo)
  43. 0
  44. Notes
  45. =====
  46. First we try some heuristics for easy and frequent cases like "x", "1/x",
  47. "x**2" and similar, so that it's fast. For all other cases, we use the
  48. Gruntz algorithm (see the gruntz() function).
  49. See Also
  50. ========
  51. limit_seq : returns the limit of a sequence.
  52. """
  53. return Limit(e, z, z0, dir).doit(deep=False)
  54. def heuristics(e, z, z0, dir):
  55. """Computes the limit of an expression term-wise.
  56. Parameters are the same as for the ``limit`` function.
  57. Works with the arguments of expression ``e`` one by one, computing
  58. the limit of each and then combining the results. This approach
  59. works only for simple limits, but it is fast.
  60. """
  61. rv = None
  62. if z0 is S.Infinity:
  63. rv = limit(e.subs(z, 1/z), z, S.Zero, "+")
  64. if isinstance(rv, Limit):
  65. return
  66. elif (e.is_Mul or e.is_Add or e.is_Pow or (e.is_Function and not isinstance(e, AppliedUndef))):
  67. r = []
  68. from sympy.simplify.simplify import together
  69. for a in e.args:
  70. l = limit(a, z, z0, dir)
  71. if l.has(S.Infinity) and l.is_finite is None:
  72. if isinstance(e, Add):
  73. m = factor_terms(e)
  74. if not isinstance(m, Mul): # try together
  75. m = together(m)
  76. if not isinstance(m, Mul): # try factor if the previous methods failed
  77. m = factor(e)
  78. if isinstance(m, Mul):
  79. return heuristics(m, z, z0, dir)
  80. return
  81. return
  82. elif isinstance(l, Limit):
  83. return
  84. elif l is S.NaN:
  85. return
  86. else:
  87. r.append(l)
  88. if r:
  89. rv = e.func(*r)
  90. if rv is S.NaN and e.is_Mul and any(isinstance(rr, AccumBounds) for rr in r):
  91. r2 = []
  92. e2 = []
  93. for ii, rval in enumerate(r):
  94. if isinstance(rval, AccumBounds):
  95. r2.append(rval)
  96. else:
  97. e2.append(e.args[ii])
  98. if len(e2) > 0:
  99. e3 = Mul(*e2).simplify()
  100. l = limit(e3, z, z0, dir)
  101. rv = l * Mul(*r2)
  102. if rv is S.NaN:
  103. try:
  104. from sympy.simplify.ratsimp import ratsimp
  105. rat_e = ratsimp(e)
  106. except PolynomialError:
  107. return
  108. if rat_e is S.NaN or rat_e == e:
  109. return
  110. return limit(rat_e, z, z0, dir)
  111. return rv
  112. class Limit(Expr):
  113. """Represents an unevaluated limit.
  114. Examples
  115. ========
  116. >>> from sympy import Limit, sin
  117. >>> from sympy.abc import x
  118. >>> Limit(sin(x)/x, x, 0)
  119. Limit(sin(x)/x, x, 0, dir='+')
  120. >>> Limit(1/x, x, 0, dir="-")
  121. Limit(1/x, x, 0, dir='-')
  122. """
  123. def __new__(cls, e, z, z0, dir="+"):
  124. e = sympify(e)
  125. z = sympify(z)
  126. z0 = sympify(z0)
  127. if z0 in (S.Infinity, S.ImaginaryUnit*S.Infinity):
  128. dir = "-"
  129. elif z0 in (S.NegativeInfinity, S.ImaginaryUnit*S.NegativeInfinity):
  130. dir = "+"
  131. if(z0.has(z)):
  132. raise NotImplementedError("Limits approaching a variable point are"
  133. " not supported (%s -> %s)" % (z, z0))
  134. if isinstance(dir, str):
  135. dir = Symbol(dir)
  136. elif not isinstance(dir, Symbol):
  137. raise TypeError("direction must be of type basestring or "
  138. "Symbol, not %s" % type(dir))
  139. if str(dir) not in ('+', '-', '+-'):
  140. raise ValueError("direction must be one of '+', '-' "
  141. "or '+-', not %s" % dir)
  142. obj = Expr.__new__(cls)
  143. obj._args = (e, z, z0, dir)
  144. return obj
  145. @property
  146. def free_symbols(self):
  147. e = self.args[0]
  148. isyms = e.free_symbols
  149. isyms.difference_update(self.args[1].free_symbols)
  150. isyms.update(self.args[2].free_symbols)
  151. return isyms
  152. def pow_heuristics(self, e):
  153. _, z, z0, _ = self.args
  154. b1, e1 = e.base, e.exp
  155. if not b1.has(z):
  156. res = limit(e1*log(b1), z, z0)
  157. return exp(res)
  158. ex_lim = limit(e1, z, z0)
  159. base_lim = limit(b1, z, z0)
  160. if base_lim is S.One:
  161. if ex_lim in (S.Infinity, S.NegativeInfinity):
  162. res = limit(e1*(b1 - 1), z, z0)
  163. return exp(res)
  164. if base_lim is S.NegativeInfinity and ex_lim is S.Infinity:
  165. return S.ComplexInfinity
  166. def doit(self, **hints):
  167. """Evaluates the limit.
  168. Parameters
  169. ==========
  170. deep : bool, optional (default: True)
  171. Invoke the ``doit`` method of the expressions involved before
  172. taking the limit.
  173. hints : optional keyword arguments
  174. To be passed to ``doit`` methods; only used if deep is True.
  175. """
  176. e, z, z0, dir = self.args
  177. if str(dir) == '+-':
  178. r = limit(e, z, z0, dir='+')
  179. l = limit(e, z, z0, dir='-')
  180. if isinstance(r, Limit) and isinstance(l, Limit):
  181. if r.args[0] == l.args[0]:
  182. return self
  183. if r == l:
  184. return l
  185. if r.is_infinite and l.is_infinite:
  186. return S.ComplexInfinity
  187. raise ValueError("The limit does not exist since "
  188. "left hand limit = %s and right hand limit = %s"
  189. % (l, r))
  190. if z0 is S.ComplexInfinity:
  191. raise NotImplementedError("Limits at complex "
  192. "infinity are not implemented")
  193. if z0.is_infinite:
  194. cdir = sign(z0)
  195. cdir = cdir/abs(cdir)
  196. e = e.subs(z, cdir*z)
  197. dir = "-"
  198. z0 = S.Infinity
  199. if hints.get('deep', True):
  200. e = e.doit(**hints)
  201. z = z.doit(**hints)
  202. z0 = z0.doit(**hints)
  203. if e == z:
  204. return z0
  205. if not e.has(z):
  206. return e
  207. if z0 is S.NaN:
  208. return S.NaN
  209. if e.has(*_illegal):
  210. return self
  211. if e.is_Order:
  212. return Order(limit(e.expr, z, z0), *e.args[1:])
  213. cdir = S.Zero
  214. if str(dir) == "+":
  215. cdir = S.One
  216. elif str(dir) == "-":
  217. cdir = S.NegativeOne
  218. def set_signs(expr):
  219. if not expr.args:
  220. return expr
  221. newargs = tuple(set_signs(arg) for arg in expr.args)
  222. if newargs != expr.args:
  223. expr = expr.func(*newargs)
  224. abs_flag = isinstance(expr, Abs)
  225. arg_flag = isinstance(expr, arg)
  226. sign_flag = isinstance(expr, sign)
  227. if abs_flag or sign_flag or arg_flag:
  228. try:
  229. sig = limit(expr.args[0], z, z0, dir)
  230. if sig.is_zero:
  231. sig = limit(1/expr.args[0], z, z0, dir)
  232. except NotImplementedError:
  233. return expr
  234. else:
  235. if sig.is_extended_real:
  236. if (sig < 0) == True:
  237. return (-expr.args[0] if abs_flag else
  238. S.NegativeOne if sign_flag else S.Pi)
  239. elif (sig > 0) == True:
  240. return (expr.args[0] if abs_flag else
  241. S.One if sign_flag else S.Zero)
  242. return expr
  243. if e.has(Float):
  244. # Convert floats like 0.5 to exact SymPy numbers like S.Half, to
  245. # prevent rounding errors which can lead to unexpected execution
  246. # of conditional blocks that work on comparisons
  247. # Also see comments in https://github.com/sympy/sympy/issues/19453
  248. from sympy.simplify.simplify import nsimplify
  249. e = nsimplify(e)
  250. e = set_signs(e)
  251. if e.is_meromorphic(z, z0):
  252. if z0 is S.Infinity:
  253. newe = e.subs(z, 1/z)
  254. # cdir changes sign as oo- should become 0+
  255. cdir = -cdir
  256. else:
  257. newe = e.subs(z, z + z0)
  258. try:
  259. coeff, ex = newe.leadterm(z, cdir=cdir)
  260. except ValueError:
  261. pass
  262. else:
  263. if ex > 0:
  264. return S.Zero
  265. elif ex == 0:
  266. return coeff
  267. if cdir == 1 or not(int(ex) & 1):
  268. return S.Infinity*sign(coeff)
  269. elif cdir == -1:
  270. return S.NegativeInfinity*sign(coeff)
  271. else:
  272. return S.ComplexInfinity
  273. if z0 is S.Infinity:
  274. if e.is_Mul:
  275. e = factor_terms(e)
  276. dummy = Dummy('z', positive=z.is_positive, negative=z.is_negative, real=z.is_real)
  277. newe = e.subs(z, 1/dummy)
  278. # cdir changes sign as oo- should become 0+
  279. cdir = -cdir
  280. newz = dummy
  281. else:
  282. newe = e.subs(z, z + z0)
  283. newz = z
  284. try:
  285. coeff, ex = newe.leadterm(newz, cdir=cdir)
  286. except (ValueError, NotImplementedError, PoleError):
  287. # The NotImplementedError catching is for custom functions
  288. from sympy.simplify.powsimp import powsimp
  289. e = powsimp(e)
  290. if e.is_Pow:
  291. r = self.pow_heuristics(e)
  292. if r is not None:
  293. return r
  294. try:
  295. coeff = newe.as_leading_term(newz, cdir=cdir)
  296. if coeff != newe and (coeff.has(exp) or coeff.has(S.Exp1)):
  297. return gruntz(coeff, newz, 0, "-" if re(cdir).is_negative else "+")
  298. except (ValueError, NotImplementedError, PoleError):
  299. pass
  300. else:
  301. if isinstance(coeff, AccumBounds) and ex == S.Zero:
  302. return coeff
  303. if coeff.has(S.Infinity, S.NegativeInfinity, S.ComplexInfinity, S.NaN):
  304. return self
  305. if not coeff.has(newz):
  306. if ex.is_positive:
  307. return S.Zero
  308. elif ex == 0:
  309. return coeff
  310. elif ex.is_negative:
  311. if cdir == 1:
  312. return S.Infinity*sign(coeff)
  313. elif cdir == -1:
  314. return S.NegativeInfinity*sign(coeff)*S.NegativeOne**(S.One + ex)
  315. else:
  316. return S.ComplexInfinity
  317. else:
  318. raise NotImplementedError("Not sure of sign of %s" % ex)
  319. # gruntz fails on factorials but works with the gamma function
  320. # If no factorial term is present, e should remain unchanged.
  321. # factorial is defined to be zero for negative inputs (which
  322. # differs from gamma) so only rewrite for non-negative z0.
  323. if z0.is_extended_nonnegative:
  324. e = e.rewrite(factorial, gamma)
  325. l = None
  326. try:
  327. r = gruntz(e, z, z0, dir)
  328. if r is S.NaN or l is S.NaN:
  329. raise PoleError()
  330. except (PoleError, ValueError):
  331. if l is not None:
  332. raise
  333. r = heuristics(e, z, z0, dir)
  334. if r is None:
  335. return self
  336. return r