/
limits.py
196 lines (153 loc) · 5.25 KB
/
limits.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
from ..core import (Add, Dummy, Expr, Float, Mul, PoleError, Rational, nan, oo,
sympify)
from ..core.function import UndefinedFunction
from ..sets import Reals
from .gruntz import limitinf
def limit(expr, z, z0, dir=None):
"""
Compute the directional limit of ``expr`` at the point ``z0``.
Examples
========
>>> limit(1/x, x, 0)
oo
>>> limit(1/x, x, 0, dir=1)
-oo
>>> limit(1/x, x, oo)
0
See Also
========
Limit
"""
return Limit(expr, z, z0, dir).doit(deep=False)
def heuristics(e, z, z0, dir):
from ..functions import cos, sin
e = e.expand()
if (e.is_Mul or e.is_Add or e.is_Pow or
(e.is_Function and not e.is_Piecewise and
not isinstance(e.func, UndefinedFunction))):
r = []
for a in e.args:
l = limit(a, z, z0, dir)
if (l.has(oo) and (l.func not in (sin, cos, Mul) and
l.is_finite is None)) or isinstance(l, Limit):
return
r.append(l)
if (rv := e.func(*r)) is not nan:
return rv
class Limit(Expr):
r"""
Represents an unevaluated directional limit of ``expr`` at the point ``z0``.
Parameters
==========
expr : Expr
algebraic expression
z : Symbol
variable of the ``expr``
z0 : Expr
limit point, `z_0`
dir : Expr or Reals, optional
selects the direction (as ``sign(dir)``) to approach the limit point
if the ``dir`` is an Expr. For infinite ``z0``, the default value
is determined from the direction of the infinity (e.g., the limit
from the left, ``dir=1``, for ``oo``). Otherwise, the default is
the limit from the right, ``dir=-1``. If ``dir=Reals``, the limit
is the bidirectional real limit.
Examples
========
>>> Limit(1/x, x, 0, dir=1)
Limit(1/x, x, 0, dir=1)
>>> _.doit()
-oo
See Also
========
limit
"""
def __new__(cls, e, z, z0, dir=None):
from ..functions import sign
e, z, z0, dir = map(sympify, [e, z, z0, dir])
if z0.is_infinite:
dir = sign(z0).simplify()
elif dir is None:
dir = Rational(-1)
if dir == Reals:
pass
elif isinstance(dir, Expr) and dir.is_nonzero:
dir = dir/abs(dir)
else:
raise ValueError('direction must be either a nonzero expression '
f'or Reals, not {dir}')
obj = super().__new__(cls)
obj._args = (e, z, z0, dir)
return obj
@property
def free_symbols(self):
e, z, z0 = self.args[:3]
return (e.free_symbols - z.free_symbols) | z0.free_symbols
def doit(self, **hints):
"""
Evaluates limit.
Notes
=====
First we handle some trivial cases (i.e. constant), then try
Gruntz algorithm (see the :py:mod:`~diofant.calculus.gruntz` module).
"""
e, z, z0, dir = self.args
if hints.get('deep', True):
e = e.doit(**hints)
z = z.doit(**hints)
z0 = z0.doit(**hints)
if dir == Reals:
right = limit(e, z, z0)
left = limit(e, z, z0, 1)
if not left.equals(right):
raise PoleError(f'left and right limits for the expression {e} '
f'at point {z}={z0} seems to be not equal')
return right
if has_Floats := e.has(Float) or z0.has(Float):
e = e.subs({k: Rational(k) for k in e.atoms(Float)})
z0 = z0.subs({k: Rational(k) for k in z0.atoms(Float)})
if z0.has(z):
newz = z.as_dummy()
r = limit(e.subs({z: newz}), newz, z0, dir)
if isinstance(r, Limit):
r = r.subs({newz: z})
return r
# Use a fresh variable to remove assumptions on the dummy variable.
newz = Dummy('z')
e, z = e.subs({z: newz}), newz
if not e.has(z):
return e
for t in e.atoms(Add):
if (o := t.getO()) and o.var == z and o.point == z0:
e = e.xreplace({t: t.removeO()})
if e.is_Order and e.var == z and e.point == z0:
if limit(e.expr, z, z0, dir):
return e
return Rational(0)
# Convert to the limit z->oo and use Gruntz algorithm.
e = e.subs({z: dir*z})
z0 = z0/dir
if z0 != oo:
e = e.subs({z: z0 - 1/z})
# We need a fresh variable with correct assumptions.
newz = Dummy(z.name, positive=True, finite=True)
e, z = e.subs({z: newz}), newz
if e.is_Boolean or e.is_Relational:
try:
has_oo = e.as_set().closure.contains(oo)
except NotImplementedError:
return self
if has_oo.is_Boolean:
return has_oo
raise NotImplementedError
try:
r = limitinf(e, z)
except (PoleError, ValueError, NotImplementedError):
r = None
if hints.get('heuristics', True):
r = heuristics(*self.args)
if r is None:
return self
if has_Floats:
r = r.evalf()
return r