Repository: sympy/sympy

Issue:
recursion error with solveset
```python
>>> solveset(f(x)-pi/2,x)
ConditionSet(x, Eq(2*f(x) - pi, 0), Complexes((-oo, oo) x (-oo, oo), False))
>>> solveset(f(x)-pi,x,S.Reals)
ConditionSet(x, Eq(f(x) - pi, 0), (-oo, oo))
>>> solveset(f(x)-pi/2,x,S.Reals)  # change pi to pi/2 and the error occurs unlike first case above
...
  File "sympy\solvers\solveset.py", line 513, in solveset_real
    result += solveset_real(equation, symbol)
  File "sympy\solvers\solveset.py", line 513, in solveset_real
    result += solveset_real(equation, symbol)
  File "sympy\solvers\solveset.py", line 513, in solveset_real
    result += solveset_real(equation, symbol)
  File "sympy\solvers\solveset.py", line 513, in solveset_real
    result += solveset_real(equation, symbol)
  File "sympy\solvers\solveset.py", line 476, in solveset_real
    elif f.is_Mul and all([_is_finite_with_finite_vars(m) for m in f.args]):
  File "sympy\solvers\solveset.py", line 327, in _is_finite_with_finite_vars
    return f.xreplace(reps).is_finite
  File "sympy\core\assumptions.py", line 245, in getit
    return _ask(fact, self)
  File "sympy\core\assumptions.py", line 288, in _ask
    a = evaluate(obj)
  File "sympy\core\add.py", line 457, in <lambda>
    (a.is_finite for a in self.args), quick_exit=True)
  File "sympy\core\logic.py", line 78, in _fuzzy_group
    for a in args:
  File "sympy\core\add.py", line 457, in <genexpr>
    (a.is_finite for a in self.args), quick_exit=True)
  File "sympy\core\assumptions.py", line 245, in getit
    return _ask(fact, self)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 288, in _ask
    a = evaluate(obj)
  File "sympy\core\mul.py", line 1245, in _eval_is_even
    is_integer = self.is_integer
  File "sympy\core\assumptions.py", line 245, in getit
    return _ask(fact, self)
  File "sympy\core\assumptions.py", line 288, in _ask
    a = evaluate(obj)
  File "sympy\core\mul.py", line 1054, in _eval_is_integer
    is_rational = self.is_rational
  File "sympy\core\assumptions.py", line 245, in getit
    return _ask(fact, self)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 300, in _ask
    _ask(pk, obj)
  File "sympy\core\assumptions.py", line 288, in _ask
    a = evaluate(obj)
  File "sympy\core\mul.py", line 1113, in _eval_is_hermitian
    return self._eval_herm_antiherm(True)
  File "sympy\core\mul.py", line 1119, in _eval_herm_antiherm
    if not t.is_commutative:
  File "sympy\core\function.py", line 453, in is_commutative
    if all(getattr(t, 'is_commutative') for t in self.args):
  File "sympy\core\function.py", line 453, in <genexpr>
    if all(getattr(t, 'is_commutative') for t in self.args):
  File "sympy\core\assumptions.py", line 241, in getit
    return self._assumptions[fact]
RuntimeError: maximum recursion depth exceeded
```

Discussion:
smichr: A strange error in mpmath has been uncovered that I will have to investigate.
asmeurer: Ping @hargup and @aktech on solveset issues. 
smichr: - [ ] check on invert_real and invert_complex's need to use domain (b/c in complex it wasn't used)
aktech: Ping me when it's ready to go.
smichr: ping @aktech. If you can take a look, that would be appreciated.
smichr: @aktech, can you give this a look?
smichr: It's been a while since I looked at this and having looked again, I agree with you: I will move the `_first` related material to the end of the helper.
smichr: @jksuom, do you have an idea about how to handle this latest issue? I am trying to do "safe simplification" by not allowing simplify to return an expression which has lost a singularity defined by some denominator, e.g. `1/(1/x+1)` is not allowed to simplify to `x/(1 + x)`. But then this disallows `factorial(x)/gamma(x)` from simplifying to 1. I wonder if I should only disallow denominators which are **Add** instances to not disappear.
smichr: If there are no more comments, @aktech , I will commit this in 12h.
aktech: @smichr I am sorry, for delaying this PR. Actually I am running short of time this week and next week so can't take a thorough look write now, so I would wait for @hargup to give a thumps up (he also probably have his exams write now). Sorry for the inconvenience.
jksuom: > 1/(1/x+1) is not allowed to simplify to x/(1 + x).

I'm afraid that trying to achieve this could result in problems hard to handle. In fact, there are reasons why I would rather allow this simplification.

BTW, you probably meant that ``factorial(x)/gamma(x)`` should simplify to ``x`` (for ``x`` a positive integer).
smichr: @jksuom , I added the comment to the wrong PR -- I'll move it to the other one regarding Relational simplification.
hargup: @smichr Thanks for this PR it greatly simplifies things. I have added my comments and after they are addressed it should be good to go.
smichr: Let's see if tests pass after the recent edit and if there are no final comments, I will commit this. Thanks for looking it over, @hargup.

Patch:
diff --git a/sympy/combinatorics/perm_groups.py b/sympy/combinatorics/perm_groups.py
index 37f07c0bc0..b1565fd52f 100644
--- a/sympy/combinatorics/perm_groups.py
+++ b/sympy/combinatorics/perm_groups.py
@@ -1306,14 +1306,14 @@ def generate_schreier_sims(self, af=False):
                     yield x._array_form
                 else:
                     yield x
-            return
+            raise StopIteration
         if len(u) == 1:
             for i in basic_orbits[0]:
                 if af:
                     yield u[0][i]._array_form
                 else:
                     yield u[0][i]
-            return
+            raise StopIteration
 
         u = list(reversed(u))
         basic_orbits = basic_orbits[::-1]
@@ -1327,7 +1327,7 @@ def generate_schreier_sims(self, af=False):
             # backtrack when finished iterating over coset
             if pos[h] >= posmax[h]:
                 if h == 0:
-                    return
+                    raise StopIteration
                 pos[h] = 0
                 h -= 1
                 stg.pop()
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 32954b2e7a..55ed6db91e 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -2622,7 +2622,7 @@ def _eval_lseries(self, x, logx=None):
                 yield series.removeO()
             else:
                 yield series
-            return
+            raise StopIteration
 
         while series.is_Order:
             n += 1
diff --git a/sympy/polys/subresultants_qq_zz.py b/sympy/polys/subresultants_qq_zz.py
index a8bc77be1c..1d71909588 100644
--- a/sympy/polys/subresultants_qq_zz.py
+++ b/sympy/polys/subresultants_qq_zz.py
@@ -1,214 +1,13 @@
 # -*- coding: utf-8 -*-
 """
-This module contains functions for the computation
-of Euclidean, generalized Sturmian and (modified) subresultant
-polynomial remainder sequences (prs's).
-
-The pseudo-remainder function prem() of sympy is _not_ used
-by any of the functions in the module.
-
-Instead of prem() we use the function
-
-rem_z().
-
-Included is also the function quo_z().
-
-1. Theoretical background:
-==========================
-Consider the polynomials f, g ∈ Z[x] of degrees deg(f) = n and
-deg(g) = m with n ≥ m.
-
-Definition 1:
-=============
-The sign sequence of a polynomial remainder sequence (prs) is the
-sequence of signs of the leading coefficients of its polynomials.
-
-Sign sequences can be computed with the function:
-
-sign_seq(poly_seq, x)
-
-Definition 2:
-=============
-A polynomial remainder sequence (prs) is called complete if the
-degree difference between any two consecutive polynomials is 1;
-otherwise, it called incomplete.
-
-It is understood that f, g belong to the sequences mentioned in
-the two definitions.
-
-1A. Euclidean and subresultant prs's:
-=====================================
-The subresultant prs of f, g is a sequence of polynomials in Z[x]
-analogous to the Euclidean prs, the sequence obtained by applying
-on f, g Euclid’s algorithm for polynomial greatest common divisors
-(gcd) in Q[x].
-
-The subresultant prs differs from the Euclidean prs in that the
-coefficients of each polynomial in the former sequence are determinants
---- also referred to as subresultants --- of appropriately selected
-sub-matrices of sylvester1(f, g, x), Sylvester’s matrix of 1840 of
-dimensions (n + m) × (n + m).
-
-Recall that the determinant of sylvester1(f, g, x) itself is
-called the resultant of f, g and serves as a criterion of whether
-the two polynomials have common roots or not.
-
-For complete prs’s the sign sequence of the Euclidean prs of f, g
-is identical to the sign sequence of the subresultant prs of f, g
-and the coefficients of one sequence  are easily computed from the
-coefficients of the  other.
-
-For incomplete prs’s the polynomials in the subresultant prs, generally
-differ in sign from those of the Euclidean prs, and --- unlike the
-case of complete prs’s --- it is not at all obvious how to compute
-the coefficients of one sequence from the coefficients of the  other.
-
-1B. Sturmian and modified subresultant prs's:
-=============================================
-For the same polynomials f, g ∈ Z[x] mentioned above, their ``modified''
-subresultant prs is a sequence of polynomials similar to the Sturmian
-prs, the sequence obtained by applying in Q[x] Sturm’s algorithm on f, g.
-
-The two sequences differ in that the coefficients of each polynomial
-in the modified subresultant prs are the determinants --- also referred
-to as modified subresultants --- of appropriately selected  sub-matrices
-of sylvester2(f, g, x), Sylvester’s matrix of 1853 of dimensions 2n × 2n.
-
-The determinant of sylvester2 itself is called the modified resultant
-of f, g and it also can serve as a criterion of whether the two
-polynomials have common roots or not.
-
-For complete prs’s the  sign sequence of the Sturmian prs of f, g is
-identical to the sign sequence of the modified subresultant prs of
-f, g and the coefficients of one sequence  are easily computed from
-the coefficients of the other.
-
-For incomplete prs’s the polynomials in the modified subresultant prs,
-generally differ in sign from those of the Sturmian prs, and --- unlike
-the case of complete prs’s --- it is not at all obvious how to compute
-the coefficients of one sequence from the coefficients of the  other.
-
-As Sylvester pointed out, the coefficients of the polynomial remainders
-obtained as (modified) subresultants are the smallest possible without
-introducing rationals and without computing (integer) greatest common
-divisors.
-
-1C. On terminology:
-===================
-Whence the terminology? Well generalized Sturmian prs's are
-``modifications'' of Euclidean prs's; the hint came from the title
-of the Pell-Gordon paper of 1917.
-
-In the literature one also encounters the name ``non signed'' and
-``signed'' prs for Euclidean and Sturmian prs respectively.
-
-Likewise ``non signed'' and ``signed'' subresultant prs for
-subresultant and modified subresultant prs respectively.
-
-2. Functions in the module:
-===========================
-No function utilizes sympy's function prem().
-
-2A. Matrices:
-=============
-The functions sylvester(f, g, x, method=1) and
-sylvester(f, g, x, method=2) compute either Sylvester matrix.
-They can be used to compute (modified) subresultant prs's by
-direct determinant evaluation.
-
-The function bezout(f, g, x, method='prs') provides a matrix of
-smaller dimensions than either Sylvester matrix. It is the function
-of choice for computing (modified) subresultant prs's by direct
-determinant evaluation.
-
-sylvester(f, g, x, method=1)
-sylvester(f, g, x, method=2)
-bezout(f, g, x, method='prs')
-
-The following identity holds:
-
-bezout(f, g, x, method='prs') =
-backward_eye(deg(f))*bezout(f, g, x, method='bz')*backward_eye(deg(f))
-
-2B. Subresultant and modified subresultant prs's by
-===================================================
-determinant evaluation:
-=======================
-Instead of utilizing the Sylvester matrices, we employ
-the Bezout matrix of smaller dimensions.
-
-subresultants_bezout(f, g, x)
-modified_subresultants_bezout(f, g, x)
-
-2C. Subresultant prs's by ONE determinant evaluation:
-=====================================================
-All three functions in this section evaluate one determinant
-per remainder polynomial; this is the determinant of an
-appropriately selected sub-matrix of sylvester1(f, g, x),
-Sylvester’s matrix of 1840.
-
-To compute the remainder polynomials the function
-subresultants_rem(f, g, x) employs rem(f, g, x).
-By contrast, the other two functions implement Van Vleck’s ideas
-of 1900 and compute the remainder polynomials by trinagularizing
-sylvester2(f, g, x), Sylvester’s matrix of 1853.
-
-
-subresultants_rem(f, g, x)
-subresultants_vv(f, g, x)
-subresultants_vv_2(f, g, x).
-
-2E. Euclidean, Sturmian prs's in Q[x]:
-======================================
-euclid_q(f, g, x)
-sturm_q(f, g, x)
-
-2F. Euclidean, Sturmian and (modified) subresultant prs's P-G:
-==============================================================
-All functions in this section are based on the Pell-Gordon (P-G)
-theorem of 1917.
-Computations are done in Q[x], employing the function rem(f, g, x)
-for the computation of the remainder polynomials.
-
-euclid_pg(f, g, x)
-sturm pg(f, g, x)
-subresultants_pg(f, g, x)
-modified_subresultants_pg(f, g, x)
-
-2G. Euclidean, Sturmian and (modified) subresultant prs's A-M-V:
-================================================================
-All functions in this section are based on the Akritas-Malaschonok-
-Vigklas (A-M-V) theorem of 2015.
-Computations are done in Z[x], employing the function rem_z(f, g, x)
-for the computation of the remainder polynomials.
-
-euclid_amv(f, g, x)
-sturm_amv(f, g, x)
-subresultants_amv(f, g, x)
-modified_subresultants_amv(f, g, x)
-
-2Ga. Exception:
-===============
-subresultants_amv_q(f, g, x)
-
-This function employs rem(f, g, x) for the computation of
-the remainder polynomials, despite the fact that it implements
-the A-M-V Theorem.
-
-It is included in our module in order to show that theorems P-G
-and A-M-V can be implemented utilizing either the function
-rem(f, g, x) or the function rem_z(f, g, x).
-
-For clearly historical reasons --- since the Collins-Brown-Traub
-coefficients-reduction factor β_i was not available in 1917 ---
-we have implemented the Pell-Gordon theorem with the function
-rem(f, g, x) and the A-M-V Theorem  with the function rem_z(f, g, x).
-"""
+Created on Mon Dec 28 13:25:02 2015
 
+@author: alkis
+"""
 
 from __future__ import print_function, division
 
-from sympy import (Abs, degree, expand, eye, floor, LC, Matrix, nan, Poly, pprint)
+from sympy import (Abs, degree, expand, floor, LC, Matrix, nan, Poly, pprint)
 from sympy import (QQ, quo, rem, S, sign, simplify, summation, var, zeros)
 
 def sylvester(f, g, x, method = 1):
@@ -324,251 +123,6 @@ def sign_seq(poly_seq, x):
     """
     return [sign(LC(poly_seq[i], x)) for i in range(len(poly_seq))]
 
-def bezout(p, q, x, method='bz'):
-    """
-    The input polynomials p, q are in Z[x] or in Q[x]. It is assumed
-    that degree(p, x) >= degree(q, x).
-
-    The default option bezout(p, q, x, method='bz') returns Bezout's
-    symmetric matrix of p and q, of dimensions deg(p) x deg(p). The
-    determinant of this matrix is equal to the determinant of sylvester2,
-    Sylvester's matrix of 1853, whose dimensions are 2*deg(p) x 2*deg(p);
-    however, the subresultants of these two matrices may vary in sign.
-
-    The other option, bezout(p, q, x, 'prs'), is of interest to us
-    in this module because it returns a matrix equivalent to sylvester2.
-    In this case all subresultants of the two matrices are identical.
-
-    Both the subresultant polynomial remainder sequence (prs) and
-    the modified subresultant prs of p and q can be computed by
-    evaluating determinants of appropriately selected submatrices of
-    bezout(p, q, x, 'prs') --- one determinant per coefficient of the
-    remainder polynomials.
-
-    The matrices bezout(p, q, x, 'bz') and bezout(p, q, x, 'prs')
-    are related by the formula
-
-    bezout(p, q, x, 'prs') =
-    backward_eye(deg(p)) * bezout(p, q, x, 'bz') * backward_eye(deg(p)),
-
-    where backward_eye() is the backward identity function.
-
-    References:
-    ===========
-    1. G.M.Diaz-Toca,L.Gonzalez-Vega: Various New Expressions for Subresultants
-    and Their Applications. Appl. Algebra in Engin., Communic. and Comp.,
-    Vol. 15, 233–266, 2004.
-
-    """
-    y = var('y')
-    degP = degree(p, x)
-
-    # expr is 0 when x = y
-    expr = p * q.subs({x:y}) - p.subs({x:y}) * q
-
-    # hence expr is exactly divisible by x - y
-    poly = Poly( quo(expr, x-y), x, y)
-
-    # form Bezout matrix and store them in B as indicated to get
-    # the LC coefficient of each poly in the first position of each row
-    B = zeros(degP)
-    for i in range(degP):
-        for j in range(degP):
-            if method == 'prs':
-                B[degP - 1 - i, degP - 1 - j] = poly.nth(i, j)
-            else:
-                B[i, j] = poly.nth(i, j)
-    return B
-
-def backward_eye(n):
-    '''
-    Returns the backward identity matrix of dimensions n x n.
-
-    Needed to "turn" the Bezout matrices
-    so that the leading coefficients are first.
-    See docstring of the function bezout(p, q, x, method='bz').
-    '''
-    M = eye(n)  # identity matrix of order n
-
-    for i in range(int(M.rows / 2)):
-        M.row_swap(0 + i, M.rows - 1 - i)
-
-    return M
-
-def process_bezout_output(poly_seq, x):
-    """
-    poly_seq is a polynomial remainder sequence computed either by
-    subresultants_bezout or by modified_subresultants_bezout.
-
-    This function removes from poly_seq all zero polynomials as well
-    as all those whose degree is equal to the degree of a previous
-    polynomial in poly_seq, as we scan it from left to right.
-
-    """
-    L = poly_seq[:]  # get a copy of the input sequence
-    d = degree(L[1], x)
-    i = 2
-    while i < len(L):
-        d_i = degree(L[i], x)
-        if d_i < 0:          # zero poly
-            L.remove(L[i])
-            i = i - 1
-        if d == d_i:         # poly degree equals degree of previous poly
-            L.remove(L[i])
-            i = i - 1
-        if d_i >= 0:
-            d = d_i
-        i = i + 1
-
-    return L
-
-def subresultants_bezout(p, q, x):
-    """
-    The input polynomials p, q are in Z[x] or in Q[x]. It is assumed
-    that degree(p, x) >= degree(q, x).
-
-    Computes the subresultant polynomial remainder sequence
-    of p, q by evaluating determinants of appropriately selected
-    submatrices of bezout(p, q, x, 'prs'). The dimensions of the
-    latter are deg(p) x deg(p).
-
-    Each coefficient is computed by evaluating the determinant of the
-    corresponding submatrix of bezout(p, q, x, 'prs').
-
-    bezout(p, q, x, 'prs) is used instead of sylvester(p, q, x, 1),
-    Sylvester's matrix of 1840, because the dimensions of the latter
-    are (deg(p) + deg(q)) x (deg(p) + deg(q)).
-
-    If the subresultant prs is complete, then the output coincides
-    with the Euclidean sequence of the polynomials p, q.
-
-    References:
-    ===========
-    1. G.M.Diaz-Toca,L.Gonzalez-Vega: Various New Expressions for Subresultants
-    and Their Applications. Appl. Algebra in Engin., Communic. and Comp.,
-    Vol. 15, 233–266, 2004.
-
-    """
-    # make sure neither p nor q is 0
-    if p == 0 or q == 0:
-        return [p, q]
-
-    f, g = p, q
-    n = degF = degree(f, x)
-    m = degG = degree(g, x)
-
-    # make sure proper degrees
-    if n == 0 and m == 0:
-        return [f, g]
-    if n < m:
-        n, m, degF, degG, f, g = m, n, degG, degF, g, f
-    if n > 0 and m == 0:
-        return [f, g]
-
-    SR_L = [f, g]      # subresultant list
-    F = LC(f, x)**(degF - degG)
-
-    # form the bezout matrix
-    B = bezout(f, g, x, 'prs')
-
-    # pick appropriate submatrices of B
-    # and form subresultant polys
-    if degF > degG:
-        j = 2
-    if degF == degG:
-        j = 1
-    while j <= degF:
-        M = B[0:j, :]
-        k, coeff_L = j - 1, []
-        while k <= degF - 1:
-            coeff_L.append(M[: ,0 : j].det())
-            if k < degF - 1:
-                M.col_swap(j - 1, k + 1)
-            k = k + 1
-
-        # apply Theorem 2.1 in the paper by Toca & Vega 2004
-        # to get correct signs
-        SR_L.append((int((-1)**(j*(j-1)/2)) * Poly(coeff_L, x) / F).as_expr())
-        j = j + 1
-
-    return process_bezout_output(SR_L, x)
-
-def modified_subresultants_bezout(p, q, x):
-    """
-    The input polynomials p, q are in Z[x] or in Q[x]. It is assumed
-    that degree(p, x) >= degree(q, x).
-
-    Computes the modified subresultant polynomial remainder sequence
-    of p, q by evaluating determinants of appropriately selected
-    submatrices of bezout(p, q, x, 'prs'). The dimensions of the
-    latter are deg(p) x deg(p).
-
-    Each coefficient is computed by evaluating the determinant of the
-    corresponding submatrix of bezout(p, q, x, 'prs').
-
-    bezout(p, q, x, 'prs') is used instead of sylvester(p, q, x, 2),
-    Sylvester's matrix of 1853, because the dimensions of the latter
-    are 2*deg(p) x 2*deg(p).
-
-    If the modified subresultant prs is complete, and LC( p ) > 0, the output
-    coincides with the (generalized) Sturm's sequence of the polynomials p, q.
-
-    References:
-    ===========
-    1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
-    and Modified Subresultant Polynomial Remainder Sequences.''
-    Serdica Journal of Computing, Vol. 8, No 1, 29–46, 2014.
-
-    2. G.M.Diaz-Toca,L.Gonzalez-Vega: Various New Expressions for Subresultants
-    and Their Applications. Appl. Algebra in Engin., Communic. and Comp.,
-    Vol. 15, 233–266, 2004.
-
-
-    """
-    # make sure neither p nor q is 0
-    if p == 0 or q == 0:
-        return [p, q]
-
-    f, g = p, q
-    n = degF = degree(f, x)
-    m = degG = degree(g, x)
-
-    # make sure proper degrees
-    if n == 0 and m == 0:
-        return [f, g]
-    if n < m:
-        n, m, degF, degG, f, g = m, n, degG, degF, g, f
-    if n > 0 and m == 0:
-        return [f, g]
-
-    SR_L = [f, g]      # subresultant list
-
-    # form the bezout matrix
-    B = bezout(f, g, x, 'prs')
-
-    # pick appropriate submatrices of B
-    # and form subresultant polys
-    if degF > degG:
-        j = 2
-    if degF == degG:
-        j = 1
-    while j <= degF:
-        M = B[0:j, :]
-        k, coeff_L = j - 1, []
-        while k <= degF - 1:
-            coeff_L.append(M[: ,0 : j].det())
-            if k < degF - 1:
-                M.col_swap(j - 1, k + 1)
-            k = k + 1
-
-        ## Theorem 2.1 in the paper by Toca & Vega 2004 is _not needed_
-        ## in this case since
-        ## the bezout matrix is equivalent to sylvester2
-        SR_L.append(( Poly(coeff_L, x)).as_expr())
-        j = j + 1
-
-    return process_bezout_output(SR_L, x)
-
 def sturm_pg(p, q, x, method=0):
     """
     p, q are polynomials in Z[x] or Q[x].
diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index aa5794e9e2..492b6b6f9a 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -1784,9 +1784,11 @@ def check_type(x, y):
         return 'type3'
     r1 = eq[0].match(diff(x(t),t) - f)
     r2 = eq[1].match(diff(y(t),t) - g)
-    num, denum = ((r1[f].subs(x(t),u).subs(y(t),v))/(r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom()
+    num, den = (
+        (r1[f].subs(x(t),u).subs(y(t),v))/
+        (r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom()
     R1 = num.match(f1*g1)
-    R2 = denum.match(f2*g2)
+    R2 = den.match(f2*g2)
     phi = (r1[f].subs(x(t),u).subs(y(t),v))/num
     if R1 and R2:
         return 'type4'
@@ -3558,8 +3560,8 @@ def ode_2nd_power_series_ordinary(eq, func, order, match):
     seriesdict = {}
     recurr = Function("r")
 
-    # Generating the recurrence relation which works this way
-    # a] For the second order term the summation begins at n = 2. The coefficients
+    # Generating the recurrence relation which works this way:
+    # for the second order term the summation begins at n = 2. The coefficients
     # p is multiplied with an*(n - 1)*(n - 2)*x**n-2 and a substitution is made such that
     # the exponent of x becomes n.
     # For example, if p is x, then the second degree recurrence term is
@@ -3635,8 +3637,7 @@ def ode_2nd_power_series_ordinary(eq, func, order, match):
     # Checking how many values are already present
     tcounter = len([t for t in finaldict.values() if t])
 
-    for count in range(tcounter, terms - 3):  # Assuming c0 and c1 to be arbitrary
-    #while tcounter < terms - 2:  # Assuming c0 and c1 to be arbitrary
+    for _ in range(tcounter, terms - 3):  # Assuming c0 and c1 to be arbitrary
         check = rhs.subs(n, startiter)
         nlhs = lhs.subs(n, startiter)
         nrhs = check.subs(finaldict)
@@ -4265,7 +4266,7 @@ def abc(eq):
         if eq is a*x + b*f(x) + c, else None.
         '''
         eq = _mexpand(eq)
-        c = eq.as_independent(x, f(x), as_Add = True)[0]
+        c = eq.as_independent(x, f(x), as_Add=True)[0]
         if not c.is_Rational:
             return
         a = eq.coeff(x)
@@ -6800,10 +6801,10 @@ def sysode_linear_2eq_order2(match_):
     r = dict()
     t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
     for i in range(2):
-        eqs = 0
+        eqs = []
         for terms in Add.make_args(eq[i]):
-            eqs += terms/fc[i,func[i],2]
-        eq[i] = eqs
+            eqs.append(terms/fc[i,func[i],2])
+        eq[i] = Add(*eqs)
     # for equations Eq(diff(x(t),t,t), a1*diff(x(t),t)+b1*diff(y(t),t)+c1*x(t)+d1*y(t)+e1)
     # and Eq(a2*diff(y(t),t,t), a2*diff(x(t),t)+b2*diff(y(t),t)+c2*x(t)+d2*y(t)+e2)
     r['a1'] = -fc[0,x(t),1]/fc[0,x(t),2] ; r['a2'] = -fc[1,x(t),1]/fc[1,y(t),2]
@@ -7184,14 +7185,16 @@ def _linear_2eq_order2_type6(x, y, t, r, eq):
     C1, C2, C3, C4 = get_numbered_constants(eq, num=4)
     k = Symbol('k')
     z = Function('z')
-    num, denum = cancel((r['c1']*x(t) + r['d1']*y(t))/(r['c2']*x(t) + r['d2']*y(t))).as_numer_denom()
+    num, den = cancel(
+        (r['c1']*x(t) + r['d1']*y(t))/
+        (r['c2']*x(t) + r['d2']*y(t))).as_numer_denom()
     f = r['c1']/num.coeff(x(t))
     a1 = num.coeff(x(t))
     b1 = num.coeff(y(t))
-    a2 = denum.coeff(x(t))
-    b2 = denum.coeff(y(t))
+    a2 = den.coeff(x(t))
+    b2 = den.coeff(y(t))
     chareq = k**2 - (a1 + b2)*k + a1*b2 - a2*b1
-    [k1, k2] = [rootof(chareq, k) for k in range(Poly(chareq).degree())]
+    k1, k2 = [rootof(chareq, k) for k in range(Poly(chareq).degree())]
     z1 = dsolve(diff(z(t),t,t) - k1*f*z(t)).rhs
     z2 = dsolve(diff(z(t),t,t) - k2*f*z(t)).rhs
     sol1 = (k1*z2 - k2*z1 + a1*(z1 - z2))/(a2*(k1-k2))
@@ -7229,12 +7232,14 @@ def _linear_2eq_order2_type7(x, y, t, r, eq):
     """
     C1, C2, C3, C4 = get_numbered_constants(eq, num=4)
     k = Symbol('k')
-    num, denum = cancel((r['a1']*x(t) + r['b1']*y(t))/(r['a2']*x(t) + r['b2']*y(t))).as_numer_denom()
+    num, den = cancel(
+        (r['a1']*x(t) + r['b1']*y(t))/
+        (r['a2']*x(t) + r['b2']*y(t))).as_numer_denom()
     f = r['a1']/num.coeff(x(t))
     a1 = num.coeff(x(t))
     b1 = num.coeff(y(t))
-    a2 = denum.coeff(x(t))
-    b2 = denum.coeff(y(t))
+    a2 = den.coeff(x(t))
+    b2 = den.coeff(y(t))
     chareq = k**2 - (a1 + b2)*k + a1*b2 - a2*b1
     [k1, k2] = [rootof(chareq, k) for k in range(Poly(chareq).degree())]
     F = Integral(f, t)
@@ -7283,10 +7288,10 @@ def _linear_2eq_order2_type8(x, y, t, r, eq):
 
     """
     C1, C2, C3, C4 = get_numbered_constants(eq, num=4)
-    num, denum = cancel(r['d1']/r['c2']).as_numer_denom()
+    num, den = cancel(r['d1']/r['c2']).as_numer_denom()
     f = -r['d1']/num
     a = num
-    b = denum
+    b = den
     mul = sqrt(abs(a*b))
     Igral = Integral(t*f, t)
     if a*b > 0:
@@ -7386,8 +7391,8 @@ def _linear_2eq_order2_type10(x, y, t, r, eq):
     q = Wild('q', exclude=[t, t**2])
     s = Wild('s', exclude=[t, t**2])
     n = Wild('n', exclude=[t, t**2])
-    num, denum = r['c1'].as_numer_denom()
-    dic = denum.match((n*(p*t**2+q*t+s)**2).expand())
+    num, den = r['c1'].as_numer_denom()
+    dic = den.match((n*(p*t**2+q*t+s)**2).expand())
     eqz = dic[p]*t**2 + dic[q]*t + dic[s]
     a = num/dic[n]
     b = cancel(r['d1']*eqz**2)
@@ -7938,9 +7943,11 @@ def _nonlinear_2eq_order1_type4(x, y, t, eq):
     g2 = Wild('g2', exclude=[u,t])
     r1 = eq[0].match(diff(x(t),t) - f)
     r2 = eq[1].match(diff(y(t),t) - g)
-    num, denum = ((r1[f].subs(x(t),u).subs(y(t),v))/(r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom()
+    num, den = (
+        (r1[f].subs(x(t),u).subs(y(t),v))/
+        (r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom()
     R1 = num.match(f1*g1)
-    R2 = denum.match(f2*g2)
+    R2 = den.match(f2*g2)
     phi = (r1[f].subs(x(t),u).subs(y(t),v))/num
     F1 = R1[f1]; F2 = R2[f2]
     G1 = R1[g1]; G2 = R2[g2]
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 999d57cbda..316911ef64 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -12,7 +12,7 @@
 from sympy.core.numbers import I, Number, Rational, oo
 from sympy.core.function import (Lambda, expand, expand_complex)
 from sympy.core.relational import Eq
-from sympy.simplify.simplify import simplify, fraction, trigsimp
+from sympy.simplify.simplify import simplify, fraction, trigsimp, signsimp
 from sympy.core.symbol import Symbol
 from sympy.functions import (log, Abs, tan, cot, sin, cos, sec, csc, exp,
                              acos, asin, atan, acsc, asec, arg,
@@ -25,95 +25,125 @@
 from sympy.matrices import Matrix
 from sympy.polys import (roots, Poly, degree, together, PolynomialError,
                          RootOf)
-from sympy.solvers.solvers import checksol, denoms
+from sympy.solvers.solvers import checksol, denoms, unrad
 from sympy.solvers.inequalities import solve_univariate_inequality
 from sympy.utilities import filldedent
 
 import warnings
 
 
-def invert_real(f_x, y, x):
-    """ Inverts a real valued function
-
-    Reduces the real valued equation ``f(x) = y`` to a set of equations ``{g(x)
-    = h_1(y), g(x) = h_2(y), ..., g(x) = h_n(y) }`` where ``g(x)`` is a simpler
-    function than ``f(x)``.  The return value is a tuple ``(g(x), set_h)``,
-    where ``g(x)`` is a function of ``x`` and ``set_h`` is the set of
-    functions ``{h_1(y), h_2(y), ..., h_n(y)}``.
+def _invert(f_x, y, x, domain=S.Complexes):
+    """
+    Reduce the complex valued equation ``f(x) = y`` to a set of equations
+    ``{g(x) = h_1(y), g(x) = h_2(y), ..., g(x) = h_n(y) }`` where ``g(x)`` is
+    a simpler function than ``f(x)``.  The return value is a tuple ``(g(x),
+    set_h)``, where ``g(x)`` is a function of ``x`` and ``set_h`` is
+    the set of function ``{h_1(y), h_2(y), ..., h_n(y)}``.
     Here, ``y`` is not necessarily a symbol.
 
-    The ``set_h`` contains the functions along with the information about their
-    domain in which they are valid, through set operations. For instance, if
-    ``y = Abs(x) - n``, is inverted, then, the ``set_h`` doesn't simply
-    return `{-n, n}`, as it doesn't explicitly mentions about the nature of
-    `n` rather it will return:
+    The ``set_h`` contains the functions along with the information
+    about their domain in which they are valid, through set
+    operations. For instance, if ``y = Abs(x) - n``, is inverted
+    in the real domain, then, the ``set_h`` doesn't simply return
+    `{-n, n}`, as the nature of `n` is unknown; rather it will return:
     `Intersection([0, oo) {n}) U Intersection((-oo, 0], {-n})`
 
+    By default, the complex domain is used but note that inverting even
+    seemingly simple functions like ``exp(x)`` can give very different
+    result in the complex domain than are obtained in the real domain.
+    (In the case of ``exp(x)``, the inversion via ``log`` is multi-valued
+    in the complex domain, having infinitely many branches.)
+
+    If you are working with real values only (or you are not sure which
+    function to use) you should probably use set the domain to
+    ``S.Reals`` (or use `invert\_real` which does that automatically).
+
 
     Examples
     ========
 
-    >>> from sympy.solvers.solveset import invert_real
-    >>> from sympy import tan, Abs, exp
-    >>> from sympy.abc import x, y, n
-    >>> invert_real(exp(x), 1, x)
-    (x, {0})
-    >>> invert_real(tan(x), y, x)
-    (x, ImageSet(Lambda(_n, _n*pi + atan(y)), Integers()))
+    >>> from sympy.solvers.solveset import invert_complex, invert_real
+    >>> from sympy.abc import x, y
+    >>> from sympy import exp, log
 
+    When does exp(x) == y?
 
-    * ``set_h`` containing information about the domain
+    >>> invert_complex(exp(x), y, x)
+    (x, ImageSet(Lambda(_n, I*(2*_n*pi + arg(y)) + log(Abs(y))), Integers()))
+    >>> invert_real(exp(x), y, x)
+    (x, Intersection((-oo, oo), {log(y)}))
 
-    >>> invert_real(Abs(x**31 + x), y, x)
-    (x**31 + x, Intersection([0, oo), {y}) U Intersection((-oo, 0], {-y}))
-    >>> invert_real(exp(Abs(x)), y, x)
-    (x, Intersection([0, oo), {log(y)}) U Intersection((-oo, 0], {-log(y)}))
+    When does exp(x) == 1?
+
+    >>> invert_complex(exp(x), 1, x)
+    (x, ImageSet(Lambda(_n, 2*_n*I*pi), Integers()))
+    >>> invert_real(exp(x), 1, x)
+    (x, {0})
 
     See Also
     ========
-    invert_complex
+    invert_real, invert_complex
     """
+    x = sympify(x)
+    if not x.is_Symbol:
+        raise ValueError("x must be a symbol")
+    f_x = sympify(f_x)
+    if not f_x.has(x):
+        raise ValueError("Inverse of constant function doesn't exist")
     y = sympify(y)
-    if not y.has(x):
-        return _invert_real(f_x, FiniteSet(y), x)
+    if y.has(x):
+        raise ValueError("y should be independent of x ")
+
+    if domain.is_subset(S.Reals):
+        x, s = _invert_real(f_x, FiniteSet(y), x)
     else:
-        raise ValueError(" y should be independent of x ")
+        x, s = _invert_complex(f_x, FiniteSet(y), x)
+    return x, s.intersection(domain) if isinstance(s, FiniteSet) else s
 
 
-def _invert_real(f, g_ys, symbol):
-    """ Helper function for invert_real """
+invert_complex = _invert
+
+
+def invert_real(f_x, y, x, domain=S.Reals):
+    return _invert(f_x, y, x, domain)
 
-    if not f.has(symbol):
-        raise ValueError("Inverse of constant function doesn't exist")
 
-    if f is symbol:
+def _invert_real(f, g_ys, symbol):
+    """Helper function for _invert."""
+
+    if f == symbol:
         return (f, g_ys)
 
-    n = Dummy('n')
-    if hasattr(f, 'inverse') and not isinstance(f, TrigonometricFunction) and \
-            not isinstance(f, HyperbolicFunction):
+    n = Dummy('n', real=True)
+
+    if hasattr(f, 'inverse') and not isinstance(f, (
+            TrigonometricFunction,
+            HyperbolicFunction,
+            )):
         if len(f.args) > 1:
             raise ValueError("Only functions with one argument are supported.")
         return _invert_real(f.args[0],
-                            imageset(Lambda(n, f.inverse()(n)), g_ys), symbol)
+                            imageset(Lambda(n, f.inverse()(n)), g_ys),
+                            symbol)
 
     if isinstance(f, Abs):
+        pos = Interval(0, S.Infinity)
+        neg = Interval(S.NegativeInfinity, 0)
         return _invert_real(f.args[0],
-                            Union(imageset(Lambda(n, n), g_ys).intersect(Interval(0, oo)),
-                                  imageset(Lambda(n, -n), g_ys).intersect(Interval(-oo, 0))),
-                            symbol)
+                    Union(imageset(Lambda(n, n), g_ys).intersect(pos),
+                          imageset(Lambda(n, -n), g_ys).intersect(neg)), symbol)
 
     if f.is_Add:
         # f = g + h
         g, h = f.as_independent(symbol)
-        if g != S.Zero:
+        if g is not S.Zero:
             return _invert_real(h, imageset(Lambda(n, n - g), g_ys), symbol)
 
     if f.is_Mul:
         # f = g*h
         g, h = f.as_independent(symbol)
 
-        if g != S.One:
+        if g is not S.One:
             return _invert_real(h, imageset(Lambda(n, n/g), g_ys), symbol)
 
     if f.is_Pow:
@@ -141,116 +171,51 @@ def _invert_real(f, g_ys, symbol):
                 return _invert_real(base, res, symbol)
 
         if not base_has_sym:
-            return _invert_real(expo, imageset(Lambda(n, log(n)/log(base)),
-                                               g_ys), symbol)
-
-    if isinstance(f, sin):
-        n = Dummy('n')
-        if isinstance(g_ys, FiniteSet):
-            sin_invs = Union(*[imageset(Lambda(n, n*pi + (-1)**n*asin(g_y)), \
-                                        S.Integers) for g_y in g_ys])
-            return _invert_real(f.args[0], sin_invs, symbol)
+            return _invert_real(expo,
+                imageset(Lambda(n, log(n)/log(base)), g_ys), symbol)
 
-    if isinstance(f, csc):
-        n = Dummy('n')
+    if isinstance(f, TrigonometricFunction):
         if isinstance(g_ys, FiniteSet):
-            csc_invs = Union(*[imageset(Lambda(n, n*pi + (-1)**n*acsc(g_y)), \
-                                        S.Integers) for g_y in g_ys])
-            return _invert_real(f.args[0], csc_invs, symbol)
-
-    if isinstance(f, cos):
-        n = Dummy('n')
-        if isinstance(g_ys, FiniteSet):
-            cos_invs_f1 = Union(*[imageset(Lambda(n, 2*n*pi + acos(g_y)), \
-                                        S.Integers) for g_y in g_ys])
-            cos_invs_f2 = Union(*[imageset(Lambda(n, 2*n*pi - acos(g_y)), \
-                                        S.Integers) for g_y in g_ys])
-            cos_invs = Union(cos_invs_f1, cos_invs_f2)
-            return _invert_real(f.args[0], cos_invs, symbol)
-
-    if isinstance(f, sec):
-        n = Dummy('n')
-        if isinstance(g_ys, FiniteSet):
-            sec_invs_f1 = Union(*[imageset(Lambda(n, 2*n*pi + asec(g_y)), \
-                                        S.Integers) for g_y in g_ys])
-            sec_invs_f2 = Union(*[imageset(Lambda(n, 2*n*pi - asec(g_y)), \
-                                        S.Integers) for g_y in g_ys])
-            sec_invs = Union(sec_invs_f1, sec_invs_f2)
-            return _invert_real(f.args[0], sec_invs, symbol)
-
-    if isinstance(f, tan) or isinstance(f, cot):
-        n = Dummy('n')
-        if isinstance(g_ys, FiniteSet):
-            tan_cot_invs = Union(*[imageset(Lambda(n, n*pi + f.inverse()(g_y)), \
-                                        S.Integers) for g_y in g_ys])
-            return _invert_real(f.args[0], tan_cot_invs, symbol)
+            def inv(trig):
+                if isinstance(f, (sin, csc)):
+                    F = asin if isinstance(f, sin) else acsc
+                    return (lambda a: n*pi + (-1)**n*F(a),)
+                if isinstance(f, (cos, sec)):
+                    F = acos if isinstance(f, cos) else asec
+                    return (
+                        lambda a: 2*n*pi + F(a),
+                        lambda a: 2*n*pi - F(a),)
+                if isinstance(f, (tan, cot)):
+                    return (lambda a: n*pi + f.inverse()(a),)
+
+            n = Dummy('n', integer=True)
+            invs = S.EmptySet
+            for L in inv(f):
+                invs += Union(*[imageset(Lambda(n, L(g)), S.Integers) for g in g_ys])
+            return _invert_real(f.args[0], invs, symbol)
 
     return (f, g_ys)
 
 
-def invert_complex(f_x, y, x):
-    """ Inverts a complex valued function.
-
-    Reduces the complex valued equation ``f(x) = y`` to a set of equations
-    ``{g(x) = h_1(y), g(x) = h_2(y), ..., g(x) = h_n(y) }`` where ``g(x)`` is
-    a simpler function than ``f(x)``.  The return value is a tuple ``(g(x),
-    set_h)``, where ``g(x)`` is a function of ``x`` and ``set_h`` is
-    the set of function ``{h_1(y), h_2(y), ..., h_n(y)}``.
-    Here, ``y`` is not necessarily a symbol.
-
-    Note that `invert\_complex` and `invert\_real` don't always produce the
-    same result even for a seemingly simple function like ``exp(x)`` because
-    the complex extension of real valued ``log`` is multivariate in the complex
-    system and has infinitely many branches. If you are working with real
-    values only or you are not sure with function to use you should use
-    `invert\_real`.
-
-
-    Examples
-    ========
-
-    >>> from sympy.solvers.solveset import invert_complex
-    >>> from sympy.abc import x, y
-    >>> from sympy import exp, log
-    >>> invert_complex(log(x), y, x)
-    (x, {exp(y)})
-    >>> invert_complex(log(x), 0, x)  # Second parameter is not a symbol
-    (x, {1})
-    >>> invert_complex(exp(x), y, x)
-    (x, ImageSet(Lambda(_n, I*(2*_n*pi + arg(y)) + log(Abs(y))), Integers()))
-
-    See Also
-    ========
-    invert_real
-    """
-    y = sympify(y)
-    if not y.has(x):
-        return _invert_complex(f_x, FiniteSet(y), x)
-    else:
-        raise ValueError(" y should be independent of x ")
-
-
 def _invert_complex(f, g_ys, symbol):
-    """ Helper function for invert_complex """
+    """Helper function for _invert."""
 
-    if not f.has(symbol):
-        raise ValueError("Inverse of constant function doesn't exist")
-
-    if f is symbol:
+    if f == symbol:
         return (f, g_ys)
 
     n = Dummy('n')
+
     if f.is_Add:
         # f = g + h
         g, h = f.as_independent(symbol)
-        if g != S.Zero:
+        if g is not S.Zero:
             return _invert_complex(h, imageset(Lambda(n, n - g), g_ys), symbol)
 
     if f.is_Mul:
         # f = g*h
         g, h = f.as_independent(symbol)
 
-        if g != S.One:
+        if g is not S.One:
             return _invert_complex(h, imageset(Lambda(n, n/g), g_ys), symbol)
 
     if hasattr(f, 'inverse') and \
@@ -267,6 +232,7 @@ def _invert_complex(f, g_ys, symbol):
                                                log(Abs(g_y))), S.Integers)
                                for g_y in g_ys if g_y != 0])
             return _invert_complex(f.args[0], exp_invs, symbol)
+
     return (f, g_ys)
 
 
@@ -317,13 +283,23 @@ def _domain_check(f, symbol, p):
                     for g in f.args])
 
 
-def _is_finite_with_finite_vars(f):
+def _is_finite_with_finite_vars(f, domain=S.Complexes):
     """
-    Return True if the given expression is finite when all free symbols
-    (that are not already specified as finite) are made finite.
+    Return True if the given expression is finite. For symbols that
+    don't assign a value for `complex` and/or `real`, the domain will
+    be used to assign a value; symbols that don't assign a value
+    for `finite` will be made finite. All other assumptions are
+    left unmodified.
     """
-    reps = dict([(s, Dummy(s.name, finite=True, **s.assumptions0))
-                for s in f.free_symbols if s.is_finite is None])
+    def assumptions(s):
+        A = s.assumptions0
+        if A.get('finite', None) is None:
+            A['finite'] = True
+        A.setdefault('complex', True)
+        A.setdefault('real', domain.is_subset(S.Reals))
+        return A
+
+    reps = dict([(s, Dummy(**assumptions(s))) for s in f.free_symbols])
     return f.xreplace(reps).is_finite
 
 
@@ -377,161 +353,15 @@ def _is_function_class_equation(func_class, f, symbol):
         return False
 
 
-def solveset_real(f, symbol):
-    """ Solves a real valued equation.
-
-    Parameters
-    ==========
-
-    f : Expr
-        The target equation
-    symbol : Symbol
-        The variable for which the equation is solved
-
-    Returns
-    =======
-
-    Set
-        A set of values for `symbol` for which `f` is equal to
-        zero. An `EmptySet` is returned if no solution is found.
-        A `ConditionSet` is returned as unsolved object if algorithms
-        to evaluate complete solutions are not yet implemented.
-
-    `solveset_real` claims to be complete in the set of the solution it
-    returns.
-
-    Raises
-    ======
-
-    NotImplementedError
-        Algorithms to solve inequalities in complex domain are
-        not yet implemented.
-    ValueError
-        The input is not valid.
-    RuntimeError
-        It is a bug, please report to the github issue tracker.
-
-
-    See Also
-    =======
-
-    solveset_complex : solver for complex domain
-
-    Examples
-    ========
-
-    >>> from sympy import Symbol, exp, sin, sqrt, I
-    >>> from sympy.solvers.solveset import solveset_real
-    >>> x = Symbol('x', real=True)
-    >>> a = Symbol('a', real=True, finite=True, positive=True)
-    >>> solveset_real(x**2 - 1, x)
-    {-1, 1}
-    >>> solveset_real(sqrt(5*x + 6) - 2 - x, x)
-    {-1, 2}
-    >>> solveset_real(x - I, x)
-    EmptySet()
-    >>> solveset_real(x - a, x)
-    {a}
-    >>> solveset_real(exp(x) - a, x)
-    {log(a)}
-
-    * In case the equation has infinitely many solutions an infinitely indexed
-      `ImageSet` is returned.
-
-    >>> solveset_real(sin(x) - 1, x)
-    ImageSet(Lambda(_n, 2*_n*pi + pi/2), Integers())
-
-    * If the equation is true for any arbitrary value of the symbol a `S.Reals`
-      set is returned.
-
-    >>> solveset_real(x - x, x)
-    (-oo, oo)
-
-    """
-    if not getattr(symbol, 'is_Symbol', False):
-        raise ValueError('A Symbol must be given, not type %s: %s' %
-            (type(symbol), symbol))
-
-    f = sympify(f)
-    if not isinstance(f, (Expr, Number)):
-        raise ValueError("%s is not a valid SymPy expression" % (f))
-
-    original_eq = f
-    f = together(f)
-
-    # In this, unlike in solveset_complex, expression should only
-    # be expanded when fraction(f)[1] does not contain the symbol
-    # for which we are solving
-    if not symbol in fraction(f)[1].free_symbols and f.is_rational_function():
-        f = expand(f)
-
-    f = piecewise_fold(f)
-
-    result = EmptySet()
-
-    if f.expand().is_zero:
-        return S.Reals
-    elif not f.has(symbol):
-        return EmptySet()
-    elif f.is_Mul and all([_is_finite_with_finite_vars(m) for m in f.args]):
-        # if f(x) and g(x) are both finite we can say that the solution of
-        # f(x)*g(x) == 0 is same as Union(f(x) == 0, g(x) == 0) is not true in
-        # general. g(x) can grow to infinitely large for the values where
-        # f(x) == 0. To be sure that we are not silently allowing any
-        # wrong solutions we are using this technique only if both f and g are
-        # finite for a finite input.
-        result = Union(*[solveset_real(m, symbol) for m in f.args])
-    elif _is_function_class_equation(TrigonometricFunction, f, symbol) or \
-            _is_function_class_equation(HyperbolicFunction, f, symbol):
-        result = _solve_real_trig(f, symbol)
-    elif f.is_Piecewise:
-        result = EmptySet()
-        expr_set_pairs = f.as_expr_set_pairs()
-        for (expr, in_set) in expr_set_pairs:
-            solns = solveset_real(expr, symbol).intersect(in_set)
-            result = result + solns
-    else:
-        lhs, rhs_s = invert_real(f, 0, symbol)
-        if lhs == symbol:
-            result = rhs_s
-        elif isinstance(rhs_s, FiniteSet):
-            equations = [lhs - rhs for rhs in rhs_s]
-            for equation in equations:
-                if equation == f:
-                    if any(_has_rational_power(g, symbol)[0]
-                           for g in equation.args):
-                        result += _solve_radical(equation,
-                                                 symbol,
-                                                 solveset_real)
-                    elif equation.has(Abs):
-                        result += _solve_abs(f, symbol)
-                    else:
-                        result += _solve_as_rational(equation, symbol,
-                                                     solveset_solver=solveset_real,
-                                                     as_poly_solver=_solve_as_poly_real)
-                else:
-                    result += solveset_real(equation, symbol)
-        else:
-            result = ConditionSet(symbol, Eq(f, 0), S.Reals)
-
-    if isinstance(result, FiniteSet):
-        result = [s for s in result
-                  if isinstance(s, RootOf)
-                  or domain_check(original_eq, symbol, s)]
-        return FiniteSet(*result).intersect(S.Reals)
-    else:
-        return result.intersect(S.Reals)
-
-
-def _solve_as_rational(f, symbol, solveset_solver, as_poly_solver):
+def _solve_as_rational(f, symbol, domain):
     """ solve rational functions"""
     f = together(f, deep=True)
     g, h = fraction(f)
     if not h.has(symbol):
-        return as_poly_solver(g, symbol)
+        return _solve_as_poly(g, symbol, domain)
     else:
-        valid_solns = solveset_solver(g, symbol)
-        invalid_solns = solveset_solver(h, symbol)
+        valid_solns = _solveset(g, symbol, domain)
+        invalid_solns = _solveset(h, symbol, domain)
         return valid_solns - invalid_solns
 
 
@@ -559,14 +389,13 @@ def _solve_real_trig(f, symbol):
         return ConditionSet(symbol, Eq(f_original, 0), S.Reals)
 
 
-def _solve_as_poly(f, symbol, solveset_solver, invert_func):
+def _solve_as_poly(f, symbol, domain=S.Complexes):
     """
     Solve the equation using polynomial techniques if it already is a
     polynomial equation or, with a change of variables, can be made so.
     """
     result = None
     if f.is_polynomial(symbol):
-
         solns = roots(f, symbol, cubics=True, quartics=True,
                       quintics=True, domain='EX')
         num_roots = sum(solns.values())
@@ -578,11 +407,11 @@ def _solve_as_poly(f, symbol, solveset_solver, invert_func):
             if poly.degree() <= len(solns):
                 result = FiniteSet(*solns)
             else:
-                result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
+                result = ConditionSet(symbol, Eq(f, 0), domain)
     else:
         poly = Poly(f)
         if poly is None:
-            result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
+            result = ConditionSet(symbol, Eq(f, 0), domain)
         gens = [g for g in poly.gens if g.has(symbol)]
 
         if len(gens) == 1:
@@ -594,17 +423,18 @@ def _solve_as_poly(f, symbol, solveset_solver, invert_func):
                                           quintics=True).keys())
 
             if len(poly_solns) < deg:
-                result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
+                result = ConditionSet(symbol, Eq(f, 0), domain)
 
             if gen != symbol:
                 y = Dummy('y')
-                lhs, rhs_s = invert_func(gen, y, symbol)
-                if lhs is symbol:
+                inverter = invert_real if domain.is_subset(S.Reals) else invert_complex
+                lhs, rhs_s = inverter(gen, y, symbol)
+                if lhs == symbol:
                     result = Union(*[rhs_s.subs(y, s) for s in poly_solns])
                 else:
-                    result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
+                    result = ConditionSet(symbol, Eq(f, 0), domain)
         else:
-            result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
+            result = ConditionSet(symbol, Eq(f, 0), domain)
 
     if result is not None:
         if isinstance(result, FiniteSet):
@@ -616,29 +446,11 @@ def _solve_as_poly(f, symbol, solveset_solver, invert_func):
                     for s in result]):
                 s = Dummy('s')
                 result = imageset(Lambda(s, expand_complex(s)), result)
+        if isinstance(result, FiniteSet):
+            result = result.intersection(domain)
         return result
     else:
-        return ConditionSet(symbol, Eq(f, 0), S.Complexes)
-
-
-def _solve_as_poly_real(f, symbol):
-    """
-    Solve real valued equation with methods to solve polynomial
-    equations.
-    """
-    return _solve_as_poly(f, symbol,
-                          solveset_solver=solveset_real,
-                          invert_func=invert_real)
-
-
-def _solve_as_poly_complex(f, symbol):
-    """
-    Solve complex valued equation with methods to solve polynomial
-    equations.
-    """
-    return _solve_as_poly(f, symbol,
-                          solveset_solver=solveset_complex,
-                          invert_func=invert_complex)
+        return ConditionSet(symbol, Eq(f, 0), domain)
 
 
 def _has_rational_power(expr, symbol):
@@ -678,7 +490,6 @@ def _has_rational_power(expr, symbol):
 
 def _solve_radical(f, symbol, solveset_solver):
     """ Helper function to solve equations with radicals """
-    from sympy.solvers.solvers import unrad
     eq, cov = unrad(f)
     if not cov:
         result = solveset_solver(eq, symbol) - \
@@ -698,8 +509,12 @@ def _solve_radical(f, symbol, solveset_solver):
     return FiniteSet(*[s for s in result if checksol(f, symbol, s) is True])
 
 
-def _solve_abs(f, symbol):
+def _solve_abs(f, symbol, domain):
     """ Helper function to solve equation involving absolute value function """
+    if not domain.is_subset(S.Reals):
+        raise ValueError(filldedent('''
+            Absolute values cannot be inverted in the
+            complex domain.'''))
     p, q, r = Wild('p'), Wild('q'), Wild('r')
     pattern_match = f.match(p*Abs(q) + r) or {}
     if not pattern_match.get(p, S.Zero).is_zero:
@@ -715,117 +530,114 @@ def _solve_abs(f, symbol):
                                            symbol).intersect(q_neg_cond)
         return Union(sols_q_pos, sols_q_neg)
     else:
-        return ConditionSet(symbol, Eq(f, 0), S.Complexes)
-
-
-def solveset_complex(f, symbol):
-    """ Solve a complex valued equation.
-
-    Parameters
-    ==========
-
-    f : Expr
-        The target equation
-    symbol : Symbol
-        The variable for which the equation is solved
-
-    Returns
-    =======
-
-    Set
-        A set of values for `symbol` for which `f` equal to
-        zero. An `EmptySet` is returned if no solution is found.
-        A `ConditionSet` is returned as an unsolved object if algorithms
-        to evaluate complete solutions are not yet implemented.
-
-    `solveset_complex` claims to be complete in the solution set that
-    it returns.
-
-    Raises
-    ======
+        return ConditionSet(symbol, Eq(f, 0), domain)
 
-    NotImplementedError
-        The algorithms to solve inequalities in complex domain  are
-        not yet implemented.
-    ValueError
-        The input is not valid.
-    RuntimeError
-        It is a bug, please report to the github issue tracker.
 
-    See Also
-    ========
 
-    solveset_real: solver for real domain
+def _solveset(f, symbol, domain, _check=False):
+    """Helper for solveset to return a result from an expression
+    that has already been sympify'ed and is known to contain the
+    given symbol."""
+    # _check controls whether the answer is checked or not
 
-    Examples
-    ========
-
-    >>> from sympy import Symbol, exp
-    >>> from sympy.solvers.solveset import solveset_complex
-    >>> from sympy.abc import x, a, b, c
-    >>> solveset_complex(a*x**2 + b*x +c, x)
-    {-b/(2*a) - sqrt(-4*a*c + b**2)/(2*a), -b/(2*a) + sqrt(-4*a*c + b**2)/(2*a)}
-
-    * Due to the fact that complex extension of my real valued functions are
-      multivariate even some simple equations can have infinitely many
-      solution.
-
-    >>> solveset_complex(exp(x) - 1, x)
-    ImageSet(Lambda(_n, 2*_n*I*pi), Integers())
-
-    """
-    if not getattr(symbol, 'is_Symbol', False):
-        raise ValueError('A Symbol must be given, not type %s: %s' %
-            (type(symbol), symbol))
+    orig_f = f
+    f = together(f)
+    if f.is_Mul:
+        _, f = f.as_independent(symbol, as_Add=False)
+    if f.is_Add:
+        a, h = f.as_independent(symbol)
+        m, h = h.as_independent(symbol, as_Add=False)
+        f = a/m + h  # XXX condition `m != 0` should be added to soln
+    f = piecewise_fold(f)
 
-    f = sympify(f)
-    original_eq = f
-    if not isinstance(f, (Expr, Number)):
-        raise ValueError(" %s is not a valid sympy expression" % (f))
+    # assign the solvers to use
+    solver = lambda f, x, domain=domain: _solveset(f, x, domain)
+    if domain.is_subset(S.Reals):
+        inverter_func = invert_real
+    else:
+        inverter_func = invert_complex
+    inverter = lambda f, rhs, symbol: inverter_func(f, rhs, symbol, domain)
 
-    f = together(f)
-    # Without this equations like a + 4*x**2 - E keep oscillating
-    # into form  a/4 + x**2 - E/4 and (a + 4*x**2 - E)/4
-    if not fraction(f)[1].has(symbol):
-        f = expand(f)
+    result = EmptySet()
 
-    if f.is_zero:
-        return S.Complexes
+    if f.expand().is_zero:
+        return domain
     elif not f.has(symbol):
+        return EmptySet()
+    elif f.is_Mul and all(_is_finite_with_finite_vars(m, domain)
+            for m in f.args):
+        # if f(x) and g(x) are both finite we can say that the solution of
+        # f(x)*g(x) == 0 is same as Union(f(x) == 0, g(x) == 0) is not true in
+        # general. g(x) can grow to infinitely large for the values where
+        # f(x) == 0. To be sure that we are not silently allowing any
+        # wrong solutions we are using this technique only if both f and g are
+        # finite for a finite input.
+        result = Union(*[solver(m, symbol) for m in f.args])
+    elif _is_function_class_equation(TrigonometricFunction, f, symbol) or \
+            _is_function_class_equation(HyperbolicFunction, f, symbol):
+        result = _solve_real_trig(f, symbol)
+    elif f.is_Piecewise:
+        dom = domain
         result = EmptySet()
-    elif f.is_Mul and all([_is_finite_with_finite_vars(m) for m in f.args]):
-        result = Union(*[solveset_complex(m, symbol) for m in f.args])
+        expr_set_pairs = f.as_expr_set_pairs()
+        for (expr, in_set) in expr_set_pairs:
+            if in_set.is_Relational:
+                in_set = in_set.as_set()
+            if in_set.is_Interval:
+                dom -= in_set
+            solns = solver(expr, symbol, in_set)
+            result += solns
     else:
-        lhs, rhs_s = invert_complex(f, 0, symbol)
+        lhs, rhs_s = inverter(f, 0, symbol)
         if lhs == symbol:
+            # do some very minimal simplification since
+            # repeated inversion may have left the result
+            # in a state that other solvers (e.g. poly)
+            # would have simplified; this is done here
+            # rather than in the inverter since here it
+            # is only done once whereas there it would
+            # be repeated for each step of the inversion
+            if isinstance(rhs_s, FiniteSet):
+                rhs_s = FiniteSet(*[Mul(*
+                    signsimp(i).as_content_primitive())
+                    for i in rhs_s])
             result = rhs_s
         elif isinstance(rhs_s, FiniteSet):
-            equations = [lhs - rhs for rhs in rhs_s]
-            result = EmptySet()
-            for equation in equations:
+            for equation in [lhs - rhs for rhs in rhs_s]:
                 if equation == f:
                     if any(_has_rational_power(g, symbol)[0]
                            for g in equation.args) or _has_rational_power(
                            equation, symbol)[0]:
                         result += _solve_radical(equation,
                                                  symbol,
-                                                 solveset_complex)
+                                                 solver)
+                    elif equation.has(Abs):
+                        result += _solve_abs(f, symbol, domain)
                     else:
-                        result += _solve_as_rational(equation, symbol,
-                                                 solveset_solver=solveset_complex,
-                                                 as_poly_solver=_solve_as_poly_complex)
+                        result += _solve_as_rational(equation, symbol, domain)
                 else:
-                    result += solveset_complex(equation, symbol)
+                    result += solver(equation, symbol)
         else:
-            result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
+            result = ConditionSet(symbol, Eq(f, 0), domain)
 
-    if isinstance(result, FiniteSet):
-        result = [s for s in result
-                  if isinstance(s, RootOf)
-                  or domain_check(original_eq, symbol, s)]
-        return FiniteSet(*result)
-    else:
-        return result
+    if _check:
+        if isinstance(result, ConditionSet):
+            # it wasn't solved or has enumerated all conditions
+            # -- leave it alone
+            return result
+
+        # whittle away all but the symbol-containing core
+        # to use this for testing
+        fx = orig_f.as_independent(symbol, as_Add=True)[1]
+        fx = fx.as_independent(symbol, as_Add=False)[1]
+
+        if isinstance(result, FiniteSet):
+            # check the result for invalid solutions
+            result = FiniteSet(*[s for s in result
+                      if isinstance(s, RootOf)
+                      or domain_check(fx, symbol, s)])
+
+    return result
 
 
 def solveset(f, symbol=None, domain=S.Complexes):
@@ -864,12 +676,6 @@ def solveset(f, symbol=None, domain=S.Complexes):
         It is a bug, please report to the github issue tracker.
 
 
-    `solveset` uses two underlying functions `solveset_real` and
-    `solveset_complex` to solve equations. They are the solvers for real and
-    complex domain respectively. `solveset` ignores the assumptions on the
-    variable being solved for and instead, uses the `domain` parameter to
-    decide which solver to use.
-
     Notes
     =====
 
@@ -889,33 +695,49 @@ def solveset(f, symbol=None, domain=S.Complexes):
     Examples
     ========
 
-    >>> from sympy import exp, Symbol, Eq, pprint, S, solveset
-    >>> from sympy.abc import x
+    >>> from sympy import exp, sin, Symbol, pprint, S
+    >>> from sympy.solvers.solveset import solveset, solveset_real
+
+    * The default domain is complex. Not specifying a domain will lead
+      to the solving of the equation in the complex domain (and this
+      is not affected by the assumptions on the symbol):
 
-    * The default domain is complex. Not specifying a domain will lead to the
-      solving of the equation in the complex domain.
+    >>> x = Symbol('x')
+    >>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
+    {2*n*I*pi | n in Integers()}
 
+    >>> x = Symbol('x', real=True)
     >>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
     {2*n*I*pi | n in Integers()}
 
-    * If you want to solve equation in real domain by the `solveset`
-      interface, then specify that the domain is real. Alternatively use
-      `solveset\_real`.
+    * If you want to use `solveset` to solve the equation in the
+      real domain, provide a real domain. (Using `solveset\_real`
+      does this automatically.)
 
+    >>> R = S.Reals
     >>> x = Symbol('x')
-    >>> solveset(exp(x) - 1, x, S.Reals)
+    >>> solveset(exp(x) - 1, x, R)
     {0}
-    >>> solveset(Eq(exp(x), 1), x, S.Reals)
+    >>> solveset_real(exp(x) - 1, x)
     {0}
 
+    The solution is mostly unaffected by assumptions on the symbol,
+    but there may be some slight difference:
+
+    >>> pprint(solveset(sin(x)/x,x), use_unicode=False)
+    ({2*n*pi | n in Integers()} \ {0}) U ({2*n*pi + pi | n in Integers()} \ {0})
+
+    >>> p = Symbol('p', positive=True)
+    >>> pprint(solveset(sin(p)/p, p), use_unicode=False)
+    {2*n*pi | n in Integers()} U {2*n*pi + pi | n in Integers()}
+
     * Inequalities can be solved over the real domain only. Use of a complex
       domain leads to a NotImplementedError.
 
-    >>> solveset(exp(x) > 1, x, S.Reals)
+    >>> solveset(exp(x) > 1, x, R)
     (0, oo)
 
     """
-
     f = sympify(f)
 
     if f is S.true:
@@ -924,6 +746,9 @@ def solveset(f, symbol=None, domain=S.Complexes):
     if f is S.false:
         return S.EmptySet
 
+    if not isinstance(f, (Expr, Number)):
+        raise ValueError("%s is not a valid SymPy expression" % (f))
+
     free_symbols = f.free_symbols
 
     if not free_symbols:
@@ -950,8 +775,7 @@ def solveset(f, symbol=None, domain=S.Complexes):
     if isinstance(f, Eq):
         from sympy.core import Add
         f = Add(f.lhs, - f.rhs, evaluate=False)
-
-    if f.is_Relational:
+    elif f.is_Relational:
         if not domain.is_subset(S.Reals):
             raise NotImplementedError(filldedent('''
                 Inequalities in the complex domain are
@@ -959,20 +783,28 @@ def solveset(f, symbol=None, domain=S.Complexes):
                 setting domain=S.Reals'''))
         try:
             result = solve_univariate_inequality(
-            f, symbol, relational=False).intersection(domain)
+            f, symbol, relational=False) - _invalid_solutions(
+            f, symbol, domain)
         except NotImplementedError:
             result = ConditionSet(symbol, f, domain)
         return result
 
-    if isinstance(f, (Expr, Number)):
-        if domain is S.Reals:
-            return solveset_real(f, symbol)
-        elif domain is S.Complexes:
-            return solveset_complex(f, symbol)
-        elif domain.is_subset(S.Reals):
-            return Intersection(solveset_real(f, symbol), domain)
-        else:
-            return Intersection(solveset_complex(f, symbol), domain)
+    return _solveset(f, symbol, domain, _check=True)
+
+
+def _invalid_solutions(f, symbol, domain):
+    bad = S.EmptySet
+    for d in denoms(f):
+        bad += _solveset(d, symbol, domain, _check=False)
+    return bad
+
+
+def solveset_real(f, symbol):
+    return solveset(f, symbol, S.Reals)
+
+
+def solveset_complex(f, symbol):
+    return solveset(f, symbol, S.Complexes)
 
 
 ###############################################################################


Test patch:
diff --git a/sympy/polys/tests/test_subresultants_qq_zz.py b/sympy/polys/tests/test_subresultants_qq_zz.py
index a377a936c2..737bdf7752 100644
--- a/sympy/polys/tests/test_subresultants_qq_zz.py
+++ b/sympy/polys/tests/test_subresultants_qq_zz.py
@@ -1,8 +1,6 @@
 from sympy import var, sturm, subresultants, prem, pquo
-from sympy.matrices import Matrix, eye
-from sympy.polys.subresultants_qq_zz import (sylvester, bezout,
-    subresultants_bezout, modified_subresultants_bezout,
-    process_bezout_output, backward_eye,
+from sympy.matrices import Matrix
+from sympy.polys.subresultants_qq_zz import (sylvester,
     sturm_pg, sturm_q, sturm_amv, euclid_pg, euclid_q,
     euclid_amv, modified_subresultants_pg, subresultants_pg,
     subresultants_amv_q, quo_z, rem_z, subresultants_amv,
@@ -34,44 +32,6 @@ def test_sylvester():
     assert sylvester(x**3 - 7*x + 7, 3*x**2 - 7, x, 2) == Matrix([
 [1, 0, -7,  7,  0,  0], [0, 3,  0, -7,  0,  0], [0, 1,  0, -7,  7,  0], [0, 0,  3,  0, -7,  0], [0, 0,  1,  0, -7,  7], [0, 0,  0,  3,  0, -7]])
 
-def test_bezout():
-    x = var('x')
-
-    p = -2*x**5+7*x**3+9*x**2-3*x+1
-    q = -10*x**4+21*x**2+18*x-3
-    assert bezout(p, q, x, 'bz').det() == sylvester(p, q, x, 2).det()
-    assert bezout(p, q, x, 'bz').det() != sylvester(p, q, x, 1).det()
-    assert bezout(p, q, x, 'prs') == backward_eye(5) * bezout(p, q, x, 'bz') * backward_eye(5)
-
-def test_subresultants_bezout():
-    x = var('x')
-
-    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
-    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
-    assert subresultants_bezout(p, q, x) == subresultants(p, q, x)
-    assert subresultants_bezout(p, q, x)[-1] == sylvester(p, q, x).det()
-    assert subresultants_bezout(p, q, x) != euclid_amv(p, q, x)
-    amv_factors = [1, 1, -1, 1, -1, 1]
-    assert subresultants_bezout(p, q, x) == [i*j for i, j in zip(amv_factors, modified_subresultants_amv(p, q, x))]
-
-    p = x**3 - 7*x + 7
-    q = 3*x**2 - 7
-    assert subresultants_bezout(p, q, x) == euclid_amv(p, q, x)
-
-def test_modified_subresultants_bezout():
-    x = var('x')
-
-    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
-    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
-    amv_factors = [1, 1, -1, 1, -1, 1]
-    assert modified_subresultants_bezout(p, q, x) == [i*j for i, j in zip(amv_factors, subresultants_amv(p, q, x))]
-    assert modified_subresultants_bezout(p, q, x)[-1] != sylvester(p + x**8, q, x).det()
-    assert modified_subresultants_bezout(p, q, x) != sturm_amv(p, q, x)
-
-    p = x**3 - 7*x + 7
-    q = 3*x**2 - 7
-    assert modified_subresultants_bezout(p, q, x) == sturm_amv(p, q, x)
-    assert modified_subresultants_bezout(-p, q, x) != sturm_amv(-p, q, x)
 
 def test_sturm_pg():
     x = var('x')
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index e8b15d0532..ee49f9cc72 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1,5 +1,5 @@
 from sympy import (
-    Abs, Dummy, Eq, Gt,
+    Abs, Dummy, Eq, Gt, Function,
     LambertW, Piecewise, Poly, Rational, S, Symbol, Matrix,
     asin, acos, acsc, asec, atan, atanh, cos, csc, erf, erfinv, erfc, erfcinv,
     exp, log, pi, sin, sinh, sec, sqrt, symbols,
@@ -15,7 +15,7 @@
 
 from sympy.polys.rootoftools import CRootOf
 
-from sympy.sets import (FiniteSet, ConditionSet)
+from sympy.sets import (FiniteSet, ConditionSet, Complement, ImageSet)
 
 from sympy.utilities.pytest import XFAIL, raises, skip, slow
 from sympy.utilities.randtest import verify_numerically as tn
@@ -39,13 +39,19 @@
 
 
 def test_invert_real():
-    x = Dummy(real=True)
+    x = Symbol('x', real=True)
+    y = Symbol('y')
     n = Symbol('n')
+    def ireal(x, s=S.Reals):
+        return Intersection(s, x)
 
     minus_n = Intersection(Interval(-oo, 0), FiniteSet(-n))
     plus_n = Intersection(Interval(0, oo), FiniteSet(n))
     assert solveset(abs(x) - n, x, S.Reals) == Union(minus_n, plus_n)
 
+    assert invert_real(exp(x), y, x) == (x, ireal(FiniteSet(log(y))))
+
+    y = Symbol('y', positive=True)
     n = Symbol('n', real=True)
     assert invert_real(x + 3, y, x) == (x, FiniteSet(y - 3))
     assert invert_real(x*3, y, x) == (x, FiniteSet(y / 3))
@@ -54,7 +60,7 @@ def test_invert_real():
     assert invert_real(exp(3*x), y, x) == (x, FiniteSet(log(y) / 3))
     assert invert_real(exp(x + 3), y, x) == (x, FiniteSet(log(y) - 3))
 
-    assert invert_real(exp(x) + 3, y, x) == (x, FiniteSet(log(y - 3)))
+    assert invert_real(exp(x) + 3, y, x) == (x, ireal(FiniteSet(log(y - 3))))
     assert invert_real(exp(x)*3, y, x) == (x, FiniteSet(log(y / 3)))
 
     assert invert_real(log(x), y, x) == (x, FiniteSet(exp(y)))
@@ -66,7 +72,7 @@ def test_invert_real():
     assert invert_real(Abs(x), y, x) == (x, Union(minus_y, plus_y))
 
     assert invert_real(2**x, y, x) == (x, FiniteSet(log(y)/log(2)))
-    assert invert_real(2**exp(x), y, x) == (x, FiniteSet(log(log(y)/log(2))))
+    assert invert_real(2**exp(x), y, x) == (x, ireal(FiniteSet(log(log(y)/log(2)))))
 
     assert invert_real(x**2, y, x) == (x, FiniteSet(sqrt(y), -sqrt(y)))
     assert invert_real(x**Rational(1, 2), y, x) == (x, FiniteSet(y**2))
@@ -157,8 +163,9 @@ def test_invert_complex():
 
     assert invert_complex(log(x), y, x) == (x, FiniteSet(exp(y)))
 
-    raises(ValueError, lambda: invert_real(S.One, y, x))
+    raises(ValueError, lambda: invert_real(1, y, x))
     raises(ValueError, lambda: invert_complex(x, x, x))
+    raises(ValueError, lambda: invert_complex(x, x, 1))
 
 
 def test_domain_check():
@@ -248,9 +255,7 @@ def test_solve_invert():
     assert solveset_real(3**(x + 2), x) == FiniteSet()
     assert solveset_real(3**(2 - x), x) == FiniteSet()
 
-    b = Symbol('b', positive=True)
-    y = Symbol('y', positive=True)
-    assert solveset_real(y - b*exp(a/x), x) == FiniteSet(a/log(y/b))
+    assert solveset_real(y - b*exp(a/x), x) == Intersection(S.Reals, FiniteSet(a/log(y/b)))
     # issue 4504
     assert solveset_real(2**x - 10, x) == FiniteSet(log(10)/log(2))
 
@@ -495,6 +500,7 @@ def test_no_sol():
 
 def test_sol_zero_real():
     assert solveset_real(0, x) == S.Reals
+    assert solveset(0, x, Interval(1, 2)) == Interval(1, 2)
     assert solveset_real(-x**2 - 2*x + (x + 1)**2 - 1, x) == S.Reals
 
 
@@ -566,11 +572,12 @@ def test_solve_abs():
         ) == Interval(-1, 2)
 
     # issue #10069
-    assert solveset_real(abs(1/(x - 1)) - 1 > 0, x) == \
-        ConditionSet(x, Eq((1 - Abs(x - 1))/Abs(x - 1) > 0, 0),
-            S.Reals)
-    assert solveset(abs(1/(x - 1)) - 1 > 0, x, domain=S.Reals
-        ) == Union(Interval.open(0, 1), Interval.open(1, 2))
+    eq = abs(1/(x - 1)) - 1 > 0
+    u = Union(Interval.open(0, 1), Interval.open(1, 2))
+    assert solveset_real(eq, x) == u
+    assert solveset(eq, x, domain=S.Reals) == u
+
+    raises(ValueError, lambda: solveset(abs(x) - 1, x))
 
 
 @XFAIL
@@ -620,16 +627,20 @@ def test_atan2():
 
 def test_piecewise():
     eq = Piecewise((x - 2, Gt(x, 2)), (2 - x, True)) - 3
-    f = Piecewise(((x - 2)**2, x >= 0), (0, True))
     assert set(solveset_real(eq, x)) == set(FiniteSet(-1, 5))
+
     absxm3 = Piecewise(
         (x - 3, S(0) <= x - 3),
-        (3 - x, S(0) > x - 3)
-    )
+        (3 - x, S(0) > x - 3))
     y = Symbol('y', positive=True)
     assert solveset_real(absxm3 - y, x) == FiniteSet(-y + 3, y + 3)
+
+    f = Piecewise(((x - 2)**2, x >= 0), (0, True))
     assert solveset(f, x, domain=S.Reals) == Union(FiniteSet(2), Interval(-oo, 0, True, True))
 
+    assert solveset(Piecewise((x + 1, x > 0), (I, True)) - I, x) == \
+        Interval(-oo, 0)
+
 
 def test_solveset_complex_polynomial():
     from sympy.abc import x, a, b, c
@@ -726,6 +737,16 @@ def test_solve_trig():
 
     assert solveset_real(sin(x)**2 + cos(x)**2, x) == S.EmptySet
 
+    assert solveset_complex(cos(x) - S.Half, x) == \
+        Union(imageset(Lambda(n, 2*n*pi + pi/3), S.Integers),
+              imageset(Lambda(n, 2*n*pi - pi/3), S.Integers))
+
+    y, a = symbols('y,a')
+    assert solveset(sin(y + a) - sin(y), a, domain=S.Reals) == \
+        Union(imageset(Lambda(n, 2*n*pi), S.Integers),
+        imageset(Lambda(n,
+        -I*(I*(2*n*pi +arg(-exp(-2*I*y))) + 2*im(y))), S.Integers))
+
 
 @XFAIL
 def test_solve_trig_abs():
@@ -739,12 +760,6 @@ def test_solve_invalid_sol():
     assert 0 not in solveset_complex((exp(x) - 1)/x, x)
 
 
-def test_solve_complex_unsolvable():
-    unsolved_object = ConditionSet(x, Eq(2*cos(x) - 1, 0), S.Complexes)
-    solution = solveset_complex(cos(x) - S.Half, x)
-    assert solution == unsolved_object
-
-
 @XFAIL
 def test_solve_trig_simplified():
     from sympy.abc import n
@@ -873,7 +888,7 @@ def test_solveset():
                                                   S.Integers)
 
 
-def test_conditonset():
+def test_conditionset():
     assert solveset(Eq(sin(x)**2 + cos(x)**2, 1), x, domain=S.Reals) == \
         ConditionSet(x, True, S.Reals)
 
@@ -889,10 +904,6 @@ def test_conditonset():
     assert solveset(x + sin(x) > 1, x, domain=S.Reals) == \
         ConditionSet(x, x + sin(x) > 1, S.Reals)
 
-    y,a = symbols('y,a')
-    assert solveset(sin(y + a) - sin(y), a, domain=S.Reals) == \
-        ConditionSet(a, Eq(-sin(y) + sin(y + a), 0), S.Reals)
-
 
 @XFAIL
 def test_conditionset_equality():
@@ -1070,3 +1081,23 @@ def test_issue_9913():
 
 def test_issue_10397():
     assert solveset(sqrt(x), x, S.Complexes) == FiniteSet(0)
+
+
+def test_simplification():
+    eq = x + (a - b)/(-2*a + 2*b)
+    assert solveset(eq, x) == FiniteSet(S.Half)
+    assert solveset(eq, x, S.Reals) == FiniteSet(S.Half)
+
+
+def test_issue_10555():
+    f = Function('f')
+    assert solveset(f(x) - pi/2, x, S.Reals) == \
+        ConditionSet(x, Eq(2*f(x) - pi, 0), S.Reals)
+
+
+def test_issue_8715():
+    eq = x + 1/x > -2 + 1/x
+    assert solveset(eq, x, S.Reals) == \
+        (Interval.open(-2, oo) - FiniteSet(0))
+    assert solveset(eq.subs(x,log(x)), x, S.Reals) == \
+        Interval.open(exp(-2), oo) - FiniteSet(1)
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index e2ae710a33..ae2cb74e29 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -41,7 +41,7 @@ def test_single_normal():
     pdf = density(Y)
     x = Symbol('x')
     assert (pdf(x) ==
-            2**S.Half*exp(-(x - mu)**2/(2*sigma**2))/(2*pi**S.Half*sigma))
+            2**S.Half*exp(-(mu - x)**2/(2*sigma**2))/(2*pi**S.Half*sigma))
 
     assert P(X**2 < 1) == erf(2**S.Half/2)