Repository: sympy/sympy

Issue:
AttributeError: 'list' object has no attribute 'xreplace'
The following code used to work in the previous version of sympy:
```
import sympy as sy
x, y = sy.Dummy(), sy.Dummy()
func = sy.lambdify([x, y], [y, x], 'sympy')

```

However, with the latest version (1.2), I get an error:

```
Traceback (most recent call last):
  File "/home/foo/.local/lib/python3.6/site-packages/sympy/utilities/lambdify.py", line 750, in _subexpr
    expr = sympify(expr).xreplace(dummies_dict)
AttributeError: 'list' object has no attribute 'xreplace'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "trial.py", line 3, in <module>
    func = sy.lambdify([x, y], [y, x], 'sympy')
  File "/home/foo/.local/lib/python3.6/site-packages/sympy/utilities/lambdify.py", line 456, in lambdify
    funcstr = funcprinter.doprint(funcname, args, expr)
  File "/home/foo/.local/lib/python3.6/site-packages/sympy/utilities/lambdify.py", line 661, in doprint
    argstrs, expr = self._preprocess(args, expr)
  File "/home/foo/.local/lib/python3.6/site-packages/sympy/utilities/lambdify.py", line 732, in _preprocess
    expr = self._subexpr(expr, {arg: dummy})
  File "/home/foo/.local/lib/python3.6/site-packages/sympy/utilities/lambdify.py", line 761, in _subexpr
    expr = [sub_expr(sympify(a), dummies_dict) for a in expr]
  File "/home/foo/.local/lib/python3.6/site-packages/sympy/utilities/lambdify.py", line 761, in <listcomp>
    expr = [sub_expr(sympify(a), dummies_dict) for a in expr]
NameError: name 'sub_expr' is not defined

```

Any idea what is causing this error?


Discussion:

Patch:
diff --git a/doc/src/modules/solvers/solveset.rst b/doc/src/modules/solvers/solveset.rst
index b82a6c9ae2..76894c7171 100644
--- a/doc/src/modules/solvers/solveset.rst
+++ b/doc/src/modules/solvers/solveset.rst
@@ -591,16 +591,6 @@ nonlinsolve
 .. autofunction:: sympy.solvers.solveset.nonlinsolve
 
 
-transolve
----------
-
-.. autofunction:: sympy.solvers.solveset._transolve
-
-.. autofunction:: sympy.solvers.solveset._is_exponential
-
-.. autofunction:: sympy.solvers.solveset._solve_expo
-
-
 Diophantine Equations (DEs)
 ---------------------------
 
diff --git a/sympy/discrete/__init__.py b/sympy/discrete/__init__.py
index 9870d8f534..443f52d6dc 100644
--- a/sympy/discrete/__init__.py
+++ b/sympy/discrete/__init__.py
@@ -3,9 +3,10 @@
 Transforms - fft, ifft, ntt, intt, fwht, ifwht,
     mobius_transform, inverse_mobius_transform
 Convolution - convolution, convolution_fft, convolution_ntt, convolution_fwht,
-    convolution_subset, covering_product
+    convolution_subset, covering_product, intersecting_product
+Recurrence - linrec
 """
 
 from .transforms import (fft, ifft, ntt, intt, fwht, ifwht,
     mobius_transform, inverse_mobius_transform)
-from .convolution import convolution, covering_product
+from .convolution import convolution
diff --git a/sympy/discrete/convolution.py b/sympy/discrete/convolution.py
index 98a2798632..d8a7f5d542 100644
--- a/sympy/discrete/convolution.py
+++ b/sympy/discrete/convolution.py
@@ -1,5 +1,6 @@
 """
-Convolution (using FFT, NTT, FWHT), Subset Convolution
+Convolution (using FFT, NTT, FWHT), Subset Convolution,
+Covering Product, Intersecting Product
 """
 from __future__ import print_function, division
 
@@ -7,20 +8,15 @@
 from sympy.core.compatibility import range, as_int, iterable
 from sympy.core.function import expand_mul
 from sympy.discrete.transforms import (
-    fft, ifft, ntt, intt, fwht, ifwht,
-    mobius_transform, inverse_mobius_transform)
+    fft, ifft, ntt, intt, fwht, ifwht)
 
 
-def convolution(a, b, cycle=0, dps=None, prime=None, dyadic=None, subset=None):
+def convolution(a, b, **hints):
     """
     Performs convolution by determining the type of desired
     convolution using hints.
 
-    Exactly one of `dps`, `prime`, `dyadic`, `subset` arguments should be
-    specified explicitly for identifying the type of convolution, and the
-    argument `cycle` can be specified optionally.
-
-    For the default arguments, linear convolution is performed using
+    If no hints are given, linear convolution is performed using
     FFT.
 
     Parameters
@@ -28,19 +24,22 @@ def convolution(a, b, cycle=0, dps=None, prime=None, dyadic=None, subset=None):
 
     a, b : iterables
         The sequences for which convolution is performed.
-    cycle : Integer
-        Specifies the length for doing cyclic convolution.
-    dps : Integer
-        Specifies the number of decimal digits for precision for
-        performing FFT on the sequence.
-    prime : Integer
-        Prime modulus of the form (m*2**k + 1) to be used for
-        performing NTT on the sequence.
-    dyadic : bool
-        Identifies the convolution type as dyadic (XOR)
-        convolution, which is performed using FWHT.
-    subset : bool
-        Identifies the convolution type as subset convolution.
+    hints : dict
+        Specifies the type of convolution to be performed.
+        The following hints can be given as keyword arguments.
+        dps : Integer
+            Specifies the number of decimal digits for precision for
+            performing FFT on the sequence.
+        prime : Integer
+            Prime modulus of the form (m*2**k + 1) to be used for
+            performing NTT on the sequence.
+        cycle : Integer
+            Specifies the length for doing cyclic convolution.
+        dyadic : bool
+            Identifies the convolution type as dyadic (XOR)
+            convolution, which is performed using FWHT.
+        subset : bool
+            Identifies the convolution type as subset convolution.
 
     Examples
     ========
@@ -70,18 +69,29 @@ def convolution(a, b, cycle=0, dps=None, prime=None, dyadic=None, subset=None):
 
     """
 
-    c = as_int(cycle)
+    fft = hints.pop('fft', None)
+    dps = hints.pop('dps', None)
+    p = hints.pop('prime', None)
+    c = as_int(hints.pop('cycle', 0))
+    dyadic = hints.pop('dyadic', None)
+    subset = hints.pop('subset', None)
+
     if c < 0:
         raise ValueError("The length for cyclic convolution must be non-negative")
 
+    fft = True if fft else None
     dyadic = True if dyadic else None
     subset = True if subset else None
-    if sum(x is not None for x in (prime, dps, dyadic, subset)) > 1:
-        raise TypeError("Ambiguity in determining the type of convolution")
+    if sum(x is not None for x in (p, dps, dyadic, subset)) > 1 or \
+            sum(x is not None for x in (fft, dyadic, subset)) > 1:
+        raise TypeError("Ambiguity in determining the convolution type")
 
-    if prime is not None:
-        ls = convolution_ntt(a, b, prime=prime)
-        return ls if not c else [sum(ls[i::c]) % prime for i in range(c)]
+    if p is not None:
+        ls = convolution_ntt(a, b, prime=p)
+        return ls if not c else [sum(ls[i::c]) % p for i in range(c)]
+
+    elif hints.pop('ntt', False):
+        raise TypeError("Prime modulus must be specified for performing NTT")
 
     if dyadic:
         ls = convolution_fwht(a, b)
@@ -345,73 +355,3 @@ def convolution_subset(a, b):
         c[mask] += expand_mul(a[smask] * b[mask^smask])
 
     return c
-
-
-#----------------------------------------------------------------------------#
-#                                                                            #
-#                              Covering Product                              #
-#                                                                            #
-#----------------------------------------------------------------------------#
-
-def covering_product(a, b):
-    """
-    Returns the covering product of given sequences.
-
-    The indices of each argument, considered as bit strings, correspond to
-    subsets of a finite set.
-
-    The covering product of given sequences is a sequence which contains
-    sum of products of the elements of the given sequences grouped by
-    bitwise OR of the corresponding indices.
-
-    The sequence is automatically padded to the right with zeros, as the
-    definition of subset based on bitmasks (indices) requires the size of
-    sequence to be a power of 2.
-
-    Parameters
-    ==========
-
-    a, b : iterables
-        The sequences for which covering product is to be obtained.
-
-    Examples
-    ========
-
-    >>> from sympy import symbols, S, I, covering_product
-    >>> u, v, x, y, z = symbols('u v x y z')
-
-    >>> covering_product([u, v], [x, y])
-    [u*x, u*y + v*x + v*y]
-    >>> covering_product([u, v, x], [y, z])
-    [u*y, u*z + v*y + v*z, x*y, x*z]
-
-    >>> covering_product([1, S(2)/3], [3, 4 + 5*I])
-    [3, 26/3 + 25*I/3]
-    >>> covering_product([1, 3, S(5)/7], [7, 8])
-    [7, 53, 5, 40/7]
-
-    References
-    ==========
-
-    .. [1] https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf
-
-    """
-
-    if not a or not b:
-        return []
-
-    a, b = a[:], b[:]
-    n = max(len(a), len(b))
-
-    if n&(n - 1): # not a power of 2
-        n = 2**n.bit_length()
-
-    # padding with zeros
-    a += [S.Zero]*(n - len(a))
-    b += [S.Zero]*(n - len(b))
-
-    a, b = mobius_transform(a), mobius_transform(b)
-    a = [expand_mul(x*y) for x, y in zip(a, b)]
-    a = inverse_mobius_transform(a)
-
-    return a
diff --git a/sympy/functions/combinatorial/factorials.py b/sympy/functions/combinatorial/factorials.py
index ec6ddea6d1..610a16e783 100644
--- a/sympy/functions/combinatorial/factorials.py
+++ b/sympy/functions/combinatorial/factorials.py
@@ -894,7 +894,7 @@ def eval(cls, n, k):
     def _eval_Mod(self, q):
         n, k = self.args
 
-        if any(x.is_integer is False for x in (n, k, q)):
+        if not all(x.is_integer for x in (n, k, q)):
             raise ValueError("Integers expected for binomial Mod")
 
         if all(x.is_Integer for x in (n, k, q)):
diff --git a/sympy/physics/units/prefixes.py b/sympy/physics/units/prefixes.py
index 42f1500ef8..184c6d44d7 100644
--- a/sympy/physics/units/prefixes.py
+++ b/sympy/physics/units/prefixes.py
@@ -29,7 +29,6 @@ class Prefix(Expr):
       class).
     """
     _op_priority = 13.0
-    is_commutative = True
 
     def __new__(cls, name, abbrev, exponent, base=sympify(10)):
 
diff --git a/sympy/physics/units/util.py b/sympy/physics/units/util.py
index 9b520c34ee..cc15b65393 100644
--- a/sympy/physics/units/util.py
+++ b/sympy/physics/units/util.py
@@ -9,7 +9,7 @@
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 
 from sympy import Add, Function, Mul, Pow, Rational, Tuple, sympify
-from sympy.core.compatibility import reduce, Iterable, ordered
+from sympy.core.compatibility import reduce, Iterable
 from sympy.physics.units.dimensions import Dimension, dimsys_default
 from sympy.physics.units.quantities import Quantity
 from sympy.physics.units.prefixes import Prefix
@@ -131,46 +131,30 @@ def get_total_scale_factor(expr):
 
 
 def quantity_simplify(expr):
-    """Return an equivalent expression in which prefixes are replaced
-    with numerical values and products of related quantities are
-    unified in a canonical manner.
-
-    Examples
-    ========
-
-    >>> from sympy.physics.units.util import quantity_simplify
-    >>> from sympy.physics.units.prefixes import kilo
-    >>> from sympy.physics.units import foot, inch
-    >>> quantity_simplify(kilo*foot*inch)
-    250*foot**2/3
-    """
-
     if expr.is_Atom:
         return expr
-
-    if isinstance(expr, Prefix):
-        return expr.scale_factor
-
-    expr = expr.func(*map(quantity_simplify, expr.args))
-
-    if not expr.is_Mul or not expr.has(Quantity):
-        return expr
-
-    args_pow = [arg.as_base_exp() for arg in expr.args]
-    quantity_pow, other_pow = sift(
-        args_pow, lambda x: isinstance(x[0], Quantity), binary=True)
-    coeff = Mul.fromiter([
-        Pow(b, e, evaluate=False) for b, e in other_pow])
+    if not expr.is_Mul:
+        return expr.func(*map(quantity_simplify, expr.args))
+
+    if expr.has(Prefix):
+        coeff, args = expr.as_coeff_mul(Prefix)
+        args = list(args)
+        for arg in args:
+            if isinstance(arg, Pow):
+                coeff = coeff * (arg.base.scale_factor ** arg.exp)
+            else:
+                coeff = coeff * arg.scale_factor
+        expr = coeff
+
+    coeff, args = expr.as_coeff_mul(Quantity)
+    args_pow = [arg.as_base_exp() for arg in args]
+    quantity_pow, other_pow = sift(args_pow, lambda x: isinstance(x[0], Quantity), binary=True)
     quantity_pow_by_dim = sift(quantity_pow, lambda x: x[0].dimension)
-    new_quantities = []
-    for _, bp in quantity_pow_by_dim.items():
-        if len(bp) == 1:
-            new_quantities.append(bp[0][0]**bp[0][1])
-        else:
-            # just let reference quantity be the first quantity,
-            # picked from an ordered list
-            bp = list(ordered(bp))
-            ref = bp[0][0]/bp[0][0].scale_factor
-            new_quantities.append(Mul.fromiter((
-                (ref*b.scale_factor)**p for b, p in bp)))
-    return coeff*Mul.fromiter(new_quantities)
+    # Just pick the first quantity:
+    ref_quantities = [i[0][0] for i in quantity_pow_by_dim.values()]
+    new_quantities = [
+        Mul.fromiter(
+            (quantity*i.scale_factor/quantity.scale_factor)**p for i, p in v)
+            if len(v) > 1 else v[0][0]**v[0][1]
+        for quantity, (k, v) in zip(ref_quantities, quantity_pow_by_dim.items())]
+    return coeff*Mul.fromiter(other_pow)*Mul.fromiter(new_quantities)
diff --git a/sympy/plotting/plot.py b/sympy/plotting/plot.py
index a890927f81..bb36d16ded 100644
--- a/sympy/plotting/plot.py
+++ b/sympy/plotting/plot.py
@@ -819,8 +819,9 @@ def get_meshes(self):
 ### Contours
 class ContourSeries(BaseSeries):
     """Representation for a contour plot."""
-    # The code is mostly repetition of SurfaceOver2DRange.
-    # Presently used in contour_plot function
+    #The code is mostly repetition of SurfaceOver2DRange.
+    #XXX: Presently not used in any of those functions.
+    #XXX: Add contour plot and use this seties.
 
     is_contour = True
 
@@ -1731,84 +1732,6 @@ def plot3d_parametric_surface(*args, **kwargs):
         plots.show()
     return plots
 
-def plot_contour(*args, **kwargs):
-    """
-    Draws contour plot of a function
-
-    Usage
-    =====
-
-    Single plot
-
-    ``plot_contour(expr, range_x, range_y, **kwargs)``
-
-    If the ranges are not specified, then a default range of (-10, 10) is used.
-
-    Multiple plot with the same range.
-
-    ``plot_contour(expr1, expr2, range_x, range_y, **kwargs)``
-
-    If the ranges are not specified, then a default range of (-10, 10) is used.
-
-    Multiple plots with different ranges.
-
-    ``plot_contour((expr1, range_x, range_y), (expr2, range_x, range_y), ..., **kwargs)``
-
-    Ranges have to be specified for every expression.
-
-    Default range may change in the future if a more advanced default range
-    detection algorithm is implemented.
-
-    Arguments
-    =========
-
-    ``expr`` : Expression representing the function along x.
-
-    ``range_x``: (x, 0, 5), A 3-tuple denoting the range of the x
-    variable.
-
-    ``range_y``: (y, 0, 5), A 3-tuple denoting the range of the y
-     variable.
-
-    Keyword Arguments
-    =================
-
-    Arguments for ``ContourSeries`` class:
-
-    ``nb_of_points_x``: int. The x range is sampled uniformly at
-    ``nb_of_points_x`` of points.
-
-    ``nb_of_points_y``: int. The y range is sampled uniformly at
-    ``nb_of_points_y`` of points.
-
-    Aesthetics:
-
-    ``surface_color``: Function which returns a float. Specifies the color for
-    the surface of the plot. See ``sympy.plotting.Plot`` for more details.
-
-    If there are multiple plots, then the same series arguments are applied to
-    all the plots. If you want to set these options separately, you can index
-    the returned ``Plot`` object and set it.
-
-    Arguments for ``Plot`` class:
-
-    ``title`` : str. Title of the plot.
-
-    See Also
-    ========
-    Plot, ContourSeries
-    """
-
-    args = list(map(sympify, args))
-    show = kwargs.pop('show', True)
-    plot_expr = check_arguments(args, 1, 2)
-    series = [ContourSeries(*arg) for arg in plot_expr]
-    plot_contours = Plot(*series, **kwargs)
-    if len(plot_expr[0].free_symbols) > 2:
-        raise ValueError('Contour Plot cannot Plot for more than two variables.')
-    if show:
-        plot_contours.show()
-    return plot_contours
 
 def check_arguments(args, expr_len, nb_of_free_symbols):
     """
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 4c8a1cad69..b310bda4ae 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -3,9 +3,6 @@
 
     - solve a single equation for a single variable, in any domain either real or complex.
 
-    - solve a single transcendental equation for a single variable in any domain either real or complex.
-      (currently supports solving in real domain only)
-
     - solve a system of linear equations with N variables and M equations.
 
     - solve a system of Non Linear Equations with N variables and M equations
@@ -13,17 +10,14 @@
 from __future__ import print_function, division
 
 from sympy.core.sympify import sympify
-from sympy.core import (S, Pow, Dummy, pi, Expr, Wild, Mul, Equality,
-                        Add)
+from sympy.core import S, Pow, Dummy, pi, Expr, Wild, Mul, Equality
 from sympy.core.containers import Tuple
 from sympy.core.facts import InconsistentAssumptions
 from sympy.core.numbers import I, Number, Rational, oo
-from sympy.core.function import (Lambda, expand_complex, AppliedUndef,
-                                expand_log)
+from sympy.core.function import (Lambda, expand_complex, AppliedUndef, Function)
 from sympy.core.relational import Eq
 from sympy.core.symbol import Symbol
 from sympy.simplify.simplify import simplify, fraction, trigsimp
-from sympy.simplify import powdenest
 from sympy.functions import (log, Abs, tan, cot, sin, cos, sec, csc, exp,
                              acos, asin, acsc, asec, arg,
                              piecewise_fold, Piecewise)
@@ -36,7 +30,7 @@
 from sympy.sets.sets import Set
 from sympy.matrices import Matrix
 from sympy.polys import (roots, Poly, degree, together, PolynomialError,
-                         RootOf, factor)
+                         RootOf)
 from sympy.solvers.solvers import (checksol, denoms, unrad,
     _simple_dens, recast_to_symbols)
 from sympy.solvers.polysys import solve_poly_system
@@ -238,16 +232,16 @@ def _invert_real(f, g_ys, symbol):
 
         if not base_has_sym:
             rhs = g_ys.args[0]
-            if base.is_positive:
+            if base > S.Zero:
                 return _invert_real(expo,
-                    imageset(Lambda(n, log(n, base, evaluate=False)), g_ys), symbol)
-            elif base.is_negative:
+                    imageset(Lambda(n, log(n)/log(base)), g_ys), symbol)
+            elif base < S.Zero:
                 from sympy.core.power import integer_log
                 s, b = integer_log(rhs, base)
                 if b:
-                    return _invert_real(expo, FiniteSet(s), symbol)
+                    return expo, FiniteSet(s)
             elif rhs is S.One:
-                # special case: 0**x - 1
+                #special case: 0**x - 1
                 return (expo, FiniteSet(0))
             return (expo, S.EmptySet)
 
@@ -906,6 +900,7 @@ def _solveset(f, symbol, domain, _check=False):
             solns = solver(expr, symbol, in_set)
             result += solns
     elif isinstance(f, Eq):
+        from sympy.core import Add
         result = solver(Add(f.lhs, - f.rhs, evaluate=False), symbol, domain)
 
     elif f.is_Relational:
@@ -948,12 +943,7 @@ def _solveset(f, symbol, domain, _check=False):
                     elif equation.has(Abs):
                         result += _solve_abs(f, symbol, domain)
                     else:
-                        result_rational = _solve_as_rational(equation, symbol, domain)
-                        if isinstance(result_rational, ConditionSet):
-                            # may be a transcendental type equation
-                            result += _transolve(equation, symbol, domain)
-                        else:
-                            result += result_rational
+                        result += _solve_as_rational(equation, symbol, domain)
                 else:
                     result += solver(equation, symbol)
 
@@ -988,386 +978,6 @@ def _solveset(f, symbol, domain, _check=False):
     return result
 
 
-def _term_factors(f):
-    """
-    Iterator to get the factors of all terms present
-    in the given equation.
-
-    Parameters
-    ==========
-
-    f : Expr
-        Equation that needs to be addressed
-
-    Returns
-    =======
-
-    Factors of all terms present in the equation.
-
-    Examples
-    ========
-
-    >>> from sympy import symbols
-    >>> from sympy.solvers.solveset import _term_factors
-    >>> x = symbols('x')
-    >>> list(_term_factors(-2 - x**2 + x*(x + 1)))
-    [-2, -1, x**2, x, x + 1]
-    """
-    for add_arg in Add.make_args(f):
-        for mul_arg in Mul.make_args(add_arg):
-            yield mul_arg
-
-
-def _solve_expo(f, symbol):
-    r"""
-    Helper function for solving (supported) exponential equations.
-
-    Parameters
-    ==========
-
-    f : Expr
-        The exponential equation to be solved
-
-    symbol : Symbol
-        The variable in which the equation is solved
-
-    Returns
-    =======
-
-    An equation in `log` form that `solveset` might better handle.
-
-    `None`:
-        If the equation is not of supported type
-
-    Notes
-    =====
-
-    Exponential equations are the sum of (currently) at most
-    two terms with one or both of them having a power with a
-    symbol-dependent exponent.
-
-    For example
-
-    .. math:: 5^{2x + 3} - 5^{3x - 1}
-
-    .. math:: 4^{5 - 9x} - e^{2 - x}
-
-    Examples
-    ========
-
-    >>> from sympy.solvers.solveset import _solve_expo as solve_expo
-    >>> from sympy import symbols
-    >>> x = symbols('x', real=True)
-    >>> solve_expo(3**(2*x) - 2**(x + 3), x)
-    2*x*log(3) - (x + 3)*log(2)
-
-    * Proof of correctness of the method
-
-    The logarithm function is the inverse of the exponential function.
-    The defining relation between exponentiation and logarithm is:
-
-    .. math:: {log_b x} = y \enspace if \enspace b^y = x
-
-    Therefore if we are given an equation with exponent terms, we can
-    convert every term to its corresponding log form. This is achieved by
-    taking logarithms and expanding the equation using log identities
-    so that it can easily be handled by `solveset`.
-
-    For example:
-
-    .. math:: 3^{2x} = 2^{x + 3}
-
-    Taking log both sides will reduce the equation to
-
-    .. math:: (2x)*log(3) = (x + 3)*log(2)
-
-    This form can be easily handed by `solveset`.
-    """
-
-    if not (isinstance(f, Add) and len(f.args) == 2):
-        # solving for the sum of more than two powers is possible
-        # but not yet implemented
-        return None
-
-    lhs, rhs = list(ordered(f.args))
-    log_type_equation = expand_log(log(lhs)) - expand_log(log(-rhs))
-
-    return log_type_equation
-
-
-def _is_exponential(f, symbol):
-    r"""
-    Helper to check whether an equation is exponential or not.
-
-    Parameters
-    ==========
-
-    f : Expr
-        The equation to be checked
-
-    symbol : Symbol
-        The variable in which the equation is checked
-
-    Returns
-    =======
-
-    `True` if the equation is in exponential form otherwise `False`.
-
-    Examples
-    ========
-
-    >>> from sympy import symbols, cos, exp
-    >>> from sympy.solvers.solveset import _is_exponential as check
-    >>> x, y = symbols('x y')
-    >>> check(x**y - x, y)
-    True
-    >>> check(x**y - x, x)
-    False
-    >>> check(exp(x + 3), x)
-    True
-    >>> check(cos(2**x), x)
-    False
-
-    * Philosophy behind the helper
-
-    The function extracts each term of the equation and checks if it is
-    of exponential form w.r.t `symbol`.
-    """
-
-    expr_args = _term_factors(f)
-    for expr_arg in expr_args:
-        if isinstance(expr_arg, (Pow, exp)) and (
-                symbol in expr_arg.exp.free_symbols):
-            return True
-    return False
-
-
-def _transolve(f, symbol, domain):
-    r"""
-    Function to solve transcendental equations. It is a helper to
-    `solveset` and should be used internally. `\_transolve`
-    currently supports the following class of equations:
-
-        - Exponential equations
-
-    Parameters
-    ==========
-
-    f : Any transcendental equation that needs to be solved.
-        This needs to be an expression, which is assumed
-        to be equal to 0.
-
-    symbol : The variable for which the equation is solved.
-        This needs to be of class `Symbol`.
-
-    domain : A set over which the equation is solved.
-        This needs to be of class `Set`.
-
-    Returns
-    =======
-
-    Set
-        A set of values for `symbol` for which `f` is equal to
-        zero. An `EmptySet` is returned if `f` does not have solutions
-        in respective domain. A `ConditionSet` is returned as unsolved
-        object if algorithms to evaluate complete solution are not
-        yet implemented.
-
-    How to use `\_transolve`
-    ========================
-
-    `\_transolve` should not be used as an independent function, because
-    it assumes that the equation (`f`) and the `symbol` comes from
-    `solveset` and might have undergone a few modification(s).
-    To use `\_transolve` as an independent function the equation (`f`)
-    and the `symbol` should be passed as they would have been by
-    `solveset`.
-
-    Examples
-    ========
-
-    >>> from sympy.solvers.solveset import _transolve as transolve
-    >>> from sympy.solvers.solvers import _tsolve as tsolve
-    >>> from sympy import symbols, S, pprint
-    >>> x = symbols('x', real=True) # assumption added
-    >>> transolve(5**(x - 3) - 3**(2*x + 1), x, S.Reals)
-    {-log(375)/(-log(5) + 2*log(3))}
-
-    How `\_transolve` works
-    =======================
-
-    `\_transolve` uses two types of helper functions to solve equations
-    of a particular class:
-
-    Identifying helpers: To determine whether a given equation
-    belongs to a certain class of equation or not. Returns either
-    True or False.
-
-    Solving helpers: Once an equation is identified, a corresponding
-    helper either solves the equation or returns a form of the equation
-    that `solveset` might better be able to handle.
-
-    * Philosophy behind the module
-
-    The purpose of `\_transolve` is to take equations which are not
-    already polynomial in their generator(s) and to either recast them
-    as such through a valid transformation or to solve them outright.
-    A pair of helper functions for each class of supported
-    transcendental functions are employed for this purpose. One
-    identifies the transcendental form of an equation and the other
-    either solves it or recasts it into a tractable form that can be
-    solved by  `solveset`.
-    For example, an equation in the form `a*b**f(x) - c*d**g(x) = 0`
-    can be transformed to
-    `log(a) + f(x)*log(b) - log(c) - g(x)*log(d) = 0`
-    (under certain assumptions) and this can be solved with `solveset`
-    if `f(x)` and `g(x)` are in polynomial form.
-
-    How `\_transolve` is better than `\_tsolve`
-    ===========================================
-
-    1) Better output
-
-    `\_transolve` provides expressions in a more simplified form.
-
-    Consider a simple exponential equation
-
-    >>> f = 3**(2*x) - 2**(x + 3)
-    >>> pprint(transolve(f, x, S.Reals), use_unicode=False)
-        -3*log(2)
-    {------------------}
-     -2*log(3) + log(2)
-    >>> pprint(tsolve(f, x), use_unicode=False)
-         /   3     \
-         | --------|
-         | log(2/9)|
-    [-log\2         /]
-
-    2) Extensible
-
-    The API of `\_transolve` is designed such that it is easily
-    extensible, i.e. the code that solves a given class of
-    equations is encapsulated in a helper and not mixed in with
-    the code of \_transolve itself.
-
-    3) Modular
-
-    `\_transolve` is designed to be modular i.e, for every class of
-    equation a separate helper for identification and solving is
-    implemented. This makes it easy to change or modify any of the
-    method implemented directly in the helpers without interfering
-    with the actual structure of the API.
-
-    4) Faster Computation
-
-    Solving equation via `\_transolve` is much faster as compared to
-    `\_tsolve`. In `solve` attempts are made computing every possibility
-    to get the solutions. This series of attempts makes solving a bit
-    slow. Whereas in `\_transolve` computation begins only when the
-    equation is identified of being a particular type.
-
-    How to add new class of equations
-    =================================
-
-    Adding a new class of equation solver is a three-step procedure:
-
-    - Identify the type of the equations
-
-      Determine the type of the class of equations to which they belong,
-      it could be of `Add`, `Pow`, etc. types. Separate internal functions
-      are used for each type, if it is already present include
-      identification and solving helpers within the routine otherwise add
-      a new internal function and then include identification and solving
-      helpers within it. Something like:
-
-      .. code-block:: python
-
-        def add_type(eq, x):
-            ....
-            if _is_exponential(eq, x):
-                new_eq = _solve_expo(eq, x)
-        ....
-        if eq.is_Add:
-            result = add_type(eq, x)
-
-    - Define the identification helper.
-
-    - Define the solving helper.
-
-    Apart from this, a few other things needs to be taken care while
-    adding an equation solver:
-
-    - Naming conventions:
-      Name of the *identification* *helper* should be as
-      `\_is\_class` where `class` will be the name or abbreviation
-      of the class of equation. The *solving* *helper* will be named as
-      `\_solve\_class`.
-      For example: for exponential equations it becomes
-      `\_is\_exponential` and `\_solve\_expo`.
-    - The helpers should take two input parameters, the equation to be
-      checked and the variable for which a solution is being sought.
-    - Be sure to consider corner cases.
-    - Add tests for each helper.
-    - Add a docstring to your helper that describes the method
-      implemented.
-      The following things needs to be included while writing the
-      documentation for the helpers:
-
-      - What is the purpose of the helper.
-      - How it solves and identifies the equation.
-      - Examples should be included with its proof of correctness.
-      - What does the helper returns.
-    """
-
-    def add_type(eq, symbol, domain):
-        """
-        Helper for `_transolve` to handle equations of
-        `Add` type, i.e. equations taking the form as
-        `a*f(x) + b*g(x) + .... + c = 0`.
-        For example: 4**x + 8**x = 0
-        """
-        result = ConditionSet(symbol, Eq(eq, 0), domain)
-        new_eq = None
-        # try factoring the equation;
-        # powdenest is used to try to get powers in the standard form
-        # for better factoring
-        simplified_equation = factor(powdenest(eq))
-
-        if simplified_equation.is_Mul:
-            new_eq = simplified_equation
-
-        # check if it is exponential type equation
-        elif _is_exponential(simplified_equation, symbol):
-            new_eq = _solve_expo(simplified_equation, symbol)
-
-        if new_eq:
-            result = _solveset(new_eq, symbol, domain)
-
-        return result
-
-    unsolved_result = ConditionSet(symbol, Eq(f, 0), domain)
-    result = unsolved_result
-
-    # invert_complex handles the call to the desired inverter based
-    # on the domain specified.
-    lhs, rhs_s = invert_complex(f, 0, symbol, domain)
-
-    if isinstance(rhs_s, FiniteSet):
-        assert (len(rhs_s.args)) == 1
-        rhs = rhs_s.args[0]
-        equation = Add(lhs, -rhs, evaluate=False)
-
-        if lhs.is_Add:
-            result = add_type(equation, symbol, domain)
-    else:
-        result = rhs_s
-
-    if isinstance(result, ConditionSet):
-        result = unsolved_result
-
-    return result
-
-
 def solveset(f, symbol=None, domain=S.Complexes):
     r"""Solves a given inequality or equation with set as output
 
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index c297cc2208..7140512739 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -134,8 +134,7 @@ def ContinuousRV(symbol, density, set=Interval(-oo, oo)):
     >>> P(X>0)
     1/2
     """
-    pdf = Piecewise((density, set.as_relational(symbol)), (0, True))
-    pdf = Lambda(symbol, pdf)
+    pdf = Lambda(symbol, density)
     dist = ContinuousDistributionHandmade(pdf, set)
     return SingleContinuousPSpace(symbol, dist).value
 
diff --git a/sympy/stats/joint_rv.py b/sympy/stats/joint_rv.py
index fdb247d0c8..f82a980fba 100644
--- a/sympy/stats/joint_rv.py
+++ b/sympy/stats/joint_rv.py
@@ -18,7 +18,8 @@
 from sympy.integrals.integrals import Integral, integrate
 from sympy.stats.rv import (ProductPSpace, NamedArgsMixin,
      ProductDomain, RandomSymbol)
-from sympy.matrices import ImmutableMatrix
+from sympy.core.containers import Tuple
+
 class JointPSpace(ProductPSpace):
     """
     Represents a joint probability space. Represented using symbols for
@@ -89,7 +90,7 @@ def __new__(cls, *args):
         args = list(map(sympify, args))
         for i in range(len(args)):
             if isinstance(args[i], list):
-                args[i] = ImmutableMatrix(args[i])
+                args[i] = Tuple.fromiter(j for j in args[i])
         return Basic.__new__(cls, *args)
 
     @property
diff --git a/sympy/stats/joint_rv_types.py b/sympy/stats/joint_rv_types.py
index a0c6b3d5c7..73af5a8031 100644
--- a/sympy/stats/joint_rv_types.py
+++ b/sympy/stats/joint_rv_types.py
@@ -1,7 +1,7 @@
 from sympy import sympify, S, pi, sqrt, exp, Lambda, Indexed, Symbol, Gt
 from sympy.stats.rv import _value_check
 from sympy.stats.joint_rv import JointDistribution, JointPSpace
-from sympy.matrices import ImmutableMatrix
+from sympy.matrices.dense import Matrix
 from sympy.matrices.expressions.determinant import det
 
 # __all__ = ['MultivariateNormal',
@@ -14,7 +14,6 @@ def multivariate_rv(cls, sym, *args):
     sym = sympify(sym)
     args = list(map(sympify, args))
     dist = cls(*args)
-    args = dist.args
     dist.check(*args)
     return JointPSpace(sym, dist).value
 
@@ -32,6 +31,7 @@ def set(self):
         return S.Reals**k
 
     def check(self, mu, sigma):
+        mu, sigma = Matrix([mu]), Matrix(sigma)
         _value_check(len(mu) == len(sigma.col(0)),
             "Size of the mean vector and covariance matrix are incorrect.")
         #check if covariance matrix is positive definite or not.
@@ -39,23 +39,23 @@ def check(self, mu, sigma):
             "The covariance matrix must be positive definite. ")
 
     def pdf(self, *args):
-        mu, sigma = self.mu, self.sigma
+        mu, sigma = Matrix(self.mu), Matrix(self.sigma)
         k = len(mu)
-        args = ImmutableMatrix(args)
+        args = Matrix(args)
         x = args - mu
         return  S(1)/sqrt((2*pi)**(k)*det(sigma))*exp(
             -S(1)/2*x.transpose()*(sigma.inv()*\
                 x))[0]
 
     def marginal_distribution(self, indices, sym):
-        sym = ImmutableMatrix([Symbol(str(Indexed(sym, i))) for i in indices])
-        _mu, _sigma = self.mu, self.sigma
+        sym = Matrix([Symbol(str(Indexed(sym, i))) for i in indices])
+        _mu, _sigma = Matrix(self.mu), Matrix(self.sigma)
         k = len(self.mu)
         for i in range(k):
             if i not in indices:
-                _mu = _mu.row_del(i)
-                _sigma = _sigma.col_del(i)
-                _sigma = _sigma.row_del(i)
+                _mu.row_del(i)
+                _sigma.col_del(i)
+                _sigma.row_del(i)
         return Lambda(sym, S(1)/sqrt((2*pi)**(len(_mu))*det(_sigma))*exp(
             -S(1)/2*(_mu - sym).transpose()*(_sigma.inv()*\
                 (_mu - sym)))[0])
@@ -73,6 +73,7 @@ def set(self):
         return S.Reals**k
 
     def check(self, mu, sigma):
+        mu, sigma = Matrix([mu]), Matrix(sigma)
         _value_check(len(mu) == len(sigma.col(0)),
             "Size of the mean vector and covariance matrix are incorrect.")
         #check if covariance matrix is positive definite or not.
@@ -81,11 +82,11 @@ def check(self, mu, sigma):
 
     def pdf(self, *args):
         from sympy.functions.special.bessel import besselk
-        mu, sigma = self.mu, self.sigma
+        mu, sigma = Matrix(self.mu), Matrix(self.sigma)
         mu_T = mu.transpose()
         k = S(len(mu))
         sigma_inv = sigma.inv()
-        args = ImmutableMatrix(args)
+        args = Matrix(args)
         args_T = args.transpose()
         x = (mu_T*sigma_inv*mu)[0]
         y = (args_T*sigma_inv*args)[0]
@@ -108,6 +109,7 @@ def set(self):
         return S.Reals**k
 
     def check(self, mu, sigma, v):
+        mu, sigma = Matrix([mu]), Matrix(sigma)
         _value_check(len(mu) == len(sigma.col(0)),
             "Size of the location vector and shape matrix are incorrect.")
         #check if covariance matrix is positive definite or not.
@@ -116,11 +118,11 @@ def check(self, mu, sigma, v):
 
     def pdf(self, *args):
         from sympy.functions.special.gamma_functions import gamma
-        mu, sigma = self.mu, self.shape_mat
+        mu, sigma = Matrix(self.mu), Matrix(self.shape_mat)
         v = S(self.dof)
         k = S(len(mu))
         sigma_inv = sigma.inv()
-        args = ImmutableMatrix(args)
+        args = Matrix(args)
         x = args - mu
         return gamma((k + v)/2)/(gamma(v/2)*(v*pi)**(k/2)*sqrt(det(sigma)))\
         *(1 + 1/v*(x.transpose()*sigma_inv*x)[0])**((-v - k)/2)
diff --git a/sympy/tensor/index_methods.py b/sympy/tensor/index_methods.py
index 947c5b5ad9..808466e100 100644
--- a/sympy/tensor/index_methods.py
+++ b/sympy/tensor/index_methods.py
@@ -12,57 +12,38 @@
 
 from __future__ import print_function, division
 
-
 from sympy.core.function import Function
 from sympy.functions import exp, Piecewise
 from sympy.tensor.indexed import Idx, Indexed
 
 
 from sympy.core.compatibility import reduce
-from collections import OrderedDict
-from sympy.utilities import sift
-class IndexConformanceException(Exception):
-    pass
-
-def _unique_and_repeated(inds):
-    """
-    Returns the unique and repeated indices. Also note, from the examples given below
-    that the order of indices is maintained as given in the input.
 
-    Examples :
-    ========
 
-    >>> from sympy.tensor.index_methods import _unique_and_repeated
+class IndexConformanceException(Exception):
+    pass
 
-    >>> _unique_and_repeated([2, 3, 1, 3, 0, 4, 0])
-    ([2, 1, 4], [3, 0])
-    """
-    uniq = OrderedDict()
-    for i in inds:
-        if i in uniq:
-            uniq[i] = 0
-        else:
-            uniq[i] = 1
-    return sift(uniq, lambda x: uniq[x], binary=True)
 
 def _remove_repeated(inds):
-    """
-    Removes repeated objects from sequences
+    """Removes repeated objects from sequences
 
     Returns a set of the unique objects and a tuple of all that have been
     removed.
 
-    Examples :
-    ========
-
     >>> from sympy.tensor.index_methods import _remove_repeated
     >>> l1 = [1, 2, 3, 2]
     >>> _remove_repeated(l1)
     ({1, 3}, (2,))
 
     """
-    u, r = _unique_and_repeated(inds)
-    return set(u), tuple(r)
+    sum_index = {}
+    for i in inds:
+        if i in sum_index:
+            sum_index[i] += 1
+        else:
+            sum_index[i] = 0
+    inds = [x for x in inds if not sum_index[x]]
+    return set(inds), tuple([ i for i in sum_index if sum_index[i] ])
 
 
 def _get_indices_Mul(expr, return_dummies=False):
diff --git a/sympy/utilities/lambdify.py b/sympy/utilities/lambdify.py
index e75f8c7867..36606386be 100644
--- a/sympy/utilities/lambdify.py
+++ b/sympy/utilities/lambdify.py
@@ -741,8 +741,7 @@ def _preprocess(self, args, expr):
 
         return argstrs, expr
 
-    @staticmethod
-    def _subexpr(expr, dummies_dict):
+    def _subexpr(self, expr, dummies_dict):
         from sympy.matrices import DeferredVector
         from sympy import sympify
 
@@ -752,13 +751,13 @@ def _subexpr(expr, dummies_dict):
             if isinstance(expr, DeferredVector):
                 pass
             elif isinstance(expr, dict):
-                k = [sub_expr(sympify(a), dummies_dict) for a in expr.keys()]
-                v = [sub_expr(sympify(a), dummies_dict) for a in expr.values()]
+                k = [self._subexpr(sympify(a), dummies_dict) for a in expr.keys()]
+                v = [self._subexpr(sympify(a), dummies_dict) for a in expr.values()]
                 expr = dict(zip(k, v))
             elif isinstance(expr, tuple):
-                expr = tuple(sub_expr(sympify(a), dummies_dict) for a in expr)
+                expr = tuple(self._subexpr(sympify(a), dummies_dict) for a in expr)
             elif isinstance(expr, list):
-                expr = [sub_expr(sympify(a), dummies_dict) for a in expr]
+                expr = [self._subexpr(sympify(a), dummies_dict) for a in expr]
         return expr
 
     def _print_funcargwrapping(self, args):


Test patch:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index 60e4bdbd72..4158202230 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -1458,8 +1458,7 @@ def test_sympy__stats__joint_rv__JointDistribution():
 
 def test_sympy__stats__joint_rv_types__MultivariateNormalDistribution():
     from sympy.stats.joint_rv_types import MultivariateNormalDistribution
-    assert _test_args(
-        MultivariateNormalDistribution([0, 1], [[1, 0],[0, 1]]))
+    assert _test_args(MultivariateNormalDistribution([0, 1], [[1, 0],[0, 1]]))
 
 def test_sympy__stats__joint_rv_types__MultivariateLaplaceDistribution():
     from sympy.stats.joint_rv_types import MultivariateLaplaceDistribution
diff --git a/sympy/discrete/tests/test_convolution.py b/sympy/discrete/tests/test_convolution.py
index 746f85e98d..1b1bf6d568 100644
--- a/sympy/discrete/tests/test_convolution.py
+++ b/sympy/discrete/tests/test_convolution.py
@@ -5,7 +5,7 @@
 from sympy.core.compatibility import range
 from sympy.discrete.convolution import (
     convolution, convolution_fft, convolution_ntt, convolution_fwht,
-    convolution_subset, covering_product)
+    convolution_subset)
 from sympy.utilities.pytest import raises
 from sympy.abc import x, y
 
@@ -16,9 +16,9 @@ def test_convolution():
     c = [3, 5, 3, 7, 8]
     d = [1422, 6572, 3213, 5552]
 
-    assert convolution(a, b) == convolution_fft(a, b)
-    assert convolution(a, b, dps=9) == convolution_fft(a, b, dps=9)
-    assert convolution(a, d, dps=7) == convolution_fft(d, a, dps=7)
+    assert convolution(a, b, fft=True) == convolution_fft(a, b)
+    assert convolution(a, b, dps=9, fft=True) == convolution_fft(a, b, dps=9)
+    assert convolution(a, d, fft=True, dps=7) == convolution_fft(d, a, dps=7)
     assert convolution(a, d[1:], dps=3) == convolution_fft(d[1:], a, dps=3)
 
     # prime moduli of the form (m*2**k + 1), sequence length
@@ -27,26 +27,35 @@ def test_convolution():
     q = 19*2**10 + 1
 
     # ntt
-    assert convolution(d, b, prime=q) == convolution_ntt(b, d, prime=q)
+    assert convolution(d, b, ntt=True, prime=q) == convolution_ntt(b, d, prime=q)
     assert convolution(c, b, prime=p) == convolution_ntt(b, c, prime=p)
-    assert convolution(d, c, prime=p) == convolution_ntt(c, d, prime=p)
+    assert convolution(d, c, prime=p, ntt=True) == convolution_ntt(c, d, prime=p)
+    raises(TypeError, lambda: convolution(b, d, ntt=True))
+    raises(TypeError, lambda: convolution(b, d, ntt=True, cycle=0))
     raises(TypeError, lambda: convolution(b, d, dps=5, prime=q))
-    raises(TypeError, lambda: convolution(b, d, dps=6, prime=q))
+    raises(TypeError, lambda: convolution(b, d, dps=6, ntt=True, prime=q))
+    raises(TypeError, lambda: convolution(b, d, fft=True, dps=7, ntt=True, prime=q))
+    # ntt is a specialized variant of fft, TypeError should not be raised
+    assert convolution(b, d, fft=True, ntt=True, prime=q) == \
+            convolution_ntt(b, d, prime=q)
 
     # fwht
     assert convolution(a, b, dyadic=True) == convolution_fwht(a, b)
     assert convolution(a, b, dyadic=False) == convolution(a, b)
-    raises(TypeError, lambda: convolution(b, d, dps=2, dyadic=True))
-    raises(TypeError, lambda: convolution(b, d, prime=p, dyadic=True))
+    raises(TypeError, lambda: convolution(b, d, fft=True, dps=2, dyadic=True))
+    raises(TypeError, lambda: convolution(b, d, ntt=True, prime=p, dyadic=True))
+    raises(TypeError, lambda: convolution(b, d, fft=True, dyadic=True))
     raises(TypeError, lambda: convolution(a, b, dps=2, dyadic=True))
     raises(TypeError, lambda: convolution(b, c, prime=p, dyadic=True))
 
     # subset
     assert convolution(a, b, subset=True) == convolution_subset(a, b) == \
             convolution(a, b, subset=True, dyadic=False) == \
-                convolution(a, b, subset=True)
+                convolution(a, b, subset=True, fft=False) == \
+                    convolution(a, b, subset=True, fft=False, ntt=False)
     assert convolution(a, b, subset=False) == convolution(a, b)
     raises(TypeError, lambda: convolution(a, b, subset=True, dyadic=True))
+    raises(TypeError, lambda: convolution(b, c, subset=True, fft=True))
     raises(TypeError, lambda: convolution(c, d, subset=True, dps=6))
     raises(TypeError, lambda: convolution(a, c, subset=True, prime=q))
 
@@ -63,6 +72,12 @@ def test_cyclic_convolution():
 
     assert convolution([1, 2, 3], [4, 5, 6], cycle=3) == [31, 31, 28]
 
+    assert convolution(a, b, fft=True, cycle=4) == \
+            convolution(a, b, cycle=4)
+
+    assert convolution(a, b, fft=True, dps=3, cycle=4) == \
+            convolution(a, b, dps=3, cycle=4)
+
     a = [S(1)/3, S(7)/3, S(5)/9, S(2)/7, S(5)/8]
     b = [S(3)/5, S(4)/7, S(7)/8, S(8)/9]
 
@@ -280,47 +295,3 @@ def test_convolution_subset():
 
     raises(TypeError, lambda: convolution_subset(x, z))
     raises(TypeError, lambda: convolution_subset(S(7)/3, u))
-
-
-def test_covering_product():
-    assert covering_product([], []) == []
-    assert covering_product([], [S(1)/3]) == []
-    assert covering_product([6 + 3*I/7], [S(2)/3]) == [4 + 2*I/7]
-
-    a = [1, S(5)/8, sqrt(7), 4 + 9*I]
-    b = [66, 81, 95, 49, 37, 89, 17]
-    c = [3 + 2*I/3, 51 + 72*I, 7, S(7)/15, 91]
-
-    assert covering_product(a, b) == [66, 1383/8, 95 + 161*sqrt(7),
-                                        130*sqrt(7) + 1303 + 2619*I, 37,
-                                        671/4, 17 + 54*sqrt(7),
-                                        89*sqrt(7) + 4661/8 + 1287*I]
-
-    assert covering_product(b, c) == [198 + 44*I, 7740 + 10638*I,
-                                        1412 + 190*I/3, 42684/5 + 31202*I/3,
-                                        9484 + 74*I/3, 22163 + 27394*I/3,
-                                        10621 + 34*I/3, 90236/15 + 1224*I]
-
-    assert covering_product(a, c) == covering_product(c, a)
-    assert covering_product(b, c[:-1]) == [198 + 44*I, 7740 + 10638*I,
-                                         1412 + 190*I/3, 42684/5 + 31202*I/3,
-                                         111 + 74*I/3, 6693 + 27394*I/3,
-                                         429 + 34*I/3, 23351/15 + 1224*I]
-
-    assert covering_product(a, c[:-1]) == [3 + 2*I/3,
-                            339/4 + 1409*I/12, 7 + 10*sqrt(7) + 2*sqrt(7)*I/3,
-                            -403 + 772*sqrt(7)/15 + 72*sqrt(7)*I + 12658*I/15]
-
-    u, v, w, x, y, z = symbols('u v w x y z')
-
-    assert covering_product([u, v, w], [x, y]) == \
-                            [u*x, u*y + v*x + v*y, w*x, w*y]
-
-    assert covering_product([u, v, w, x], [y, z]) == \
-                            [u*y, u*z + v*y + v*z, w*y, w*z + x*y + x*z]
-
-    assert covering_product([u, v], [x, y, z]) == \
-                    covering_product([x, y, z], [u, v])
-
-    raises(TypeError, lambda: covering_product(x, z))
-    raises(TypeError, lambda: covering_product(S(7)/3, u))
diff --git a/sympy/functions/combinatorial/tests/test_comb_factorials.py b/sympy/functions/combinatorial/tests/test_comb_factorials.py
index 10ecf2dc0d..b02b126673 100644
--- a/sympy/functions/combinatorial/tests/test_comb_factorials.py
+++ b/sympy/functions/combinatorial/tests/test_comb_factorials.py
@@ -401,11 +401,6 @@ def test_binomial_Mod():
     p, q = 10**5 + 3, 10**9 + 33 # prime modulo
     r, s = 10**7 + 5, 33333333 # composite modulo
 
-    n, k, m = symbols('n k m')
-    assert (binomial(n, k) % q).subs({n: s, k: p}) == Mod(binomial(s, p), q)
-    assert (binomial(n, k) % m).subs({n: 8, k: 5, m: 13}) == 4
-    assert (binomial(9, k) % 7).subs(k, 2) == 1
-
     # Lucas Theorem
     assert Mod(binomial(156675, 4433, evaluate=False), p) == Mod(binomial(156675, 4433), p)
     assert Mod(binomial(123456, 43253, evaluate=False), p) == Mod(binomial(123456, 43253), p)
diff --git a/sympy/physics/continuum_mechanics/tests/test_beam.py b/sympy/physics/continuum_mechanics/tests/test_beam.py
index 677ef7d84d..0c0a46c3c4 100644
--- a/sympy/physics/continuum_mechanics/tests/test_beam.py
+++ b/sympy/physics/continuum_mechanics/tests/test_beam.py
@@ -2,7 +2,6 @@
 from sympy.physics.continuum_mechanics.beam import Beam
 from sympy.functions import SingularityFunction, Piecewise
 from sympy.utilities.pytest import raises
-from sympy.physics.units import meter, newton, kilo, giga, milli
 
 x = Symbol('x')
 y = Symbol('y')
@@ -197,31 +196,6 @@ def test_statically_indeterminate():
     assert p == q
 
 
-def test_beam_units():
-    E = Symbol('E')
-    I = Symbol('I')
-    R1, R2 = symbols('R1, R2')
-
-    b = Beam(8*meter, 200*giga*newton/meter**2, 400*1000000*(milli*meter)**4)
-    b.apply_load(5*kilo*newton, 2*meter, -1)
-    b.apply_load(R1, 0*meter, -1)
-    b.apply_load(R2, 8*meter, -1)
-    b.apply_load(10*kilo*newton/meter, 4*meter, 0, end=8*meter)
-    b.bc_deflection = [(0*meter, 0*meter), (8*meter, 0*meter)]
-    b.solve_for_reaction_loads(R1, R2)
-    assert b.reaction_loads == {R1: -13750*newton, R2: -31250*newton}
-
-    b = Beam(3*meter, E*newton/meter**2, I*meter**4)
-    b.apply_load(8*kilo*newton, 1*meter, -1)
-    b.apply_load(R1, 0*meter, -1)
-    b.apply_load(R2, 3*meter, -1)
-    b.apply_load(12*kilo*newton*meter, 2*meter, -2)
-    b.bc_deflection = [(0*meter, 0*meter), (3*meter, 0*meter)]
-    b.solve_for_reaction_loads(R1, R2)
-    assert b.reaction_loads == {R1: -28000*newton/3, R2: 4000*newton/3}
-    assert b.deflection().subs(x, 1*meter) == 62000*meter/(9*E*I)
-
-
 def test_composite_beam():
     E = Symbol('E')
     I = Symbol('I')
diff --git a/sympy/physics/units/tests/test_util.py b/sympy/physics/units/tests/test_util.py
index 3d961f1b0f..7394706abe 100644
--- a/sympy/physics/units/tests/test_util.py
+++ b/sympy/physics/units/tests/test_util.py
@@ -115,34 +115,20 @@ def test_convert_to_tuples_of_quantities():
 
 
 def test_eval_simplify():
-    from sympy.physics.units import cm, mm, km, m, K, Quantity, kilo, foot
+    from sympy.physics.units import cm, mm, km, m, K, Quantity, kilo
     from sympy.simplify.simplify import simplify
     from sympy.core.symbol import symbols
     from sympy.utilities.pytest import raises
     from sympy.core.function import Lambda
 
-    x, y = symbols('x y')
+    x = symbols('x')
 
     assert ((cm/mm).simplify()) == 10
     assert ((km/m).simplify()) == 1000
     assert ((km/cm).simplify()) == 100000
     assert ((10*x*K*km**2/m/cm).simplify()) == 1000000000*x*kelvin
-    assert ((cm/km/m).simplify()) == 1/(10000000*centimeter)
+    assert ((cm/km/m).simplify()) == 1/(100*kilometer)
 
     assert (3*kilo*meter).simplify() == 3000*meter
     assert (4*kilo*meter/(2*kilometer)).simplify() == 2
     assert (4*kilometer**2/(kilo*meter)**2).simplify() == 4
-
-
-def test_quantity_simplify():
-    from sympy.physics.units.util import quantity_simplify
-    from sympy.physics.units import kilo, foot
-    from sympy.core.symbol import symbols
-
-    x, y = symbols('x y')
-
-    assert quantity_simplify(x*(8*kilo*newton*meter + y)) == x*(8000*meter*newton + y)
-    assert quantity_simplify(foot*inch*(foot + inch)) == foot**2*(foot + inch)/12
-    assert quantity_simplify(foot*inch*(foot*foot + inch*(foot + inch))) == foot**2*(foot**2 + inch*(foot + inch))/12
-    assert quantity_simplify(2**(foot/inch*kilo/1000)*inch) == 4096*inch
-    assert quantity_simplify(foot**2*inch + inch**2*foot) == 13*foot**3/144
diff --git a/sympy/plotting/tests/test_plot.py b/sympy/plotting/tests/test_plot.py
index 58897db510..0c8449e9f1 100644
--- a/sympy/plotting/tests/test_plot.py
+++ b/sympy/plotting/tests/test_plot.py
@@ -2,7 +2,7 @@
                    oo, LambertW, I, meijerg, exp_polar, Max, Piecewise, And)
 from sympy.plotting import (plot, plot_parametric, plot3d_parametric_line,
                             plot3d, plot3d_parametric_surface)
-from sympy.plotting.plot import unset_show, plot_contour
+from sympy.plotting.plot import unset_show
 from sympy.utilities import lambdify as lambdify_
 from sympy.utilities.pytest import skip, raises
 from sympy.plotting.experimental_lambdify import lambdify
@@ -177,21 +177,6 @@ def plot_and_save(name):
     p.save(tmp_file('%s_parametric_surface' % name))
     p._backend.close()
 
-    # Single Contour plot.
-    p = plot_contour(sin(x)*sin(y), (x, -5, 5), (y, -5, 5))
-    p.save(tmp_file('%s_contour_plot' % name))
-    p._backend.close()
-
-    # Multiple Contour plots with same range.
-    p = plot_contour(x**2 + y**2, x**3 + y**3, (x, -5, 5), (y, -5, 5))
-    p.save(tmp_file('%s_contour_plot' % name))
-    p._backend.close()
-
-    # Multiple Contour plots with different range.
-    p = plot_contour((x**2 + y**2, (x, -5, 5), (y, -5, 5)), (x**3 + y**3, (x, -3, 3), (y, -3, 3)))
-    p.save(tmp_file('%s_contour_plot' % name))
-    p._backend.close()
-
     ###
     # Examples from the 'colors' notebook
     ###
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 3fc41d623f..bc3a7992d8 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -1952,11 +1952,6 @@ def test_issue_13849():
     assert solve((t*(sqrt(5) + sqrt(2)) - sqrt(2), t), t) == []
 
 
-def test_issue_14860():
-    from sympy.physics.units import newton, kilo
-    assert solve(8*kilo*newton + x + y, x) == [-8000*newton - y]
-
-
 def test_issue_14721():
     k, h, a, b = symbols(':4')
     assert solve([
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 92a4c39233..eb468dda80 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -38,8 +38,7 @@
     solveset_real, domain_check, solveset_complex, linear_eq_to_matrix,
     linsolve, _is_function_class_equation, invert_real, invert_complex,
     solveset, solve_decomposition, substitution, nonlinsolve, solvify,
-    _is_finite_with_finite_vars, _transolve, _is_exponential,
-    _solve_expo, _term_factors)
+    _is_finite_with_finite_vars)
 
 
 a = Symbol('a', real=True)
@@ -270,7 +269,7 @@ def test_solve_invert():
         S.Reals, FiniteSet(a/log(y/b)))
 
     # issue 4504
-    assert solveset_real(2**x - 10, x) == FiniteSet(1 + log(5)/log(2))
+    assert solveset_real(2**x - 10, x) == FiniteSet(log(10)/log(2))
 
 
 def test_errorinverses():
@@ -560,6 +559,23 @@ def test_poly_gens():
         FiniteSet(-Rational(3, 2), S.Half)
 
 
+@XFAIL
+def test_uselogcombine_1():
+    assert solveset_real(log(x - 3) + log(x + 3), x) == \
+        FiniteSet(sqrt(10))
+    assert solveset_real(log(x + 1) - log(2*x - 1), x) == FiniteSet(2)
+    assert solveset_real(log(x + 3) + log(1 + 3/x) - 3) == FiniteSet(
+        -3 + sqrt(-12 + exp(3))*exp(S(3)/2)/2 + exp(3)/2,
+        -sqrt(-12 + exp(3))*exp(S(3)/2)/2 - 3 + exp(3)/2)
+
+
+@XFAIL
+def test_uselogcombine_2():
+    eq = z - log(x) + log(y/(x*(-1 + y**2/x**2)))
+    assert solveset_real(eq, x) == \
+        FiniteSet(-sqrt(y*(y - exp(z))), sqrt(y*(y - exp(z))))
+
+
 def test_solve_abs():
     x = Symbol('x')
     n = Dummy('n')
@@ -634,6 +650,14 @@ def test_solve_only_exp_1():
     assert solveset_real(exp(x) + exp(-x) - y, x) != S.EmptySet
 
 
+@XFAIL
+def test_solve_only_exp_2():
+    assert solveset_real(exp(x/y)*exp(-z/y) - 2, y) == \
+        FiniteSet((x - z)/log(2))
+    assert solveset_real(sqrt(exp(x)) + sqrt(exp(-x)) - 4, x) == \
+        FiniteSet(2*log(-sqrt(3) + 2), 2*log(sqrt(3) + 2))
+
+
 def test_atan2():
     # The .inverse() method on atan2 works only if x.is_real is True and the
     # second argument is a real constant
@@ -725,8 +749,6 @@ def test_solveset_real_exp():
     assert solveset(Eq((-3)**(x - 3), -3**39), x, S.Reals) == FiniteSet(42)
     assert solveset(Eq(2**x, y), x, S.Reals) == Intersection(S.Reals, FiniteSet(log(y)/log(2)))
 
-    assert invert_real((-2)**(2*x) - 16, 0, x) == (x, FiniteSet(2))
-
 
 def test_solve_complex_log():
     assert solveset_complex(log(x), x) == FiniteSet(1)
@@ -820,11 +842,11 @@ def test_solve_trig_simplified():
 @XFAIL
 def test_solve_lambert():
     assert solveset_real(x*exp(x) - 1, x) == FiniteSet(LambertW(1))
-    assert solveset_real(exp(x) + x, x) == FiniteSet(-LambertW(1))
     assert solveset_real(x + 2**x, x) == \
         FiniteSet(-LambertW(log(2))/log(2))
 
     # issue 4739
+    assert solveset_real(exp(log(5)*x) - 2**x, x) == FiniteSet(0)
     ans = solveset_real(3*x + 5 + 2**(-5*x + 3), x)
     assert ans == FiniteSet(-Rational(5, 3) +
                             LambertW(-10240*2**(S(1)/3)*log(2)/3)/(5*log(2)))
@@ -845,6 +867,9 @@ def test_solve_lambert():
     assert solveset_real(3*x + log(4*x), x) == \
         FiniteSet(LambertW(Rational(3, 4))/3)
 
+    assert solveset_complex(x**z*y**z - 2, z) == \
+        FiniteSet(log(2)/(log(x) + log(y)))
+
     assert solveset_real(x**x - 2) == FiniteSet(exp(LambertW(log(2))))
 
     a = Symbol('a')
@@ -855,8 +880,14 @@ def test_solve_lambert():
     assert solveset_real((a/x + exp(x/2)).diff(x), x) == \
         FiniteSet(4*LambertW(sqrt(2)*sqrt(a)/4))
 
+    assert solveset_real(1/(1/x - y + exp(y)), x) == EmptySet()
     # coverage test
+    p = Symbol('p', positive=True)
+    w = Symbol('w')
+    assert solveset_real((1/p + 1)**(p + 1), p) == EmptySet()
     assert solveset_real(tanh(x + 3)*tanh(x - 3) - 1, x) == EmptySet()
+    assert solveset_real(2*x**w - 4*y**w, w) == \
+        solveset_real((x/y)**w - 2, w)
 
     assert solveset_real((x**2 - 2*x + 1).subs(x, log(x) + 3*x), x) == \
         FiniteSet(LambertW(3*S.Exp1)/3)
@@ -890,25 +921,15 @@ def test_solve_lambert():
 
     assert solveset_real(x**3 - 3**x, x) == \
         FiniteSet(-3/log(3)*LambertW(-log(3)/3))
+    assert solveset_real(x**2 - 2**x, x) == FiniteSet(2)
+    assert solveset_real(-x**2 + 2**x, x) == FiniteSet(2)
     assert solveset_real(3**cos(x) - cos(x)**3) == FiniteSet(
         acos(-3*LambertW(-log(3)/3)/log(3)))
 
-    assert solveset_real(x**2 - 2**x, x) == \
-        solveset_real(-x**2 + 2**x, x)
-
-    assert solveset_real(3*log(x) - x*log(3)) == FiniteSet(
-        -3*LambertW(-log(3)/3)/log(3),
-        -3*LambertW(-log(3)/3, -1)/log(3))
-
-    assert solveset_real(LambertW(2*x) - y) == FiniteSet(
-        y*exp(y)/2)
-
-
-@XFAIL
-def test_other_lambert():
-    a = S(6)/5
-    assert solveset_real(x**a - a**x, x) == FiniteSet(
-        a, -a*LambertW(-log(a)/a)/log(a))
+    assert solveset_real(4**(x/2) - 2**(x/3), x) == FiniteSet(0)
+    assert solveset_real(5**(x/2) - 2**(x/3), x) == FiniteSet(0)
+    b = sqrt(6)*sqrt(log(2))/sqrt(log(5))
+    assert solveset_real(5**(x/2) - 2**(3/x), x) == FiniteSet(-b, b)
 
 
 def test_solveset():
@@ -1124,11 +1145,6 @@ def test_linsolve():
     f = Function('f')
     assert linsolve([f(x) - 5], f(x)) == FiniteSet((5,))
 
-    # Issue #14860
-    from sympy.physics.units import meter, newton, kilo
-    Eqns = [8*kilo*newton + x + y, 28*kilo*newton*meter + 3*x*meter]
-    assert linsolve(Eqns, x, y) == {(-28000*newton/3, 4000*newton/3)}
-
 
 def test_solve_decomposition():
     x = Symbol('x')
@@ -1597,9 +1613,9 @@ def test_issue_10555():
     f = Function('f')
     g = Function('g')
     assert solveset(f(x) - pi/2, x, S.Reals) == \
-        ConditionSet(x, Eq(f(x) - pi/2, 0), S.Reals)
+        ConditionSet(x, Eq(2*f(x) - pi, 0), S.Reals)
     assert solveset(f(g(x)) - pi/2, g(x), S.Reals) == \
-        ConditionSet(g(x), Eq(f(g(x)) - pi/2, 0), S.Reals)
+        ConditionSet(g(x), Eq(2*f(g(x)) - pi, 0), S.Reals)
 
 
 def test_issue_8715():
@@ -1733,184 +1749,3 @@ def test_issue_14454():
     number = CRootOf(x**4 + x - 1, 2)
     raises(ValueError, lambda: invert_real(number, 0, x, S.Reals))
     assert invert_real(x**2, number, x, S.Reals)  # no error
-
-
-def test_term_factors():
-    assert list(_term_factors(3**x - 2)) == [-2, 3**x]
-    expr = 4**(x + 1) + 4**(x + 2) + 4**(x - 1) - 3**(x + 2) - 3**(x + 3)
-    assert set(_term_factors(expr)) == set([
-        3**(x + 2), 4**(x + 2), 3**(x + 3), 4**(x - 1), -1, 4**(x + 1)])
-
-
-#################### tests for transolve and its helpers ###############
-
-def test_transolve():
-    from sympy.abc import x
-
-    assert _transolve(3**x, x, S.Reals) == S.EmptySet
-    assert _transolve(3**x - 9**(x + 5), x, S.Reals) == FiniteSet(-10)
-
-
-# exponential tests
-def test_exponential_real():
-    from sympy.abc import x, y, z
-
-    e1 = 3**(2*x) - 2**(x + 3)
-    e2 = 4**(5 - 9*x) - 8**(2 - x)
-    e3 = 2**x + 4**x
-    e4 = exp(log(5)*x) - 2**x
-    e5 = exp(x/y)*exp(-z/y) - 2
-    e6 = 5**(x/2) - 2**(x/3)
-    e7 = 2**x + 4**x + 8**x - 84
-    e8 = 4**(x + 1) + 4**(x + 2) + 4**(x - 1) - 3**(x + 2) - 3**(x + 3)
-    e9 = -9*exp(-2*x + 5) + 4*exp(3*x + 1)
-
-    assert solveset(e1, x, S.Reals) == FiniteSet(
-        -3*log(2)/(-2*log(3) + log(2)))
-    assert solveset(e2, x, S.Reals) == FiniteSet(4/S(15))
-    assert solveset(e3, x, S.Reals) == S.EmptySet
-    assert solveset(e4, x, S.Reals) == FiniteSet(0)
-    assert solveset(e5, x, S.Reals) == Intersection(
-        S.Reals, FiniteSet(y*log(2*exp(z/y))))
-    assert solveset(e6, x, S.Reals) == FiniteSet(0)
-    assert solveset(e7, x, S.Reals) == FiniteSet(2)
-    assert solveset(e8, x, S.Reals) == FiniteSet(2)
-    assert solveset(e9, x, S.Reals) == FiniteSet(log(9*exp(4)/4)/5)
-
-    assert solveset_real(-9*exp(-2*x + 5) + 2**(x + 1), x) == FiniteSet(
-        -((-5 - 2*log(3) + log(2))/(log(2) + 2)))
-    assert solveset_real(4**(x/2) - 2**(x/3), x) == FiniteSet(0)
-    b = sqrt(6)*sqrt(log(2))/sqrt(log(5))
-    assert solveset_real(5**(x/2) - 2**(3/x), x) == FiniteSet(-b, b)
-
-    # coverage test
-    C1, C2 = symbols('C1 C2')
-    f = Function('f')
-    assert solveset_real(C1 + C2/x**2 - exp(-f(x)), f(x)) == Intersection(
-        S.Reals, FiniteSet(-log(C1 + C2/x**2)))
-    y = symbols('y', positive=True)
-    assert solveset_real(x**2 - y**2/exp(x), y) == Intersection(
-        S.Reals, FiniteSet(-sqrt(x**2*exp(x)), sqrt(x**2*exp(x))))
-    p = Symbol('p', positive=True)
-    assert solveset_real((1/p + 1)**(p + 1), p) == EmptySet()
-
-
-@XFAIL
-def test_exponential_complex():
-    from sympy.abc import x
-    from sympy import Dummy
-    n = Dummy('n')
-
-    assert solveset_complex(2**x + 4**x, x) == imageset(
-        Lambda(n, I*(2*n*pi + pi)/log(2)), S.Integers)
-    assert solveset_complex(x**z*y**z - 2, z) == FiniteSet(
-        log(2)/(log(x) + log(y)))
-    assert solveset_complex(4**(x/2) - 2**(x/3), x) == imageset(
-        Lambda(n, 3*n*I*pi/log(2)), S.Integers)
-    assert solveset(2**x + 32, x) == imageset(
-        Lambda(n, (I*(2*n*pi + pi) + 5*log(2))/log(2)), S.Integers)
-
-    eq = (2**exp(y**2/x) + 2)/(x**2 + 15)
-    a = sqrt(x)*sqrt(-log(log(2)) + log(log(2) + 2*n*I*pi))
-    assert solveset_complex(eq, y) == FiniteSet(-a, a)
-
-    union1 = imageset(Lambda(n, I*(2*n*pi - 2*pi/3)/log(2)), S.Integers)
-    union2 = imageset(Lambda(n, I*(2*n*pi + 2*pi/3)/log(2)), S.Integers)
-    assert solveset(2**x + 4**x + 8**x, x) == Union(union1, union2)
-
-    eq = 4**(x + 1) + 4**(x + 2) + 4**(x - 1) - 3**(x + 2) - 3**(x + 3)
-    res = solveset(eq, x)
-    num = 2*n*I*pi - 4*log(2) + 2*log(3)
-    den = -2*log(2) + log(3)
-    ans = imageset(Lambda(n, num/den), S.Integers)
-    assert res == ans
-
-
-def test_expo_conditionset():
-    from sympy.abc import x, y
-
-    f1 = (exp(x) + 1)**x - 2
-    f2 = (x + 2)**y*x - 3
-    f3 = 2**x - exp(x) - 3
-    f4 = log(x) - exp(x)
-
-    assert solveset(f1, x, S.Reals) == ConditionSet(
-        x, Eq((exp(x) + 1)**x - 2, 0), S.Reals)
-    assert solveset(f2, x, S.Reals) == ConditionSet(
-        x, Eq(x*(x + 2)**y - 3, 0), S.Reals)
-    assert solveset(f3, x, S.Reals) == ConditionSet(
-        x, Eq(2**x - exp(x) - 3, 0), S.Reals)
-    assert solveset(f4, x, S.Reals) == ConditionSet(
-        x, Eq(-exp(x) + log(x), 0), S.Reals)
-
-
-def test_exponential_symbols():
-    x, y, z = symbols('x y z', positive=True)
-    from sympy import simplify
-
-    assert solveset(z**x - y, x, S.Reals) == Intersection(
-        S.Reals, FiniteSet(log(y)/log(z)))
-
-    w = symbols('w')
-    f1 = 2*x**w - 4*y**w
-    f2 = (x/y)**w - 2
-    ans1 = solveset(f1, w, S.Reals)
-    ans2 = solveset(f2, w, S.Reals)
-    assert ans1 == simplify(ans2)
-
-    assert solveset(x**x, x, S.Reals) == S.EmptySet
-    assert solveset(x**y - 1, y, S.Reals) == FiniteSet(0)
-    assert solveset(exp(x/y)*exp(-z/y) - 2, y, S.Reals) == FiniteSet(
-        (x - z)/log(2)) - FiniteSet(0)
-
-
-@XFAIL
-def test_issue_10864():
-    assert solveset(x**(y*z) - x, x, S.Reals) == FiniteSet(1)
-
-
-@XFAIL
-def test_solve_only_exp_2():
-    assert solveset_real(sqrt(exp(x)) + sqrt(exp(-x)) - 4, x) == \
-        FiniteSet(2*log(-sqrt(3) + 2), 2*log(sqrt(3) + 2))
-
-
-def test_is_exponential():
-    x, y, z = symbols('x y z')
-
-    assert _is_exponential(3**x - 2, x) is True
-    assert _is_exponential(5**x - 7**(2 - x), x) is True
-    assert _is_exponential(sin(2**x) - 4*x, x) is False
-    assert _is_exponential(x**y - z, y) is True
-    assert _is_exponential(2**x + 4**x - 1, x) is True
-    assert _is_exponential(x**(y*z) - x, x) is False
-
-
-def test_solve_expo():
-    assert _solve_expo(3**(2*x) - 2**(x + 3), x) == \
-        2*x*log(3) - (x + 3)*log(2)
-    y = symbols('y')
-    assert _solve_expo(2**y + 4**y, y) == \
-        log(2**y) - log(-4**y)
-    assert _solve_expo(2**x + 3**x + 5**x, x) is None
-
-# end of exponential tests
-
-
-# logarithmic tests
-@XFAIL
-def test_logarithmic():
-    assert solveset_real(log(x - 3) + log(x + 3), x) == FiniteSet(
-        sqrt(10))
-    assert solveset_real(log(x + 1) - log(2*x - 1), x) == FiniteSet(2)
-    assert solveset_real(log(x + 3) + log(1 + 3/x) - 3) == FiniteSet(
-        -3 + sqrt(-12 + exp(3))*exp(S(3)/2)/2 + exp(3)/2,
-        -sqrt(-12 + exp(3))*exp(S(3)/2)/2 - 3 + exp(3)/2)
-
-    eq = z - log(x) + log(y/(x*(-1 + y**2/x**2)))
-    assert solveset_real(eq, x) == \
-        FiniteSet(-sqrt(y*(y - exp(z))), sqrt(y*(y - exp(z))))
-    eq = log(exp(2*x) + 1) + log(-tanh(x) + 1) - log(2)
-    assert solveset_real(eq, x) == EmptySet()
-
-# end of logarithmic tests
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index c70acb1e03..a06485aa69 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -165,7 +165,6 @@ def test_sample():
     assert sample(Z) in Z.pspace.domain.set
     sym, val = list(Z.pspace.sample().items())[0]
     assert sym == Z and val in Interval(0, oo)
-    assert density(Z)(-1) == 0
 
 
 def test_ContinuousRV():
diff --git a/sympy/utilities/tests/test_lambdify.py b/sympy/utilities/tests/test_lambdify.py
index 108e47a45b..149a7e3136 100644
--- a/sympy/utilities/tests/test_lambdify.py
+++ b/sympy/utilities/tests/test_lambdify.py
@@ -883,3 +883,18 @@ def test_lambdify_inspect():
     # Test that inspect.getsource works but don't hard-code implementation
     # details
     assert 'x**2' in inspect.getsource(f)
+
+def test_issue_14941():
+    x, y = Dummy(), Dummy()
+
+    # test dict
+    f1 = lambdify([x, y], {x: 3, y: 3}, 'sympy')
+    assert f1(2, 3) == {2: 3, 3: 3}
+
+    # test tuple
+    f2 = lambdify([x, y], (y, x), 'sympy')
+    assert f2(2, 3) == (3, 2)
+
+    # test list
+    f3 = lambdify([x, y], [y, x], 'sympy')
+    assert f3(2, 3) == [3, 2]