intersection.py 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533
  1. from sympy.core.basic import _aresame
  2. from sympy.core.function import Lambda, expand_complex
  3. from sympy.core.mul import Mul
  4. from sympy.core.numbers import ilcm, Float
  5. from sympy.core.relational import Eq
  6. from sympy.core.singleton import S
  7. from sympy.core.symbol import (Dummy, symbols)
  8. from sympy.core.sorting import ordered
  9. from sympy.functions.elementary.complexes import sign
  10. from sympy.functions.elementary.integers import floor, ceiling
  11. from sympy.sets.fancysets import ComplexRegion
  12. from sympy.sets.sets import (FiniteSet, Intersection, Interval, Set, Union)
  13. from sympy.multipledispatch import Dispatcher
  14. from sympy.sets.conditionset import ConditionSet
  15. from sympy.sets.fancysets import (Integers, Naturals, Reals, Range,
  16. ImageSet, Rationals)
  17. from sympy.sets.sets import EmptySet, UniversalSet, imageset, ProductSet
  18. from sympy.simplify.radsimp import numer
  19. intersection_sets = Dispatcher('intersection_sets')
  20. @intersection_sets.register(ConditionSet, ConditionSet)
  21. def _(a, b):
  22. return None
  23. @intersection_sets.register(ConditionSet, Set)
  24. def _(a, b):
  25. return ConditionSet(a.sym, a.condition, Intersection(a.base_set, b))
  26. @intersection_sets.register(Naturals, Integers)
  27. def _(a, b):
  28. return a
  29. @intersection_sets.register(Naturals, Naturals)
  30. def _(a, b):
  31. return a if a is S.Naturals else b
  32. @intersection_sets.register(Interval, Naturals)
  33. def _(a, b):
  34. return intersection_sets(b, a)
  35. @intersection_sets.register(ComplexRegion, Set)
  36. def _(self, other):
  37. if other.is_ComplexRegion:
  38. # self in rectangular form
  39. if (not self.polar) and (not other.polar):
  40. return ComplexRegion(Intersection(self.sets, other.sets))
  41. # self in polar form
  42. elif self.polar and other.polar:
  43. r1, theta1 = self.a_interval, self.b_interval
  44. r2, theta2 = other.a_interval, other.b_interval
  45. new_r_interval = Intersection(r1, r2)
  46. new_theta_interval = Intersection(theta1, theta2)
  47. # 0 and 2*Pi means the same
  48. if ((2*S.Pi in theta1 and S.Zero in theta2) or
  49. (2*S.Pi in theta2 and S.Zero in theta1)):
  50. new_theta_interval = Union(new_theta_interval,
  51. FiniteSet(0))
  52. return ComplexRegion(new_r_interval*new_theta_interval,
  53. polar=True)
  54. if other.is_subset(S.Reals):
  55. new_interval = []
  56. x = symbols("x", cls=Dummy, real=True)
  57. # self in rectangular form
  58. if not self.polar:
  59. for element in self.psets:
  60. if S.Zero in element.args[1]:
  61. new_interval.append(element.args[0])
  62. new_interval = Union(*new_interval)
  63. return Intersection(new_interval, other)
  64. # self in polar form
  65. elif self.polar:
  66. for element in self.psets:
  67. if S.Zero in element.args[1]:
  68. new_interval.append(element.args[0])
  69. if S.Pi in element.args[1]:
  70. new_interval.append(ImageSet(Lambda(x, -x), element.args[0]))
  71. if S.Zero in element.args[0]:
  72. new_interval.append(FiniteSet(0))
  73. new_interval = Union(*new_interval)
  74. return Intersection(new_interval, other)
  75. @intersection_sets.register(Integers, Reals)
  76. def _(a, b):
  77. return a
  78. @intersection_sets.register(Range, Interval)
  79. def _(a, b):
  80. # Check that there are no symbolic arguments
  81. if not all(i.is_number for i in a.args + b.args[:2]):
  82. return
  83. # In case of null Range, return an EmptySet.
  84. if a.size == 0:
  85. return S.EmptySet
  86. # trim down to self's size, and represent
  87. # as a Range with step 1.
  88. start = ceiling(max(b.inf, a.inf))
  89. if start not in b:
  90. start += 1
  91. end = floor(min(b.sup, a.sup))
  92. if end not in b:
  93. end -= 1
  94. return intersection_sets(a, Range(start, end + 1))
  95. @intersection_sets.register(Range, Naturals)
  96. def _(a, b):
  97. return intersection_sets(a, Interval(b.inf, S.Infinity))
  98. @intersection_sets.register(Range, Range)
  99. def _(a, b):
  100. # Check that there are no symbolic range arguments
  101. if not all(all(v.is_number for v in r.args) for r in [a, b]):
  102. return None
  103. # non-overlap quick exits
  104. if not b:
  105. return S.EmptySet
  106. if not a:
  107. return S.EmptySet
  108. if b.sup < a.inf:
  109. return S.EmptySet
  110. if b.inf > a.sup:
  111. return S.EmptySet
  112. # work with finite end at the start
  113. r1 = a
  114. if r1.start.is_infinite:
  115. r1 = r1.reversed
  116. r2 = b
  117. if r2.start.is_infinite:
  118. r2 = r2.reversed
  119. # If both ends are infinite then it means that one Range is just the set
  120. # of all integers (the step must be 1).
  121. if r1.start.is_infinite:
  122. return b
  123. if r2.start.is_infinite:
  124. return a
  125. from sympy.solvers.diophantine.diophantine import diop_linear
  126. # this equation represents the values of the Range;
  127. # it's a linear equation
  128. eq = lambda r, i: r.start + i*r.step
  129. # we want to know when the two equations might
  130. # have integer solutions so we use the diophantine
  131. # solver
  132. va, vb = diop_linear(eq(r1, Dummy('a')) - eq(r2, Dummy('b')))
  133. # check for no solution
  134. no_solution = va is None and vb is None
  135. if no_solution:
  136. return S.EmptySet
  137. # there is a solution
  138. # -------------------
  139. # find the coincident point, c
  140. a0 = va.as_coeff_Add()[0]
  141. c = eq(r1, a0)
  142. # find the first point, if possible, in each range
  143. # since c may not be that point
  144. def _first_finite_point(r1, c):
  145. if c == r1.start:
  146. return c
  147. # st is the signed step we need to take to
  148. # get from c to r1.start
  149. st = sign(r1.start - c)*step
  150. # use Range to calculate the first point:
  151. # we want to get as close as possible to
  152. # r1.start; the Range will not be null since
  153. # it will at least contain c
  154. s1 = Range(c, r1.start + st, st)[-1]
  155. if s1 == r1.start:
  156. pass
  157. else:
  158. # if we didn't hit r1.start then, if the
  159. # sign of st didn't match the sign of r1.step
  160. # we are off by one and s1 is not in r1
  161. if sign(r1.step) != sign(st):
  162. s1 -= st
  163. if s1 not in r1:
  164. return
  165. return s1
  166. # calculate the step size of the new Range
  167. step = abs(ilcm(r1.step, r2.step))
  168. s1 = _first_finite_point(r1, c)
  169. if s1 is None:
  170. return S.EmptySet
  171. s2 = _first_finite_point(r2, c)
  172. if s2 is None:
  173. return S.EmptySet
  174. # replace the corresponding start or stop in
  175. # the original Ranges with these points; the
  176. # result must have at least one point since
  177. # we know that s1 and s2 are in the Ranges
  178. def _updated_range(r, first):
  179. st = sign(r.step)*step
  180. if r.start.is_finite:
  181. rv = Range(first, r.stop, st)
  182. else:
  183. rv = Range(r.start, first + st, st)
  184. return rv
  185. r1 = _updated_range(a, s1)
  186. r2 = _updated_range(b, s2)
  187. # work with them both in the increasing direction
  188. if sign(r1.step) < 0:
  189. r1 = r1.reversed
  190. if sign(r2.step) < 0:
  191. r2 = r2.reversed
  192. # return clipped Range with positive step; it
  193. # can't be empty at this point
  194. start = max(r1.start, r2.start)
  195. stop = min(r1.stop, r2.stop)
  196. return Range(start, stop, step)
  197. @intersection_sets.register(Range, Integers)
  198. def _(a, b):
  199. return a
  200. @intersection_sets.register(Range, Rationals)
  201. def _(a, b):
  202. return a
  203. @intersection_sets.register(ImageSet, Set)
  204. def _(self, other):
  205. from sympy.solvers.diophantine import diophantine
  206. # Only handle the straight-forward univariate case
  207. if (len(self.lamda.variables) > 1
  208. or self.lamda.signature != self.lamda.variables):
  209. return None
  210. base_set = self.base_sets[0]
  211. # Intersection between ImageSets with Integers as base set
  212. # For {f(n) : n in Integers} & {g(m) : m in Integers} we solve the
  213. # diophantine equations f(n)=g(m).
  214. # If the solutions for n are {h(t) : t in Integers} then we return
  215. # {f(h(t)) : t in integers}.
  216. # If the solutions for n are {n_1, n_2, ..., n_k} then we return
  217. # {f(n_i) : 1 <= i <= k}.
  218. if base_set is S.Integers:
  219. gm = None
  220. if isinstance(other, ImageSet) and other.base_sets == (S.Integers,):
  221. gm = other.lamda.expr
  222. var = other.lamda.variables[0]
  223. # Symbol of second ImageSet lambda must be distinct from first
  224. m = Dummy('m')
  225. gm = gm.subs(var, m)
  226. elif other is S.Integers:
  227. m = gm = Dummy('m')
  228. if gm is not None:
  229. fn = self.lamda.expr
  230. n = self.lamda.variables[0]
  231. try:
  232. solns = list(diophantine(fn - gm, syms=(n, m), permute=True))
  233. except (TypeError, NotImplementedError):
  234. # TypeError if equation not polynomial with rational coeff.
  235. # NotImplementedError if correct format but no solver.
  236. return
  237. # 3 cases are possible for solns:
  238. # - empty set,
  239. # - one or more parametric (infinite) solutions,
  240. # - a finite number of (non-parametric) solution couples.
  241. # Among those, there is one type of solution set that is
  242. # not helpful here: multiple parametric solutions.
  243. if len(solns) == 0:
  244. return S.EmptySet
  245. elif any(s.free_symbols for tupl in solns for s in tupl):
  246. if len(solns) == 1:
  247. soln, solm = solns[0]
  248. (t,) = soln.free_symbols
  249. expr = fn.subs(n, soln.subs(t, n)).expand()
  250. return imageset(Lambda(n, expr), S.Integers)
  251. else:
  252. return
  253. else:
  254. return FiniteSet(*(fn.subs(n, s[0]) for s in solns))
  255. if other == S.Reals:
  256. from sympy.solvers.solvers import denoms, solve_linear
  257. def _solution_union(exprs, sym):
  258. # return a union of linear solutions to i in expr;
  259. # if i cannot be solved, use a ConditionSet for solution
  260. sols = []
  261. for i in exprs:
  262. x, xis = solve_linear(i, 0, [sym])
  263. if x == sym:
  264. sols.append(FiniteSet(xis))
  265. else:
  266. sols.append(ConditionSet(sym, Eq(i, 0)))
  267. return Union(*sols)
  268. f = self.lamda.expr
  269. n = self.lamda.variables[0]
  270. n_ = Dummy(n.name, real=True)
  271. f_ = f.subs(n, n_)
  272. re, im = f_.as_real_imag()
  273. im = expand_complex(im)
  274. re = re.subs(n_, n)
  275. im = im.subs(n_, n)
  276. ifree = im.free_symbols
  277. lam = Lambda(n, re)
  278. if im.is_zero:
  279. # allow re-evaluation
  280. # of self in this case to make
  281. # the result canonical
  282. pass
  283. elif im.is_zero is False:
  284. return S.EmptySet
  285. elif ifree != {n}:
  286. return None
  287. else:
  288. # univarite imaginary part in same variable;
  289. # use numer instead of as_numer_denom to keep
  290. # this as fast as possible while still handling
  291. # simple cases
  292. base_set &= _solution_union(
  293. Mul.make_args(numer(im)), n)
  294. # exclude values that make denominators 0
  295. base_set -= _solution_union(denoms(f), n)
  296. return imageset(lam, base_set)
  297. elif isinstance(other, Interval):
  298. from sympy.solvers.solveset import (invert_real, invert_complex,
  299. solveset)
  300. f = self.lamda.expr
  301. n = self.lamda.variables[0]
  302. new_inf, new_sup = None, None
  303. new_lopen, new_ropen = other.left_open, other.right_open
  304. if f.is_real:
  305. inverter = invert_real
  306. else:
  307. inverter = invert_complex
  308. g1, h1 = inverter(f, other.inf, n)
  309. g2, h2 = inverter(f, other.sup, n)
  310. if all(isinstance(i, FiniteSet) for i in (h1, h2)):
  311. if g1 == n:
  312. if len(h1) == 1:
  313. new_inf = h1.args[0]
  314. if g2 == n:
  315. if len(h2) == 1:
  316. new_sup = h2.args[0]
  317. # TODO: Design a technique to handle multiple-inverse
  318. # functions
  319. # Any of the new boundary values cannot be determined
  320. if any(i is None for i in (new_sup, new_inf)):
  321. return
  322. range_set = S.EmptySet
  323. if all(i.is_real for i in (new_sup, new_inf)):
  324. # this assumes continuity of underlying function
  325. # however fixes the case when it is decreasing
  326. if new_inf > new_sup:
  327. new_inf, new_sup = new_sup, new_inf
  328. new_interval = Interval(new_inf, new_sup, new_lopen, new_ropen)
  329. range_set = base_set.intersect(new_interval)
  330. else:
  331. if other.is_subset(S.Reals):
  332. solutions = solveset(f, n, S.Reals)
  333. if not isinstance(range_set, (ImageSet, ConditionSet)):
  334. range_set = solutions.intersect(other)
  335. else:
  336. return
  337. if range_set is S.EmptySet:
  338. return S.EmptySet
  339. elif isinstance(range_set, Range) and range_set.size is not S.Infinity:
  340. range_set = FiniteSet(*list(range_set))
  341. if range_set is not None:
  342. return imageset(Lambda(n, f), range_set)
  343. return
  344. else:
  345. return
  346. @intersection_sets.register(ProductSet, ProductSet)
  347. def _(a, b):
  348. if len(b.args) != len(a.args):
  349. return S.EmptySet
  350. return ProductSet(*(i.intersect(j) for i, j in zip(a.sets, b.sets)))
  351. @intersection_sets.register(Interval, Interval)
  352. def _(a, b):
  353. # handle (-oo, oo)
  354. infty = S.NegativeInfinity, S.Infinity
  355. if a == Interval(*infty):
  356. l, r = a.left, a.right
  357. if l.is_real or l in infty or r.is_real or r in infty:
  358. return b
  359. # We can't intersect [0,3] with [x,6] -- we don't know if x>0 or x<0
  360. if not a._is_comparable(b):
  361. return None
  362. empty = False
  363. if a.start <= b.end and b.start <= a.end:
  364. # Get topology right.
  365. if a.start < b.start:
  366. start = b.start
  367. left_open = b.left_open
  368. elif a.start > b.start:
  369. start = a.start
  370. left_open = a.left_open
  371. else:
  372. start = a.start
  373. if not _aresame(a.start, b.start):
  374. # For example Integer(2) != Float(2)
  375. # Prefer the Float boundary because Floats should be
  376. # contagious in calculations.
  377. if b.start.has(Float) and not a.start.has(Float):
  378. start = b.start
  379. elif a.start.has(Float) and not b.start.has(Float):
  380. start = a.start
  381. else:
  382. #this is to ensure that if Eq(a.start, b.start) but
  383. #type(a.start) != type(b.start) the order of a and b
  384. #does not matter for the result
  385. start = list(ordered([a,b]))[0].start
  386. left_open = a.left_open or b.left_open
  387. if a.end < b.end:
  388. end = a.end
  389. right_open = a.right_open
  390. elif a.end > b.end:
  391. end = b.end
  392. right_open = b.right_open
  393. else:
  394. # see above for logic with start
  395. end = a.end
  396. if not _aresame(a.end, b.end):
  397. if b.end.has(Float) and not a.end.has(Float):
  398. end = b.end
  399. elif a.end.has(Float) and not b.end.has(Float):
  400. end = a.end
  401. else:
  402. end = list(ordered([a,b]))[0].end
  403. right_open = a.right_open or b.right_open
  404. if end - start == 0 and (left_open or right_open):
  405. empty = True
  406. else:
  407. empty = True
  408. if empty:
  409. return S.EmptySet
  410. return Interval(start, end, left_open, right_open)
  411. @intersection_sets.register(EmptySet, Set)
  412. def _(a, b):
  413. return S.EmptySet
  414. @intersection_sets.register(UniversalSet, Set)
  415. def _(a, b):
  416. return b
  417. @intersection_sets.register(FiniteSet, FiniteSet)
  418. def _(a, b):
  419. return FiniteSet(*(a._elements & b._elements))
  420. @intersection_sets.register(FiniteSet, Set)
  421. def _(a, b):
  422. try:
  423. return FiniteSet(*[el for el in a if el in b])
  424. except TypeError:
  425. return None # could not evaluate `el in b` due to symbolic ranges.
  426. @intersection_sets.register(Set, Set)
  427. def _(a, b):
  428. return None
  429. @intersection_sets.register(Integers, Rationals)
  430. def _(a, b):
  431. return a
  432. @intersection_sets.register(Naturals, Rationals)
  433. def _(a, b):
  434. return a
  435. @intersection_sets.register(Rationals, Reals)
  436. def _(a, b):
  437. return a
  438. def _intlike_interval(a, b):
  439. try:
  440. if b._inf is S.NegativeInfinity and b._sup is S.Infinity:
  441. return a
  442. s = Range(max(a.inf, ceiling(b.left)), floor(b.right) + 1)
  443. return intersection_sets(s, b) # take out endpoints if open interval
  444. except ValueError:
  445. return None
  446. @intersection_sets.register(Integers, Interval)
  447. def _(a, b):
  448. return _intlike_interval(a, b)
  449. @intersection_sets.register(Naturals, Interval)
  450. def _(a, b):
  451. return _intlike_interval(a, b)