order.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507
  1. from sympy.core import S, sympify, Expr, Dummy, Add, Mul
  2. from sympy.core.cache import cacheit
  3. from sympy.core.containers import Tuple
  4. from sympy.core.function import Function, PoleError, expand_power_base, expand_log
  5. from sympy.core.sorting import default_sort_key
  6. from sympy.sets.sets import Complement
  7. from sympy.utilities.iterables import uniq, is_sequence
  8. class Order(Expr):
  9. r""" Represents the limiting behavior of some function.
  10. Explanation
  11. ===========
  12. The order of a function characterizes the function based on the limiting
  13. behavior of the function as it goes to some limit. Only taking the limit
  14. point to be a number is currently supported. This is expressed in
  15. big O notation [1]_.
  16. The formal definition for the order of a function `g(x)` about a point `a`
  17. is such that `g(x) = O(f(x))` as `x \rightarrow a` if and only if for any
  18. `\delta > 0` there exists a `M > 0` such that `|g(x)| \leq M|f(x)|` for
  19. `|x-a| < \delta`. This is equivalent to `\lim_{x \rightarrow a}
  20. \sup |g(x)/f(x)| < \infty`.
  21. Let's illustrate it on the following example by taking the expansion of
  22. `\sin(x)` about 0:
  23. .. math ::
  24. \sin(x) = x - x^3/3! + O(x^5)
  25. where in this case `O(x^5) = x^5/5! - x^7/7! + \cdots`. By the definition
  26. of `O`, for any `\delta > 0` there is an `M` such that:
  27. .. math ::
  28. |x^5/5! - x^7/7! + ....| <= M|x^5| \text{ for } |x| < \delta
  29. or by the alternate definition:
  30. .. math ::
  31. \lim_{x \rightarrow 0} | (x^5/5! - x^7/7! + ....) / x^5| < \infty
  32. which surely is true, because
  33. .. math ::
  34. \lim_{x \rightarrow 0} | (x^5/5! - x^7/7! + ....) / x^5| = 1/5!
  35. As it is usually used, the order of a function can be intuitively thought
  36. of representing all terms of powers greater than the one specified. For
  37. example, `O(x^3)` corresponds to any terms proportional to `x^3,
  38. x^4,\ldots` and any higher power. For a polynomial, this leaves terms
  39. proportional to `x^2`, `x` and constants.
  40. Examples
  41. ========
  42. >>> from sympy import O, oo, cos, pi
  43. >>> from sympy.abc import x, y
  44. >>> O(x + x**2)
  45. O(x)
  46. >>> O(x + x**2, (x, 0))
  47. O(x)
  48. >>> O(x + x**2, (x, oo))
  49. O(x**2, (x, oo))
  50. >>> O(1 + x*y)
  51. O(1, x, y)
  52. >>> O(1 + x*y, (x, 0), (y, 0))
  53. O(1, x, y)
  54. >>> O(1 + x*y, (x, oo), (y, oo))
  55. O(x*y, (x, oo), (y, oo))
  56. >>> O(1) in O(1, x)
  57. True
  58. >>> O(1, x) in O(1)
  59. False
  60. >>> O(x) in O(1, x)
  61. True
  62. >>> O(x**2) in O(x)
  63. True
  64. >>> O(x)*x
  65. O(x**2)
  66. >>> O(x) - O(x)
  67. O(x)
  68. >>> O(cos(x))
  69. O(1)
  70. >>> O(cos(x), (x, pi/2))
  71. O(x - pi/2, (x, pi/2))
  72. References
  73. ==========
  74. .. [1] `Big O notation <https://en.wikipedia.org/wiki/Big_O_notation>`_
  75. Notes
  76. =====
  77. In ``O(f(x), x)`` the expression ``f(x)`` is assumed to have a leading
  78. term. ``O(f(x), x)`` is automatically transformed to
  79. ``O(f(x).as_leading_term(x),x)``.
  80. ``O(expr*f(x), x)`` is ``O(f(x), x)``
  81. ``O(expr, x)`` is ``O(1)``
  82. ``O(0, x)`` is 0.
  83. Multivariate O is also supported:
  84. ``O(f(x, y), x, y)`` is transformed to
  85. ``O(f(x, y).as_leading_term(x,y).as_leading_term(y), x, y)``
  86. In the multivariate case, it is assumed the limits w.r.t. the various
  87. symbols commute.
  88. If no symbols are passed then all symbols in the expression are used
  89. and the limit point is assumed to be zero.
  90. """
  91. is_Order = True
  92. __slots__ = ()
  93. @cacheit
  94. def __new__(cls, expr, *args, **kwargs):
  95. expr = sympify(expr)
  96. if not args:
  97. if expr.is_Order:
  98. variables = expr.variables
  99. point = expr.point
  100. else:
  101. variables = list(expr.free_symbols)
  102. point = [S.Zero]*len(variables)
  103. else:
  104. args = list(args if is_sequence(args) else [args])
  105. variables, point = [], []
  106. if is_sequence(args[0]):
  107. for a in args:
  108. v, p = list(map(sympify, a))
  109. variables.append(v)
  110. point.append(p)
  111. else:
  112. variables = list(map(sympify, args))
  113. point = [S.Zero]*len(variables)
  114. if not all(v.is_symbol for v in variables):
  115. raise TypeError('Variables are not symbols, got %s' % variables)
  116. if len(list(uniq(variables))) != len(variables):
  117. raise ValueError('Variables are supposed to be unique symbols, got %s' % variables)
  118. if expr.is_Order:
  119. expr_vp = dict(expr.args[1:])
  120. new_vp = dict(expr_vp)
  121. vp = dict(zip(variables, point))
  122. for v, p in vp.items():
  123. if v in new_vp.keys():
  124. if p != new_vp[v]:
  125. raise NotImplementedError(
  126. "Mixing Order at different points is not supported.")
  127. else:
  128. new_vp[v] = p
  129. if set(expr_vp.keys()) == set(new_vp.keys()):
  130. return expr
  131. else:
  132. variables = list(new_vp.keys())
  133. point = [new_vp[v] for v in variables]
  134. if expr is S.NaN:
  135. return S.NaN
  136. if any(x in p.free_symbols for x in variables for p in point):
  137. raise ValueError('Got %s as a point.' % point)
  138. if variables:
  139. if any(p != point[0] for p in point):
  140. raise NotImplementedError(
  141. "Multivariable orders at different points are not supported.")
  142. if point[0] is S.Infinity:
  143. s = {k: 1/Dummy() for k in variables}
  144. rs = {1/v: 1/k for k, v in s.items()}
  145. ps = [S.Zero for p in point]
  146. elif point[0] is S.NegativeInfinity:
  147. s = {k: -1/Dummy() for k in variables}
  148. rs = {-1/v: -1/k for k, v in s.items()}
  149. ps = [S.Zero for p in point]
  150. elif point[0] is not S.Zero:
  151. s = {k: Dummy() + point[0] for k in variables}
  152. rs = {(v - point[0]).together(): k - point[0] for k, v in s.items()}
  153. ps = [S.Zero for p in point]
  154. else:
  155. s = ()
  156. rs = ()
  157. ps = list(point)
  158. expr = expr.subs(s)
  159. if expr.is_Add:
  160. expr = expr.factor()
  161. if s:
  162. args = tuple([r[0] for r in rs.items()])
  163. else:
  164. args = tuple(variables)
  165. if len(variables) > 1:
  166. # XXX: better way? We need this expand() to
  167. # workaround e.g: expr = x*(x + y).
  168. # (x*(x + y)).as_leading_term(x, y) currently returns
  169. # x*y (wrong order term!). That's why we want to deal with
  170. # expand()'ed expr (handled in "if expr.is_Add" branch below).
  171. expr = expr.expand()
  172. old_expr = None
  173. while old_expr != expr:
  174. old_expr = expr
  175. if expr.is_Add:
  176. lst = expr.extract_leading_order(args)
  177. expr = Add(*[f.expr for (e, f) in lst])
  178. elif expr:
  179. try:
  180. expr = expr.as_leading_term(*args)
  181. except PoleError:
  182. if isinstance(expr, Function) or\
  183. all(isinstance(arg, Function) for arg in expr.args):
  184. # It is not possible to simplify an expression
  185. # containing only functions (which raise error on
  186. # call to leading term) further
  187. pass
  188. else:
  189. orders = []
  190. pts = tuple(zip(args, ps))
  191. for arg in expr.args:
  192. try:
  193. lt = arg.as_leading_term(*args)
  194. except PoleError:
  195. lt = arg
  196. if lt not in args:
  197. order = Order(lt)
  198. else:
  199. order = Order(lt, *pts)
  200. orders.append(order)
  201. if expr.is_Add:
  202. new_expr = Order(Add(*orders), *pts)
  203. if new_expr.is_Add:
  204. new_expr = Order(Add(*[a.expr for a in new_expr.args]), *pts)
  205. expr = new_expr.expr
  206. elif expr.is_Mul:
  207. expr = Mul(*[a.expr for a in orders])
  208. elif expr.is_Pow:
  209. expr = orders[0].expr**orders[1].expr
  210. expr = expr.as_independent(*args, as_Add=False)[1]
  211. expr = expand_power_base(expr)
  212. expr = expand_log(expr)
  213. if len(args) == 1:
  214. # The definition of O(f(x)) symbol explicitly stated that
  215. # the argument of f(x) is irrelevant. That's why we can
  216. # combine some power exponents (only "on top" of the
  217. # expression tree for f(x)), e.g.:
  218. # x**p * (-x)**q -> x**(p+q) for real p, q.
  219. x = args[0]
  220. margs = list(Mul.make_args(
  221. expr.as_independent(x, as_Add=False)[1]))
  222. for i, t in enumerate(margs):
  223. if t.is_Pow:
  224. b, q = t.args
  225. if b in (x, -x) and q.is_real and not q.has(x):
  226. margs[i] = x**q
  227. elif b.is_Pow and not b.exp.has(x):
  228. b, r = b.args
  229. if b in (x, -x) and r.is_real:
  230. margs[i] = x**(r*q)
  231. elif b.is_Mul and b.args[0] is S.NegativeOne:
  232. b = -b
  233. if b.is_Pow and not b.exp.has(x):
  234. b, r = b.args
  235. if b in (x, -x) and r.is_real:
  236. margs[i] = x**(r*q)
  237. expr = Mul(*margs)
  238. expr = expr.subs(rs)
  239. if expr.is_Order:
  240. expr = expr.expr
  241. if not expr.has(*variables) and not expr.is_zero:
  242. expr = S.One
  243. # create Order instance:
  244. vp = dict(zip(variables, point))
  245. variables.sort(key=default_sort_key)
  246. point = [vp[v] for v in variables]
  247. args = (expr,) + Tuple(*zip(variables, point))
  248. obj = Expr.__new__(cls, *args)
  249. return obj
  250. def _eval_nseries(self, x, n, logx, cdir=0):
  251. return self
  252. @property
  253. def expr(self):
  254. return self.args[0]
  255. @property
  256. def variables(self):
  257. if self.args[1:]:
  258. return tuple(x[0] for x in self.args[1:])
  259. else:
  260. return ()
  261. @property
  262. def point(self):
  263. if self.args[1:]:
  264. return tuple(x[1] for x in self.args[1:])
  265. else:
  266. return ()
  267. @property
  268. def free_symbols(self):
  269. return self.expr.free_symbols | set(self.variables)
  270. def _eval_power(b, e):
  271. if e.is_Number and e.is_nonnegative:
  272. return b.func(b.expr ** e, *b.args[1:])
  273. if e == O(1):
  274. return b
  275. return
  276. def as_expr_variables(self, order_symbols):
  277. if order_symbols is None:
  278. order_symbols = self.args[1:]
  279. else:
  280. if (not all(o[1] == order_symbols[0][1] for o in order_symbols) and
  281. not all(p == self.point[0] for p in self.point)): # pragma: no cover
  282. raise NotImplementedError('Order at points other than 0 '
  283. 'or oo not supported, got %s as a point.' % self.point)
  284. if order_symbols and order_symbols[0][1] != self.point[0]:
  285. raise NotImplementedError(
  286. "Multiplying Order at different points is not supported.")
  287. order_symbols = dict(order_symbols)
  288. for s, p in dict(self.args[1:]).items():
  289. if s not in order_symbols.keys():
  290. order_symbols[s] = p
  291. order_symbols = sorted(order_symbols.items(), key=lambda x: default_sort_key(x[0]))
  292. return self.expr, tuple(order_symbols)
  293. def removeO(self):
  294. return S.Zero
  295. def getO(self):
  296. return self
  297. @cacheit
  298. def contains(self, expr):
  299. r"""
  300. Return True if expr belongs to Order(self.expr, \*self.variables).
  301. Return False if self belongs to expr.
  302. Return None if the inclusion relation cannot be determined
  303. (e.g. when self and expr have different symbols).
  304. """
  305. from sympy.simplify.powsimp import powsimp
  306. expr = sympify(expr)
  307. if expr.is_zero:
  308. return True
  309. if expr is S.NaN:
  310. return False
  311. point = self.point[0] if self.point else S.Zero
  312. if expr.is_Order:
  313. if (any(p != point for p in expr.point) or
  314. any(p != point for p in self.point)):
  315. return None
  316. if expr.expr == self.expr:
  317. # O(1) + O(1), O(1) + O(1, x), etc.
  318. return all(x in self.args[1:] for x in expr.args[1:])
  319. if expr.expr.is_Add:
  320. return all(self.contains(x) for x in expr.expr.args)
  321. if self.expr.is_Add and point.is_zero:
  322. return any(self.func(x, *self.args[1:]).contains(expr)
  323. for x in self.expr.args)
  324. if self.variables and expr.variables:
  325. common_symbols = tuple(
  326. [s for s in self.variables if s in expr.variables])
  327. elif self.variables:
  328. common_symbols = self.variables
  329. else:
  330. common_symbols = expr.variables
  331. if not common_symbols:
  332. return None
  333. if (self.expr.is_Pow and len(self.variables) == 1
  334. and self.variables == expr.variables):
  335. symbol = self.variables[0]
  336. other = expr.expr.as_independent(symbol, as_Add=False)[1]
  337. if (other.is_Pow and other.base == symbol and
  338. self.expr.base == symbol):
  339. if point.is_zero:
  340. rv = (self.expr.exp - other.exp).is_nonpositive
  341. if point.is_infinite:
  342. rv = (self.expr.exp - other.exp).is_nonnegative
  343. if rv is not None:
  344. return rv
  345. r = None
  346. ratio = self.expr/expr.expr
  347. ratio = powsimp(ratio, deep=True, combine='exp')
  348. for s in common_symbols:
  349. from sympy.series.limits import Limit
  350. l = Limit(ratio, s, point).doit(heuristics=False)
  351. if not isinstance(l, Limit):
  352. l = l != 0
  353. else:
  354. l = None
  355. if r is None:
  356. r = l
  357. else:
  358. if r != l:
  359. return
  360. return r
  361. if self.expr.is_Pow and len(self.variables) == 1:
  362. symbol = self.variables[0]
  363. other = expr.as_independent(symbol, as_Add=False)[1]
  364. if (other.is_Pow and other.base == symbol and
  365. self.expr.base == symbol):
  366. if point.is_zero:
  367. rv = (self.expr.exp - other.exp).is_nonpositive
  368. if point.is_infinite:
  369. rv = (self.expr.exp - other.exp).is_nonnegative
  370. if rv is not None:
  371. return rv
  372. obj = self.func(expr, *self.args[1:])
  373. return self.contains(obj)
  374. def __contains__(self, other):
  375. result = self.contains(other)
  376. if result is None:
  377. raise TypeError('contains did not evaluate to a bool')
  378. return result
  379. def _eval_subs(self, old, new):
  380. if old in self.variables:
  381. newexpr = self.expr.subs(old, new)
  382. i = self.variables.index(old)
  383. newvars = list(self.variables)
  384. newpt = list(self.point)
  385. if new.is_symbol:
  386. newvars[i] = new
  387. else:
  388. syms = new.free_symbols
  389. if len(syms) == 1 or old in syms:
  390. if old in syms:
  391. var = self.variables[i]
  392. else:
  393. var = syms.pop()
  394. # First, try to substitute self.point in the "new"
  395. # expr to see if this is a fixed point.
  396. # E.g. O(y).subs(y, sin(x))
  397. point = new.subs(var, self.point[i])
  398. if point != self.point[i]:
  399. from sympy.solvers.solveset import solveset
  400. d = Dummy()
  401. sol = solveset(old - new.subs(var, d), d)
  402. if isinstance(sol, Complement):
  403. e1 = sol.args[0]
  404. e2 = sol.args[1]
  405. sol = set(e1) - set(e2)
  406. res = [dict(zip((d, ), sol))]
  407. point = d.subs(res[0]).limit(old, self.point[i])
  408. newvars[i] = var
  409. newpt[i] = point
  410. elif old not in syms:
  411. del newvars[i], newpt[i]
  412. if not syms and new == self.point[i]:
  413. newvars.extend(syms)
  414. newpt.extend([S.Zero]*len(syms))
  415. else:
  416. return
  417. return Order(newexpr, *zip(newvars, newpt))
  418. def _eval_conjugate(self):
  419. expr = self.expr._eval_conjugate()
  420. if expr is not None:
  421. return self.func(expr, *self.args[1:])
  422. def _eval_derivative(self, x):
  423. return self.func(self.expr.diff(x), *self.args[1:]) or self
  424. def _eval_transpose(self):
  425. expr = self.expr._eval_transpose()
  426. if expr is not None:
  427. return self.func(expr, *self.args[1:])
  428. def __neg__(self):
  429. return self
  430. O = Order