Repository: sympy/sympy

Issue:
Intersection(S.Reals, Interval(0, oo), FiniteSet(n)) is not simplified
```python
In [10]: Intersection(S.Reals, Interval(0, oo), FiniteSet(n))
Out[10]: Intersection((-oo, oo), [0, oo), {n})
```

It should simplify to the following:
```python
In [10]: Intersection(S.Reals, Interval(0, oo), FiniteSet(n))
Out[10]: Intersection([0, oo), {n})
```

Discussion:
aktech: @hargup  @flacjacket 
leosartaj: This maybe done by implementing ``_intersect`` in ``Reals`` just like other sets?
gxyd: I am not sure if my comments matter, but I think the problem is not with the `S.Reals`  for ex. 
```
>>> Intersction(Interval(1,2), Interval(3,4), FiniteSet(n))
# returned as it is instead of EmptySet
```
As I see the code you are only apllying it to `S.Reals` so I guess my example would not be taken care of.

I think problem is in line no 1442 on sets.py, since if it finds any FiniteSet containing an element which is `None` for its containment. We should also make the `Intersection(star_here other_args)` there (I cannot find star on.my keypad) . This code may work on line 1442.
```
 other_args = list([intersection(star_here other_args)])
If s.emptyset in other_args:
     return s.emptyset
```
Please bare with me if you find anything wrong here.
aktech: @leosartaj 

> This maybe done by implementing _intersect in Reals just like other sets?

These are Global rules & all global rules should reside at one [place](https://github.com/sympy/sympy/blob/master/sympy/sets/sets.py#L1427).
aktech: @gxyd Thanks for the review. 

> I am not sure if my comments matter, but I think the problem is not with the S.Reals for ex.

Every comment matters, feel free to comment.

> As I see the code you are only apllying it to S.Reals so I guess my example would not be taken care of.

I see I am only fixing the subset of the problem you are talking about, So my fix isn't general, which is not apt. +1 for your code, it looks good.

> (I cannot find star on.my keypad)

Press Shift + 8


Patch:
diff --git a/sympy/polys/domains/pythonrational.py b/sympy/polys/domains/pythonrational.py
index fb455d829a..d3c6a39f55 100644
--- a/sympy/polys/domains/pythonrational.py
+++ b/sympy/polys/domains/pythonrational.py
@@ -8,8 +8,6 @@
 from sympy.polys.domains.domainelement import DomainElement
 
 from sympy.core.compatibility import integer_types
-from sympy.core.sympify import converter
-from sympy.core.numbers import Rational
 from sympy.printing.defaults import DefaultPrinting
 from sympy.utilities import public
 
@@ -246,8 +244,3 @@ def denom(self):
 
     numerator = numer
     denominator = denom
-
-
-def sympify_pythonrational(arg):
-    return Rational(arg.p, arg.q)
-converter[PythonRational] = sympify_pythonrational
diff --git a/sympy/polys/ring_series.py b/sympy/polys/ring_series.py
index bd93775834..74a6331df3 100644
--- a/sympy/polys/ring_series.py
+++ b/sympy/polys/ring_series.py
@@ -1,12 +1,10 @@
-from sympy.polys.domains import QQ, EX, ExpressionDomain
+from sympy.polys.domains import QQ
 from sympy.polys.rings import PolyElement, ring
 from sympy.polys.monomials import monomial_min, monomial_mul, monomial_div
-from sympy.polys.polyerrors import DomainError
 from mpmath.libmp.libintmath import ifac
 from sympy.core.numbers import Rational
 from sympy.core.compatibility import as_int, range
-from sympy.core import S, evaluate
-from sympy.functions import sin, cos, tan, atan, exp, atanh, tanh, log
+from sympy.core import S
 from mpmath.libmp.libintmath import giant_steps
 import math
 
@@ -34,8 +32,8 @@ def _invert_monoms(p1):
     terms = list(p1.items())
     terms.sort()
     deg = p1.degree()
-    R = p1.ring
-    p = R.zero
+    ring = p1.ring
+    p = ring.zero
     cv = p1.listcoeffs()
     mv = p1.listmonoms()
     for i in range(len(mv)):
@@ -70,9 +68,9 @@ def rs_trunc(p1, x, prec):
     >>> rs_trunc(p, x, 10)
     x**5 + x + 1
     """
-    R = p1.ring
-    p = R.zero
-    i = R.gens.index(x)
+    ring = p1.ring
+    p = ring.zero
+    i = ring.gens.index(x)
     for exp1 in p1:
         if exp1[i] >= prec:
             continue
@@ -97,18 +95,18 @@ def rs_mul(p1, p2, x, prec):
     >>> rs_mul(p1, p2, x, 3)
     3*x**2 + 3*x + 1
     """
-    R = p1.ring
-    p = R.zero
-    if R.__class__ != p2.ring.__class__ or R != p2.ring:
+    ring = p1.ring
+    p = ring.zero
+    if ring.__class__ != p2.ring.__class__ or ring != p2.ring:
         raise ValueError('p1 and p2 must have the same ring')
-    iv = R.gens.index(x)
+    iv = ring.gens.index(x)
     if not isinstance(p2, PolyElement):
         raise ValueError('p1 and p2 must have the same ring')
-    if R == p2.ring:
+    if ring == p2.ring:
         get = p.get
         items2 = list(p2.items())
         items2.sort(key=lambda e: e[0][iv])
-        if R.ngens == 1:
+        if ring.ngens == 1:
             for exp1, v1 in p1.items():
                 for exp2, v2 in items2:
                     exp = exp1[0] + exp2[0]
@@ -118,7 +116,7 @@ def rs_mul(p1, p2, x, prec):
                     else:
                         break
         else:
-            monomial_mul = R.monomial_mul
+            monomial_mul = ring.monomial_mul
             for exp1, v1 in p1.items():
                 for exp2, v2 in items2:
                     if exp1[iv] + exp2[iv] < prec:
@@ -145,13 +143,13 @@ def rs_square(p1, x, prec):
     >>> rs_square(p, x, 3)
     6*x**2 + 4*x + 1
     """
-    R = p1.ring
-    p = R.zero
-    iv = R.gens.index(x)
+    ring = p1.ring
+    p = ring.zero
+    iv = ring.gens.index(x)
     get = p.get
     items = list(p1.items())
     items.sort(key=lambda e: e[0][iv])
-    monomial_mul = R.monomial_mul
+    monomial_mul = ring.monomial_mul
     for i in range(len(items)):
         exp1, v1 = items[i]
         for j in range(i):
@@ -232,10 +230,10 @@ def _has_constant_term(p, x):
     >>> _has_constant_term(p, x)
     True
     """
-    R = p.ring
-    iv = R.gens.index(x)
-    zm = R.zero_monom
-    a = [0]*R.ngens
+    ring = p.ring
+    iv = ring.gens.index(x)
+    zm = ring.zero_monom
+    a = [0]*ring.ngens
     a[iv] = 1
     miv = tuple(a)
     for expv in p:
@@ -261,22 +259,17 @@ def _series_inversion1(p, x, prec):
     -x**3 + x**2 - x + 1
 
     """
-    R = p.ring
-    zm = R.zero_monom
-    c = p[zm]
+    ring = p.ring
+    zm = ring.zero_monom
     if zm not in p:
         raise ValueError('no constant term in series')
-    if _has_constant_term(p - c, x):
-        raise ValueError('p cannot contain a constant term depending on '
-            'parameters')
-    one = R(1)
-    if R.domain is EX:
-        one = 1
-    if c != one:
+    if _has_constant_term(p - p[zm], x):
+        raise ValueError('p cannot contain a constant term depending on parameters')
+    if p[zm] != ring(1):
         # TODO add check that it is a unit
-        p1 = R(1)/c
+        p1 = ring(1)/p[zm]
     else:
-        p1 = R(1)
+        p1 = ring(1)
     for precx in _giant_steps(prec):
         tmp = p1.square()
         tmp = rs_mul(tmp, p, x, precx)
@@ -299,9 +292,9 @@ def rs_series_inversion(p, x, prec):
     >>> rs_series_inversion(1 + x*y**2, y, 4)
     -x*y**2 + 1
     """
-    R = p.ring
-    zm = R.zero_monom
-    ii = R.gens.index(x)
+    ring = p.ring
+    zm = ring.zero_monom
+    ii = ring.gens.index(x)
     m = min(p, key=lambda k: k[ii])[ii]
     if m:
         raise NotImplementedError('no constant term in series')
@@ -424,10 +417,10 @@ def rs_series_from_list(p, c, x, prec, concur=1):
     sympy.polys.ring.compose
 
     """
-    R = p.ring
+    ring = p.ring
     n = len(c)
     if not concur:
-        q = R(1)
+        q = ring(1)
         s = c[0]*q
         for i in range(1, n):
             q = rs_mul(q, p, x, prec)
@@ -437,9 +430,9 @@ def rs_series_from_list(p, c, x, prec, concur=1):
     K, r = divmod(n, J)
     if r:
         K += 1
-    ax = [R(1)]
+    ax = [ring(1)]
     b = 1
-    q = R(1)
+    q = ring(1)
     if len(p) < 20:
         for i in range(1, J):
             q = rs_mul(q, p, x, prec)
@@ -453,8 +446,8 @@ def rs_series_from_list(p, c, x, prec, concur=1):
             ax.append(q)
     # optimize using rs_square
     pj = rs_mul(ax[-1], p, x, prec)
-    b = R(1)
-    s = R(0)
+    b = ring(1)
+    s = ring(0)
     for k in range(K - 1):
         r = J*k
         s1 = c[r]
@@ -468,7 +461,7 @@ def rs_series_from_list(p, c, x, prec, concur=1):
     k = K - 1
     r = J*k
     if r < n:
-        s1 = c[r]*R(1)
+        s1 = c[r]*ring(1)
         for j in range(1, J):
             if r + j >= n:
                 break
@@ -494,10 +487,10 @@ def rs_diff(p, x):
     >>> rs_diff(p, x)
     2*x*y**3 + 1
     """
-    R = p.ring
-    n = R.gens.index(x)
-    p1 = R.zero
-    mn = [0]*R.ngens
+    ring = p.ring
+    n = ring.gens.index(x)
+    p1 = ring.zero
+    mn = [0]*ring.ngens
     mn[n] = 1
     mn = tuple(mn)
     for expv in p:
@@ -521,10 +514,10 @@ def rs_integrate(p, x):
     >>> rs_integrate(p, x)
     1/3*x**3*y**3 + 1/2*x**2
     """
-    R = p.ring
-    p1 = R.zero
-    n = R.gens.index(x)
-    mn = [0]*R.ngens
+    ring = p.ring
+    p1 = ring.zero
+    n = ring.gens.index(x)
+    mn = [0]*ring.ngens
     mn[n] = 1
     mn = tuple(mn)
 
@@ -543,7 +536,7 @@ def fun(p, f, *args):
       args[-2] = iv: names of the series variables
       args[-1] = prec: list of the precisions of the series variables
 
-    The case with f method name is used to compute rs_tan and rs_nth_root
+    The case with f method name is used to compute tan and nth_root
     of a multivariate series:
 
       fun(p, tan, iv, prec)
@@ -562,11 +555,11 @@ def fun(p, f, *args):
     1/3*x**3*y**3 + 2*x**3*y**2 + x**3*y + 1/3*x**3 + x**2*y + x*y + x
 
     """
-    _R = p.ring
-    R1, _x = ring('_x', _R.domain)
+    _ring = p.ring
+    ring1, _x = ring('_x', _ring)
     h = int(args[-1])
     args1 = args[:-2] + (_x, h)
-    zm = _R.zero_monom
+    zm = _ring.zero_monom
     # separate the constant term of the series
     # compute the univariate series f(_x, .., 'x', sum(nv))
     # or _x.f(..., 'x', sum(nv)
@@ -594,8 +587,8 @@ def mul_xin(p, i, n):
     x_i is the ith variable in p
     """
     n = as_int(n)
-    R = p.ring
-    q = R(0)
+    ring = p.ring
+    q = ring(0)
     for k, v in p.items():
         k1 = list(k)
         k1[i] += n
@@ -621,34 +614,14 @@ def rs_log(p, x, prec):
     >>> rs_log(1 + x, x, 8)
     1/7*x**7 - 1/6*x**6 + 1/5*x**5 - 1/4*x**4 + 1/3*x**3 - 1/2*x**2 + x
     """
-    R = p.ring
+    ring = p.ring
     if p == 1:
         return 0
-    if _has_constant_term(p, x):
-        const = 0
-        zm = R.zero_monom
-        c = p[zm]
-        if c == 1:
-            pass
-        else:
-            c_expr = c.as_expr()
-            if R.domain is EX:
-                const = log(c_expr)
-            elif isinstance(c, PolyElement):
-                try:
-                    const = R(log(c_expr))
-                except ValueError:
-                    raise DomainError("The given series can't be expanded in this "
-                        "domain.")
-            else:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain")
-
-        dlog = p.diff(x)
-        dlog = rs_mul(dlog, _series_inversion1(p, x, prec), x, prec - 1)
-        return rs_integrate(dlog, x) + const
-    else:
-        raise NotImplementedError
+    if _has_constant_term(p - 1, x):
+        raise NotImplementedError('p - 1 must not have a constant term in the series variables')
+    dlog = p.diff(x)
+    dlog = rs_mul(dlog, _series_inversion1(p, x, prec), x, prec - 1)
+    return rs_integrate(dlog, x)
 
 def rs_LambertW(p, iv, prec):
     """
@@ -670,12 +643,12 @@ def rs_LambertW(p, iv, prec):
 
     LambertW
     """
-    R = p.ring
-    p1 = R(0)
+    ring = p.ring
+    p1 = ring(0)
     if _has_constant_term(p, iv):
         raise NotImplementedError('Polynomial must not have constant term in \
               the series variables')
-    if iv in R.gens:
+    if iv in ring.gens:
         for precx in _giant_steps(prec):
             e = rs_exp(p1, iv, precx)
             p2 = rs_mul(e, p1, iv, precx) - p
@@ -691,8 +664,8 @@ def _exp1(p, x, prec):
     """
     Helper function for ``rs_exp``
     """
-    R = p.ring
-    p1 = R(1)
+    ring = p.ring
+    p1 = ring(1)
     for precx in _giant_steps(prec):
         pt = p - rs_log(p1, x, precx)
         tmp = rs_mul(pt, p1, x, precx)
@@ -713,31 +686,12 @@ def rs_exp(p, x, prec):
     >>> rs_exp(x**2, x, 7)
     1/6*x**6 + 1/2*x**4 + x**2 + 1
     """
-    R = p.ring
+    ring = p.ring
     if _has_constant_term(p, x):
-        zm = R.zero_monom
-        c = p[zm]
-        c_expr = c.as_expr()
-        if R.domain is EX:
-            const = exp(c_expr)
-        elif isinstance(c, PolyElement):
-            try:
-                const = R(exp(c_expr))
-            except ValueError:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain.")
-        else:
-            raise DomainError("The given series can't be expanded in this "
-                "domain")
-        p1 = p - c
-
-    # Makes use of sympy fuctions to evaluate the values of the cos/sin
-    # of the constant term.
-        return const*rs_exp(p1, x, prec)
-
+        raise NotImplementedError
     if len(p) > 20:
         return _exp1(p, x, prec)
-    one = R(1)
+    one = ring(1)
     n = 1
     k = 1
     c = []
@@ -752,8 +706,8 @@ def rs_exp(p, x, prec):
 # TODO
 # Needs to be benchmarked before use in rs_atan
 def _atan_series(p, iv, prec):
-    R = p.ring
-    mo = R(-1)
+    ring = p.ring
+    mo = ring(-1)
     c = [-mo]
     p2 = rs_square(p, iv, prec)
     for k in range(1, prec):
@@ -783,32 +737,19 @@ def rs_atan(p, x, prec):
 
     atan
     """
-    R = p.ring
-    const = 0
     if _has_constant_term(p, x):
-        zm = R.zero_monom
-        c = p[zm]
-        c_expr = c.as_expr()
-        if R.domain is EX:
-            const = atan(c_expr)
-        elif isinstance(c, PolyElement):
-            try:
-                const = R(atan(c_expr))
-            except ValueError:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain.")
-        else:
-            raise DomainError("The given series can't be expanded in this "
-                "domain")
+        raise NotImplementedError('Polynomial must not have constant term in \
+              the series variables')
+    ring = p.ring
 
     # Instead of using a closed form formula, we differentiate atan(p) to get
     # `1/(1+p**2) * dp`, whose series expansion is much easier to calculate.
     # Finally we integrate to get back atan
     dp = p.diff(x)
-    p1 = rs_square(p, x, prec) + R(1)
+    p1 = rs_square(p, x, prec) + ring(1)
     p1 = rs_series_inversion(p1, x, prec - 1)
     p1 = rs_mul(dp, p1, x, prec - 1)
-    return rs_integrate(p1, x) + const
+    return rs_integrate(p1, x)
 
 def _tan1(p, x, prec):
     """
@@ -823,8 +764,8 @@ def _tan1(p, x, prec):
     Then `f(r) = 0`
     Or `y  = atan(x)` where `x = tan(y)` as required.
     """
-    R = p.ring
-    p1 = R(0)
+    ring = p.ring
+    p1 = ring(0)
     for precx in _giant_steps(prec):
         tmp = p - rs_atan(p1, x, precx)
         tmp = rs_mul(tmp, 1 + p1.square(), x, precx)
@@ -852,36 +793,13 @@ def rs_tan(p, x, prec):
 
    tan
    """
-    R = p.ring
-    const = 0
+    ring = p.ring
     if _has_constant_term(p, x):
-        zm = R.zero_monom
-        c = p[zm]
-        c_expr = c.as_expr()
-        if R.domain is EX:
-            const = tan(c_expr)
-        elif isinstance(c, PolyElement):
-            try:
-                const = R(tan(c_expr))
-            except ValueError:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain.")
-        else:
-            raise DomainError("The given series can't be expanded in this "
-                "domain")
-            raise NotImplementedError
-        p1 = p - c
-
-    # Makes use of sympy fuctions to evaluate the values of the cos/sin
-    # of the constant term.
-        t2 = rs_tan(p1, x, prec)
-        t = rs_series_inversion(1 - const*t2, x, prec)
-        return rs_mul(const + t2, t, x, prec)
-
-    if R.ngens == 1:
+        raise NotImplementedError('Polynomial must not have constant term in \
+              series variables')
+    if ring.ngens == 1:
         return _tan1(p, x, prec)
-    else:
-        return fun(p, _tan1, x, prec)
+    return fun(p, rs_tan, x, prec)
 
 def rs_sin(p, x, prec):
     """
@@ -904,38 +822,19 @@ def rs_sin(p, x, prec):
 
     sin
     """
-    R = x.ring
+    ring = x.ring
     if not p:
-        return R(0)
+        return ring(0)
+    # Support for constant term can be extended on the lines of rs_cos
     if _has_constant_term(p, x):
-        zm = R.zero_monom
-        c = p[zm]
-        c_expr = c.as_expr()
-        if R.domain is EX:
-            t1, t2 = sin(c_expr), cos(c_expr)
-        elif isinstance(c, PolyElement):
-            try:
-                t1, t2 = R(sin(c_expr)), R(cos(c_expr))
-            except ValueError:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain.")
-        else:
-            raise DomainError("The given series can't be expanded in this "
-                "domain")
-            raise NotImplementedError
-        p1 = p - c
-
-    # Makes use of sympy cos, sin fuctions to evaluate the values of the cos/sin
-    # of the constant term.
-        return rs_sin(p1, x, prec)*t2 + rs_cos(p1, x, prec)*t1
-
-    # Series is calculated in terms of tan as its evaluation is fast.
+        raise NotImplementedError
+    # Series is calcualed in terms of tan as its evaluation is fast.
     if len(p) > 20 and p.ngens == 1:
         t = rs_tan(p/2, x, prec)
         t2 = rs_square(t, x, prec)
         p1 = rs_series_inversion(1 + t2, x, prec)
         return rs_mul(p1, 2*t, x, prec)
-    one = R(1)
+    one = ring(1)
     n = 1
     c = [0]
     for k in range(2, prec + 2, 2):
@@ -965,36 +864,26 @@ def rs_cos(p, iv, prec):
 
     cos
     """
-    R = p.ring
+    ring = p.ring
     if _has_constant_term(p, iv):
-        zm = R.zero_monom
-        c = p[zm]
-        c_expr = c.as_expr()
-        if R.domain is EX:
-            t1, t2 = sin(c_expr), cos(c_expr)
-        elif isinstance(c, PolyElement):
-            try:
-                t1, t2 = R(sin(c_expr)), R(cos(c_expr))
-            except ValueError:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain.")
-        else:
-            raise DomainError("The given series can't be expanded in this "
-                "domain")
+        zm = ring.zero_monom
+        c = S(p[zm])
+        if not c.is_real:
             raise NotImplementedError
         p1 = p - c
 
     # Makes use of sympy cos, sin fuctions to evaluate the values of the cos/sin
-    # of the constant term.
-        return rs_cos(p1, iv, prec)*t2 - rs_sin(p1, iv, prec)*t1
+    # of the constant term. Should it be left unevaluated?
+        from sympy.functions import cos, sin
+        return cos(c)*rs_cos(p1, iv, prec) -  sin(c)*rs_sin(p1, iv, prec)
 
     # Series is calculated in terms of tan as its evaluation is fast.
-    if len(p) > 20 and R.ngens == 1:
+    if len(p) > 20 and ring.ngens == 1:
         t = rs_tan(p/2, iv, prec)
         t2 = rs_square(t, iv, prec)
         p1 = rs_series_inversion(1+t2, iv, prec)
         return rs_mul(p1 ,1 - t2, iv, prec)
-    one = R(1)
+    one = ring(1)
     n = 1
     c = []
     for k in range(2, prec + 2, 2):
@@ -1021,8 +910,8 @@ def check_series_var(p, iv):
 # TODO
 # Needs to be benchmarked before use in rs_atanh
 def _atanh(p, iv, prec):
-    R = p.ring
-    one = R(1)
+    ring = p.ring
+    one = ring(1)
     c = [one]
     p2 = rs_square(p, iv, prec)
     for k in range(1, prec):
@@ -1031,7 +920,7 @@ def _atanh(p, iv, prec):
     s = rs_mul(s, p, iv, prec)
     return s
 
-def rs_atanh(p, x, prec):
+def rs_atanh(p, iv, prec):
     """
     Hyperbolic arctangent of a series
 
@@ -1052,32 +941,19 @@ def rs_atanh(p, x, prec):
 
     atanh
     """
-    R = p.ring
-    const = 0
-    if _has_constant_term(p, x):
-        zm = R.zero_monom
-        c = p[zm]
-        c_expr = c.as_expr()
-        if R.domain is EX:
-            const = atanh(c_expr)
-        elif isinstance(c, PolyElement):
-            try:
-                const = R(atanh(c_expr))
-            except ValueError:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain.")
-        else:
-            raise DomainError("The given series can't be expanded in this "
-                "domain")
+    if _has_constant_term(p, iv):
+        raise NotImplementedError('Polynomial must not have constant term in \
+              the series variables')
+    ring = p.ring
 
     # Instead of using a closed form formula, we differentiate atanh(p) to get
     # `1/(1-p**2) * dp`, whose series expansion is much easier to calculate.
     # Finally we integrate to get back atanh
-    dp = rs_diff(p, x)
-    p1 = - rs_square(p, x, prec) + 1
-    p1 = rs_series_inversion(p1, x, prec - 1)
-    p1 = rs_mul(dp, p1, x, prec - 1)
-    return rs_integrate(p1, x) + const
+    dp = rs_diff(p, iv)
+    p1 = - rs_square(p, iv, prec) + 1
+    p1 = rs_series_inversion(p1, iv, prec - 1)
+    p1 = rs_mul(dp, p1, iv, prec - 1)
+    return rs_integrate(p1, iv)
 
 def rs_sinh(p, iv, prec):
     """
@@ -1144,15 +1020,15 @@ def _tanh(p, iv, prec):
 
     _tanh
     """
-    R = p.ring
-    p1 = R(0)
+    ring = p.ring
+    p1 = ring(0)
     for precx in _giant_steps(prec):
         tmp = p - rs_atanh(p1, iv, precx)
         tmp = rs_mul(tmp, 1 - p1.square(), iv, precx)
         p1 += tmp
     return p1
 
-def rs_tanh(p, x, prec):
+def rs_tanh(p, iv, prec):
     """
     Hyperbolic tangent of a series
 
@@ -1173,33 +1049,13 @@ def rs_tanh(p, x, prec):
 
     tanh
     """
-    R = p.ring
-    const = 0
-    if _has_constant_term(p, x):
-        zm = R.zero_monom
-        c = p[zm]
-        c_expr = c.as_expr()
-        if R.domain is EX:
-            const = tanh(c_expr)
-        elif isinstance(c, PolyElement):
-            try:
-                const = R(tanh(c_expr))
-            except ValueError:
-                raise DomainError("The given series can't be expanded in this "
-                    "domain.")
-        else:
-            raise DomainError("The given series can't be expanded in this "
-                "domain")
-            raise NotImplementedError
-        p1 = p - c
-        t1 = rs_tanh(p1, x, prec)
-        t = rs_series_inversion(1 + const*t1, x, prec)
-        return rs_mul(const + t1, t, x, prec)
-
-    if R.ngens == 1:
-        return _tanh(p, x, prec)
-    else:
-        return fun(p, _tanh, x, prec)
+    ring = p.ring
+    if _has_constant_term(p, iv):
+        raise NotImplementedError('Polynomial must not have constant term in \
+              the series variables')
+    if ring.ngens == 1:
+        return _tanh(p, iv, prec)
+    return fun(p, _tanh, iv, prec)
 
 def rs_newton(p, x, prec):
     """
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index cdbee1f317..1f4b4690ae 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -1439,7 +1439,10 @@ def reduce(args):
                 unk = [x for x in s
                        if any(other.contains(x) not in (True, False) for other in other_args)]
                 if unk:
-                    res += Intersection(*([s.func(*unk)] + other_args), evaluate=False)
+                    other_sets = Intersection(*other_args)
+                    if other_sets.is_EmptySet:
+                        return EmptySet()
+                    res += Intersection(s.func(*unk), other_sets, evaluate=False)
                 return res
 
         # If any of the sets are unions, return a Union of Intersections


Test patch:
diff --git a/sympy/polys/tests/test_ring_series.py b/sympy/polys/tests/test_ring_series.py
index e0c1b2e87c..0af6732caa 100644
--- a/sympy/polys/tests/test_ring_series.py
+++ b/sympy/polys/tests/test_ring_series.py
@@ -1,4 +1,4 @@
-from sympy.polys.domains import QQ, EX
+from sympy.polys.domains import QQ
 from sympy.polys.rings import ring
 from sympy.polys.ring_series import (_invert_monoms, rs_integrate,
   rs_trunc, rs_mul, rs_square, rs_pow, _has_constant_term, rs_series_inversion,
@@ -7,8 +7,6 @@
   rs_sinh, rs_cosh, rs_tanh, _tan1, fun, rs_series_reversion)
 from sympy.utilities.pytest import raises
 from sympy.core.compatibility import range
-from sympy.core.symbol import symbols
-from sympy.functions import sin, cos, exp, tan, atan, atanh, tanh, log
 
 def test_ring_series1():
     R, x = ring('x', QQ)
@@ -153,18 +151,6 @@ def test_log():
     assert p1 == (4*x**5*y**2 - 2*x**5*y - 2*x**4*y**2 + x**5/5 + 2*x**4*y -
                   x**4/4 - 2*x**3*y + x**3/3 + 2*x**2*y - x**2/2 + x)
 
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', EX)
-    assert rs_log(x + a, x, 5) == -EX(1/(4*a**4))*x**4 + EX(1/(3*a**3))*x**3 \
-        - EX(1/(2*a**2))*x**2 + EX(1/a)*x + EX(log(a))
-    assert rs_log(x + x**2*y + a, x, 4) == -EX(a**(-2))*x**3*y + \
-        EX(1/(3*a**3))*x**3 + EX(1/a)*x**2*y - EX(1/(2*a**2))*x**2 + \
-        EX(1/a)*x + EX(log(a))
-
-    p = x + x**2 + 3
-    assert rs_log(p, x, 10).compose(x, 5) == EX(log(3) + 19281291595/9920232)
-
 def test_exp():
     R, x = ring('x', QQ)
     p = x + x**4
@@ -180,23 +166,6 @@ def test_exp():
     p1 = rs_exp(p, x, prec)
     assert p1 == x + 1
 
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', QQ[exp(a), a])
-    assert rs_exp(x + a, x, 5) == exp(a)*x**4/24 + exp(a)*x**3/6 + \
-        exp(a)*x**2/2 + exp(a)*x + exp(a)
-    assert rs_exp(x + x**2*y + a, x, 5) == exp(a)*x**4*y**2/2 + \
-            exp(a)*x**4*y/2 + exp(a)*x**4/24 + exp(a)*x**3*y + \
-            exp(a)*x**3/6 + exp(a)*x**2*y + exp(a)*x**2/2 + exp(a)*x + exp(a)
-
-    R, x, y = ring('x, y', EX)
-    assert rs_exp(x + a, x, 5) ==  EX(exp(a)/24)*x**4 + EX(exp(a)/6)*x**3 + \
-        EX(exp(a)/2)*x**2 + EX(exp(a))*x + EX(exp(a))
-    assert rs_exp(x + x**2*y + a, x, 5) == EX(exp(a)/2)*x**4*y**2 + \
-        EX(exp(a)/2)*x**4*y + EX(exp(a)/24)*x**4 + EX(exp(a))*x**3*y + \
-        EX(exp(a)/6)*x**3 + EX(exp(a))*x**2*y + EX(exp(a)/2)*x**2 + \
-        EX(exp(a))*x + EX(exp(a))
-
 def test_newton():
     R, x = ring('x', QQ)
     p = x**2 - 2
@@ -223,19 +192,6 @@ def test_atan():
         2*x**7*y**9 - 1/7*x**7*y**7 - 1/3*x**6*y**9 + x**6*y**7 - x**5*y**7 + \
         1/5*x**5*y**5 - x**4*y**5 - 1/3*x**3*y**3 + x**2*y**3 + x*y
 
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', EX)
-    assert rs_atan(x + a, x, 5) == -EX((a**3 - a)/(a**8 + 4*a**6 + 6*a**4 + \
-        4*a**2 + 1))*x**4 + EX((3*a**2 - 1)/(3*a**6 + 9*a**4 + \
-        9*a**2 + 3))*x**3 - EX(a/(a**4 + 2*a**2 + 1))*x**2 + \
-        EX(1/(a**2 + 1))*x + EX(atan(a))
-    assert rs_atan(x + x**2*y + a, x, 4) == -EX(2*a/(a**4 + 2*a**2 + 1)) \
-        *x**3*y + EX((3*a**2 - 1)/(3*a**6 + 9*a**4 + 9*a**2 + 3))*x**3 + \
-        EX(1/(a**2 + 1))*x**2*y - EX(a/(a**4 + 2*a**2 + 1))*x**2 + EX(1/(a**2 \
-        + 1))*x + EX(atan(a))
-
-
 def test_tan():
     R, x, y = ring('x, y', QQ)
     assert rs_tan(x, x, 9) == \
@@ -244,29 +200,6 @@ def test_tan():
         4/3*x**7*y**9 + 17/315*x**7*y**7 + 1/3*x**6*y**9 + 2/3*x**6*y**7 + \
         x**5*y**7 + 2/15*x**5*y**5 + x**4*y**5 + 1/3*x**3*y**3 + x**2*y**3 + x*y
 
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', QQ[tan(a), a])
-    assert rs_tan(x + a, x, 5) == (tan(a)**5 + 5*tan(a)**3/3 + \
-        2*tan(a)/3)*x**4 + (tan(a)**4 + 4*tan(a)**2/3 + 1/3)*x**3 + \
-        (tan(a)**3 + tan(a))*x**2 + (tan(a)**2 + 1)*x + tan(a)
-    assert rs_tan(x + x**2*y + a, x, 4) == (2*tan(a)**3 + 2*tan(a))*x**3*y + \
-        (tan(a)**4 + 4/3*tan(a)**2 + 1/3)*x**3 + (tan(a)**2 + 1)*x**2*y + \
-        (tan(a)**3 + tan(a))*x**2 + (tan(a)**2 + 1)*x + tan(a)
-
-    R, x, y = ring('x, y', EX)
-    assert rs_tan(x + a, x, 5) == EX(tan(a)**5 + 5*tan(a)**3/3 + \
-        2*tan(a)/3)*x**4 + EX(tan(a)**4 + 4*tan(a)**2/3 + EX(1)/3)*x**3 + \
-        EX(tan(a)**3 + tan(a))*x**2 + EX(tan(a)**2 + 1)*x + EX(tan(a))
-    assert rs_tan(x + x**2*y + a, x, 4) ==  EX(2*tan(a)**3 + \
-        2*tan(a))*x**3*y + EX(tan(a)**4 + 4*tan(a)**2/3 + EX(1)/3)*x**3 + \
-        EX(tan(a)**2 + 1)*x**2*y + EX(tan(a)**3 + tan(a))*x**2 + \
-        EX(tan(a)**2 + 1)*x + EX(tan(a))
-
-    p = x + x**2 + 5
-    assert rs_atan(p, x, 10).compose(x, 10) == EX(atan(5) + 67701870330562640/ \
-        668083460499)
-
 def test_sin():
     R, x, y = ring('x, y', QQ)
     assert rs_sin(x, x, 9) == \
@@ -276,23 +209,6 @@ def test_sin():
         + 1/24*x**6*y**7 - 1/2*x**5*y**7 + 1/120*x**5*y**5 - 1/2*x**4*y**5 \
         - 1/6*x**3*y**3 + x**2*y**3 + x*y
 
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', QQ[sin(a), cos(a), a])
-    assert rs_sin(x + a, x, 5) == sin(a)*x**4/24 - cos(a)*x**3/6 - \
-        sin(a)*x**2/2 + cos(a)*x + sin(a)
-    assert rs_sin(x + x**2*y + a, x, 5) == -sin(a)*x**4*y**2/2 - \
-        cos(a)*x**4*y/2 + sin(a)*x**4/24 - sin(a)*x**3*y - cos(a)*x**3/6 + \
-        cos(a)*x**2*y - sin(a)*x**2/2 + cos(a)*x + sin(a)
-
-    R, x, y = ring('x, y', EX)
-    assert rs_sin(x + a, x, 5) == EX(sin(a)/24)*x**4 - EX(cos(a)/6)*x**3 - \
-        EX(sin(a)/2)*x**2 + EX(cos(a))*x + EX(sin(a))
-    assert rs_sin(x + x**2*y + a, x, 5) == -EX(sin(a)/2)*x**4*y**2 - \
-        EX(cos(a)/2)*x**4*y + EX(sin(a)/24)*x**4 - EX(sin(a))*x**3*y - \
-        EX(cos(a)/6)*x**3 + EX(cos(a))*x**2*y - EX(sin(a)/2)*x**2 + \
-        EX(cos(a))*x + EX(sin(a))
-
 def test_cos():
     R, x, y = ring('x, y', QQ)
     assert rs_cos(x, x, 9) == \
@@ -302,23 +218,6 @@ def test_cos():
         1/120*x**7*y**8 + 1/4*x**6*y**8 - 1/720*x**6*y**6 + 1/6*x**5*y**6 \
         - 1/2*x**4*y**6 + 1/24*x**4*y**4 - x**3*y**4 - 1/2*x**2*y**2 + 1
 
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', QQ[sin(a), cos(a), a])
-    assert rs_cos(x + a, x, 5) == cos(a)*x**4/24 + sin(a)*x**3/6 - \
-        cos(a)*x**2/2 - sin(a)*x + cos(a)
-    assert rs_cos(x + x**2*y + a, x, 5) == -cos(a)*x**4*y**2/2 + \
-        sin(a)*x**4*y/2 + cos(a)*x**4/24 - cos(a)*x**3*y + sin(a)*x**3/6 - \
-        sin(a)*x**2*y - cos(a)*x**2/2 - sin(a)*x + cos(a)
-
-    R, x, y = ring('x, y', EX)
-    assert rs_cos(x + a, x, 5) == EX(cos(a)/24)*x**4 + EX(sin(a)/6)*x**3 - \
-        EX(cos(a)/2)*x**2 - EX(sin(a))*x + EX(cos(a))
-    assert rs_cos(x + x**2*y + a, x, 5) == -EX(cos(a)/2)*x**4*y**2 + \
-        EX(sin(a)/2)*x**4*y + EX(cos(a)/24)*x**4 - EX(cos(a))*x**3*y + \
-        EX(sin(a)/6)*x**3 - EX(sin(a))*x**2*y - EX(cos(a)/2)*x**2 - \
-        EX(sin(a))*x + EX(cos(a))
-
 def test_cos_sin():
     R, x, y = ring('x, y', QQ)
     cos, sin = rs_cos_sin(x, x, 9)
@@ -335,22 +234,6 @@ def test_atanh():
         2*x**7*y**9 + 1/7*x**7*y**7 + 1/3*x**6*y**9 + x**6*y**7 + x**5*y**7 + \
         1/5*x**5*y**5 + x**4*y**5 + 1/3*x**3*y**3 + x**2*y**3 + x*y
 
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', EX)
-    assert rs_atanh(x + a, x, 5) == EX((a**3 + a)/(a**8 - 4*a**6 + 6*a**4 - \
-        4*a**2 + 1))*x**4 - EX((3*a**2 + 1)/(3*a**6 - 9*a**4 + \
-        9*a**2 - 3))*x**3 + EX(a/(a**4 - 2*a**2 + 1))*x**2 - EX(1/(a**2 - \
-        1))*x + EX(atanh(a))
-    assert rs_atanh(x + x**2*y + a, x, 4) == EX(2*a/(a**4 - 2*a**2 + \
-        1))*x**3*y - EX((3*a**2 + 1)/(3*a**6 - 9*a**4 + 9*a**2 - 3))*x**3 - \
-        EX(1/(a**2 - 1))*x**2*y + EX(a/(a**4 - 2*a**2 + 1))*x**2 - \
-        EX(1/(a**2 - 1))*x + EX(atanh(a))
-
-    p = x + x**2 + 5
-    assert rs_atanh(p, x, 10).compose(x, 10) == EX(-733442653682135/5079158784 \
-        + atanh(5))
-
 def test_sinh():
     R, x, y = ring('x, y', QQ)
     assert rs_sinh(x, x, 9) == 1/5040*x**7 + 1/120*x**5 + 1/6*x**3 + x
@@ -375,14 +258,3 @@ def test_tanh():
         17/45*x**8*y**9 + 4/3*x**7*y**9 - 17/315*x**7*y**7 - 1/3*x**6*y**9 + \
         2/3*x**6*y**7 - x**5*y**7 + 2/15*x**5*y**5 - x**4*y**5 - \
         1/3*x**3*y**3 + x**2*y**3 + x*y
-
-    # Constant term in series
-    a = symbols('a')
-    R, x, y = ring('x, y', EX)
-    assert rs_tanh(x + a, x, 5) == EX(tanh(a)**5 - 5*tanh(a)**3/3 + \
-        2*tanh(a)/3)*x**4 + EX(-tanh(a)**4 + 4*tanh(a)**2/3 - QQ(1, 3))*x**3 + \
-        EX(tanh(a)**3 - tanh(a))*x**2 + EX(-tanh(a)**2 + 1)*x + EX(tanh(a))
-
-    p = rs_tanh(x + x**2*y + a, x, 4)
-    assert (p.compose(x, 10)).compose(y, 5) == EX(-1000*tanh(a)**4 + \
-        10100*tanh(a)**3 + 2470*tanh(a)**2/3 - 10099*tanh(a) + QQ(530, 3))
diff --git a/sympy/sets/tests/test_sets.py b/sympy/sets/tests/test_sets.py
index 1ea8b30b37..fd8bc62c21 100644
--- a/sympy/sets/tests/test_sets.py
+++ b/sympy/sets/tests/test_sets.py
@@ -253,6 +253,7 @@ def test_intersect():
     assert Union(Interval(0, 1), Interval(2, 3)).intersection(Interval(1, 2)) == \
         Union(Interval(1, 1), Interval(2, 2))
 
+
 def test_intersection():
     # iterable
     i = Intersection(FiniteSet(1, 2, 3), Interval(2, 5), evaluate=False)
@@ -281,6 +282,16 @@ def test_intersection():
             Intersection(S.Integers, S.Naturals, S.Reals, evaluate=False)
 
 
+def test_issue_9623():
+    n = Symbol('n')
+
+    a = S.Reals
+    b = Interval(0, oo)
+    c = FiniteSet(n)
+
+    assert Intersection(a, b, c) == Intersection(b, c)
+    assert Intersection(Interval(1, 2), Interval(3, 4), FiniteSet(n)) == EmptySet()
+
 
 def test_is_disjoint():
     assert Interval(0, 2).is_disjoint(Interval(1, 2)) == False