Repository: sympy/sympy

Issue:
trigsimp with fu method throws an error
While the default method of `trigsimp` works, `fu` throws an error

```
from sympy import trigsimp, symbols, cos, sin, Eq
x = symbols('x')
expr = Eq(sin(x)**2+cos(x)**2,1)
trigsimp(expr, method='fu')
```

throws

```
AttributeError                            Traceback (most recent call last)
<ipython-input-67-4a087ad870b4> in <module>
      1 x = symbols('x')
      2 expr = Eq(sin(x)**2+cos(x)**2,1)
----> 3 trigsimp(expr, method='fu')

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/simplify/trigsimp.py in trigsimp(expr, **opts)
    506                    }[method]
    507 
--> 508     return trigsimpfunc(expr)
    509 
    510 

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/simplify/trigsimp.py in <lambda>(x)
    498 
    499     trigsimpfunc = {
--> 500         'fu': (lambda x: fu(x, **opts)),
    501         'matching': (lambda x: futrig(x)),
    502         'groebner': (lambda x: groebnersimp(x, **opts)),

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/simplify/fu.py in fu(rv, measure)
   1768             rv = TR2(rv)
   1769     if rv.has(sin, cos):
-> 1770         rv1 = fRL2(rv)
   1771         rv2 = TR8(TRmorrie(rv1))
   1772         rv = min([was, rv, rv1, rv2], key=measure)

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/strategies/core.py in minrule(expr)
    115     objective = kwargs.get('objective', identity)
    116     def minrule(expr):
--> 117         return min([rule(expr) for rule in rules], key=objective)
    118     return minrule

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/strategies/core.py in <listcomp>(.0)
    115     objective = kwargs.get('objective', identity)
    116     def minrule(expr):
--> 117         return min([rule(expr) for rule in rules], key=objective)
    118     return minrule

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/strategies/core.py in chain_rl(expr)
     42     def chain_rl(expr):
     43         for rule in rules:
---> 44             expr = rule(expr)
     45         return expr
     46     return chain_rl

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/strategies/core.py in minrule(expr)
    115     objective = kwargs.get('objective', identity)
    116     def minrule(expr):
--> 117         return min([rule(expr) for rule in rules], key=objective)
    118     return minrule

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/strategies/core.py in <listcomp>(.0)
    115     objective = kwargs.get('objective', identity)
    116     def minrule(expr):
--> 117         return min([rule(expr) for rule in rules], key=objective)
    118     return minrule

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/strategies/core.py in chain_rl(expr)
     42     def chain_rl(expr):
     43         for rule in rules:
---> 44             expr = rule(expr)
     45         return expr
     46     return chain_rl

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/simplify/fu.py in TR0(rv)
    225     # although it would be nice to use cancel, it doesn't work
    226     # with noncommutatives
--> 227     return rv.normal().factor().expand()
    228 
    229 

AttributeError: 'BooleanTrue' object has no attribute 'normal'
```



Discussion:

Patch:
diff --git a/sympy/core/relational.py b/sympy/core/relational.py
index 23843a9cf6..28bd38209d 100644
--- a/sympy/core/relational.py
+++ b/sympy/core/relational.py
@@ -301,6 +301,10 @@ def _eval_simplify(self, ratio, measure, rational, inverse):
         else:
             return self
 
+    def _eval_trigsimp(self, **opts):
+        from sympy.simplify import trigsimp
+        return self.func(trigsimp(self.lhs, **opts), trigsimp(self.rhs, **opts))
+
     def __nonzero__(self):
         raise TypeError("cannot determine truth value of Relational")
 
diff --git a/sympy/integrals/heurisch.py b/sympy/integrals/heurisch.py
index 78b6fdbbdd..ef190e26f6 100644
--- a/sympy/integrals/heurisch.py
+++ b/sympy/integrals/heurisch.py
@@ -564,9 +564,9 @@ def _exponent(g):
     A, B = _exponent(f), a + max(b, c)
 
     if A > 1 and B > 1:
-        monoms = tuple(itermonomials(V, A + B - 1 + degree_offset))
+        monoms = itermonomials(V, A + B - 1 + degree_offset)
     else:
-        monoms = tuple(itermonomials(V, A + B + degree_offset))
+        monoms = itermonomials(V, A + B + degree_offset)
 
     poly_coeffs = _symbols('A', len(monoms))
 
diff --git a/sympy/polys/monomials.py b/sympy/polys/monomials.py
index 9bff6c2477..b5ee66af29 100644
--- a/sympy/polys/monomials.py
+++ b/sympy/polys/monomials.py
@@ -10,30 +10,21 @@
 from sympy.polys.polyerrors import ExactQuotientFailed
 from sympy.polys.polyutils import PicklableWithSlots, dict_from_expr
 from sympy.utilities import public
-from sympy.core.compatibility import is_sequence
 
 @public
-def itermonomials(variables, max_degrees, min_degrees=None):
+def itermonomials(variables, max_degree, min_degree = 0):
     r"""
-    `max_degrees` and `min_degrees` are either both integers or both lists.
-    Unless otherwise specified, `min_degrees` is either 0 or [0,...,0].
+    Generate a set of monomials of the degree greater than or equal
+    to `min_degree` and less than or equal to `max_degree`.
 
-    A generator of all monomials `monom` is returned, such that
-    either
-    min_degree <= total_degree(monom) <= max_degree,
-    or
-    min_degrees[i] <= degree_list(monom)[i] <= max_degrees[i], for all i.
-
-    Case I:: `max_degrees` and `min_degrees` are both integers.
-    ===========================================================
     Given a set of variables `V` and a min_degree `N` and a max_degree `M`
     generate a set of monomials of degree less than or equal to `N` and greater
     than or equal to `M`. The total number of monomials in commutative
     variables is huge and is given by the following formula if `M = 0`:
 
-        .. math::
+    .. math::
 
-            \frac{(\#V + N)!}{\#V! N!}
+        \frac{(\#V + N)!}{\#V! N!}
 
     For example if we would like to generate a dense polynomial of
     a total degree `N = 50` and `M = 0`, which is the worst case, in 5
@@ -59,103 +50,44 @@ def itermonomials(variables, max_degrees, min_degrees=None):
         [1, x, y, x**2, x*y, y**2, x**3, x**2*y, x*y**2, y**3]
 
         >>> a, b = symbols('a, b', commutative=False)
-        >>> set(itermonomials([a, b, x], 2))
+        >>> itermonomials([a, b, x], 2)
         {1, a, a**2, b, b**2, x, x**2, a*b, b*a, x*a, x*b}
 
         >>> sorted(itermonomials([x, y], 2, 1), key=monomial_key('grlex', [y, x]))
         [x, y, x**2, x*y, y**2]
 
 
-    Case II:: `max_degrees` and `min_degrees` are both lists.
-    =========================================================
-    If max_degrees = [d_1, ..., d_n] and min_degrees = [e_1, ..., e_n],
-    the number of monomials generated is:
-
-        (d_1 - e_1 + 1) * ... * (d_n - e_n + 1)
-
-    Example
-    =======
-
-    Let us generate all monomials `monom` in variables `x`, and `y`
-    such that [1, 2][i] <= degree_list(monom)[i] <= [2, 4][i], i = 0, 1 ::
-
-        >>> from sympy import symbols
-        >>> from sympy.polys.monomials import itermonomials
-        >>> from sympy.polys.orderings import monomial_key
-        >>> from itertools import product
-        >>> from sympy.core import Mul
-        >>> from sympy.abc import x, y
-
-        >>> sorted(itermonomials([x, y], [2, 4], [1, 2]), reverse=True, key=monomial_key('lex', [x, y]))
-        [x**2*y**4, x**2*y**3, x**2*y**2, x*y**4, x*y**3, x*y**2]
     """
-    n = len(variables)
-    if is_sequence(max_degrees):
-        if len(max_degrees) != n:
-            raise ValueError('Argument sizes do not match')
-        if min_degrees is None:
-            min_degrees = [0]*n
-        elif not is_sequence(min_degrees):
-            raise ValueError('min_degrees is not a list')
-        else:
-            if len(min_degrees) != n:
-                raise ValueError('Argument sizes do not match')
-            if any(i < 0 for i in min_degrees):
-                raise ValueError("min_degrees can't contain negative numbers")
-        total_degree = False
-    else:
-        max_degree = max_degrees
-        if max_degree < 0:
-            raise ValueError("max_degrees can't be negative")
-        if min_degrees is None:
-            min_degree = 0
-        else:
-            if min_degrees < 0:
-                raise ValueError("min_degrees can't be negative")
-            min_degree = min_degrees
-        total_degree = True
-    if total_degree:
-        if min_degree > max_degree:
-            return
-        if not variables or max_degree == 0:
-            yield S(1)
-            return
-        # Force to list in case of passed tuple or other incompatible collection
-        variables = list(variables) + [S(1)]
-        if all(variable.is_commutative for variable in variables):
-            monomials_list_comm = []
-            for item in combinations_with_replacement(variables, max_degree):
-                powers = dict()
-                for variable in variables:
-                    powers[variable] = 0
-                for variable in item:
-                    if variable != 1:
-                        powers[variable] += 1
-                if max(powers.values()) >= min_degree:
-                    monomials_list_comm.append(Mul(*item))
-            for mon in set(monomials_list_comm):
-                yield mon
-        else:
-            monomials_list_non_comm = []
-            for item in product(variables, repeat=max_degree):
-                powers = dict()
-                for variable in variables:
-                    powers[variable] = 0
-                for variable in item:
-                    if variable != 1:
-                        powers[variable] += 1
-                if max(powers.values()) >= min_degree:
-                    monomials_list_non_comm.append(Mul(*item))
-            for mon in set(monomials_list_non_comm):
-                yield mon
+    if max_degree < 0 or min_degree > max_degree:
+        return set()
+    if not variables or max_degree == 0:
+        return {S(1)}
+    # Force to list in case of passed tuple or other incompatible collection
+    variables = list(variables) + [S(1)]
+    if all(variable.is_commutative for variable in variables):
+        monomials_list_comm = []
+        for item in combinations_with_replacement(variables, max_degree):
+            powers = dict()
+            for variable in variables:
+                powers[variable] = 0
+            for variable in item:
+                if variable != 1:
+                    powers[variable] += 1
+            if sum(powers.values()) >= min_degree:
+                monomials_list_comm.append(Mul(*item))
+        return set(monomials_list_comm)
     else:
-        if any(min_degrees[i] > max_degrees[i] for i in range(n)):
-            raise ValueError('min_degrees[i] must be <= max_degrees[i] for all i')
-        power_lists = []
-        for var, min_d, max_d in zip(variables, min_degrees, max_degrees):
-            power_lists.append([var**i for i in range(min_d, max_d + 1)])
-        for powers in product(*power_lists):
-            yield Mul(*powers)
+        monomials_list_non_comm = []
+        for item in product(variables, repeat=max_degree):
+            powers = dict()
+            for variable in variables:
+                powers[variable] = 0
+            for variable in item:
+                if variable != 1:
+                    powers[variable] += 1
+            if sum(powers.values()) >= min_degree:
+                monomials_list_non_comm.append(Mul(*item))
+        return set(monomials_list_non_comm)
 
 def monomial_count(V, N):
     r"""
@@ -179,7 +111,7 @@ def monomial_count(V, N):
     >>> monomial_count(2, 2)
     6
 
-    >>> M = list(itermonomials([x, y], 2))
+    >>> M = itermonomials([x, y], 2)
 
     >>> sorted(M, key=monomial_key('grlex', [y, x]))
     [1, x, y, x**2, x*y, y**2]
diff --git a/sympy/stats/__init__.py b/sympy/stats/__init__.py
index b25a23ff81..3f20f76c7d 100644
--- a/sympy/stats/__init__.py
+++ b/sympy/stats/__init__.py
@@ -43,8 +43,7 @@
 from .rv_interface import (
     cdf, characteristic_function, covariance, density, dependent, E, given, independent, P, pspace,
     random_symbols, sample, sample_iter, skewness, std, variance, where,
-    correlation, moment, cmoment, smoment, sampling_density, moment_generating_function, entropy, H,
-    quantile
+    correlation, moment, cmoment, smoment, sampling_density, moment_generating_function, entropy, H
 )
 __all__.extend(rv_interface.__all__)
 
diff --git a/sympy/stats/crv.py b/sympy/stats/crv.py
index 970c7106fc..2ad0240c21 100644
--- a/sympy/stats/crv.py
+++ b/sympy/stats/crv.py
@@ -10,7 +10,7 @@
 
 from __future__ import print_function, division
 
-from sympy import (Interval, Intersection, symbols, sympify, Dummy, nan,
+from sympy import (Interval, Intersection, symbols, sympify, Dummy,
         Integral, And, Or, Piecewise, cacheit, integrate, oo, Lambda,
         Basic, S, exp, I, FiniteSet, Ne, Eq, Union, poly, series, factorial)
 from sympy.functions.special.delta_functions import DiracDelta
@@ -288,32 +288,6 @@ def expectation(self, expr, var, evaluate=True, **kwargs):
         else:
             return Integral(expr * self.pdf(var), (var, self.set), **kwargs)
 
-    @cacheit
-    def compute_quantile(self, **kwargs):
-        """ Compute the Quantile from the PDF
-
-        Returns a Lambda
-        """
-        x, p = symbols('x, p', real=True, finite=True, cls=Dummy)
-        left_bound = self.set.start
-
-        pdf = self.pdf(x)
-        cdf = integrate(pdf, (x, left_bound, x), **kwargs)
-
-        quantile = solveset(cdf - p, x, S.Reals)
-        return Lambda(p, Piecewise((quantile, (p >= 0) & (p <= 1) ), (nan, True)))
-
-    def _quantile(self, x):
-        return None
-
-    def quantile(self, x, **kwargs):
-        """ Cumulative density function """
-        if len(kwargs) == 0:
-            quantile = self._quantile(x)
-            if quantile is not None:
-                return quantile
-        return self.compute_quantile(**kwargs)(x)
-
 class ContinuousDistributionHandmade(SingleContinuousDistribution):
     _argnames = ('pdf',)
 
@@ -401,20 +375,6 @@ def compute_moment_generating_function(self, expr, **kwargs):
         mgf = integrate(exp(t * x) * d(x), (x, -oo, oo), **kwargs)
         return Lambda(t, mgf)
 
-    @cacheit
-    def compute_quantile(self, expr, **kwargs):
-        if not self.domain.set.is_Interval:
-            raise ValueError(
-                "Quantile not well defined on multivariate expressions")
-
-        d = self.compute_cdf(expr, **kwargs)
-        x = symbols('x', real=True, finite=True, cls=Dummy)
-        p = symbols('x', real=True, positive=True, finite=True, cls=Dummy)
-
-        quantile = solveset(d(x) - p, x, self.set)
-
-        return Lambda(p, quantile)
-
     def probability(self, condition, **kwargs):
         z = Dummy('z', real=True, finite=True)
         cond_inv = False
@@ -554,14 +514,6 @@ def compute_density(self, expr, **kwargs):
         fy = sum(fx(g) * abs(g.diff(y)) for g in gs)
         return Lambda(y, fy)
 
-    def compute_quantile(self, expr, **kwargs):
-
-        if expr == self.value:
-            p = symbols("p", real=True, cls=Dummy)
-            return Lambda(p, self.distribution.quantile(p, **kwargs))
-        else:
-            return ContinuousPSpace.compute_quantile(self, expr, **kwargs)
-
 def _reduce_inequalities(conditions, var, **kwargs):
     try:
         return reduce_rational_inequalities(conditions, var, **kwargs)
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index 3ff50ac501..6791c91413 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -47,11 +47,11 @@
 
 from sympy import (log, sqrt, pi, S, Dummy, Interval, sympify, gamma,
                    Piecewise, And, Eq, binomial, factorial, Sum, floor, Abs,
-                   Lambda, Basic, lowergamma, erf, erfi,  erfinv, I, hyper,
-                   uppergamma, sinh, atan, Ne, expint)
+                   Lambda, Basic, lowergamma, erf, erfi, I, hyper, uppergamma,
+                   sinh, atan, Ne, expint)
 
 from sympy import beta as beta_fn
-from sympy import cos, sin, tan, atan, exp, besseli, besselj, besselk
+from sympy import cos, sin, exp, besseli, besselj, besselk
 from sympy.external import import_module
 from sympy.matrices import MatrixBase
 from sympy.stats.crv import (SingleContinuousPSpace, SingleContinuousDistribution,
@@ -491,9 +491,6 @@ def _moment_generating_function(self, t):
         raise NotImplementedError("The moment generating function for the "
                                   "Cauchy distribution does not exist.")
 
-    def _quantile(self, p):
-        return self.x0 + self.gamma*tan(pi*(p - S.Half))
-
 def Cauchy(name, x0, gamma):
     r"""
     Create a continuous random variable with a Cauchy distribution.
@@ -949,9 +946,6 @@ def _moment_generating_function(self, t):
         rate = self.rate
         return rate / (rate - t)
 
-    def _quantile(self, p):
-        return -log(1-p)/self.rate
-
 def Exponential(name, rate):
     r"""
     Create a continuous random variable with an Exponential distribution.
@@ -977,12 +971,12 @@ def Exponential(name, rate):
     ========
 
     >>> from sympy.stats import Exponential, density, cdf, E
-    >>> from sympy.stats import variance, std, skewness, quantile
-    >>> from sympy import Symbol, symbols
+    >>> from sympy.stats import variance, std, skewness
+    >>> from sympy import Symbol
 
     >>> l = Symbol("lambda", positive=True)
     >>> z = Symbol("z")
-    >>> p = Symbol("p")
+
     >>> X = Exponential("x", l)
 
     >>> density(X)(z)
@@ -991,9 +985,6 @@ def Exponential(name, rate):
     >>> cdf(X)(z)
     Piecewise((1 - exp(-lambda*z), z >= 0), (0, True))
 
-    >>> quantile(X)(p)
-    -log(1 - p)/lambda
-
     >>> E(X)
     1/lambda
 
@@ -1773,9 +1764,6 @@ def _characteristic_function(self, t):
     def _moment_generating_function(self, t):
         return exp(self.mu*t) * Beta(1 - self.s*t, 1 + self.s*t)
 
-    def _quantile(self, p):
-        return self.mu - self.s*log(-S.One + S.One/p)
-
 def Logistic(name, mu, s):
     r"""
     Create a continuous random variable with a logistic distribution.
@@ -2100,10 +2088,6 @@ def _moment_generating_function(self, t):
         mean, std = self.mean, self.std
         return exp(mean*t + std**2*t**2/2)
 
-    def _quantile(self, p):
-        mean, std = self.mean, self.std
-        return mean + std*sqrt(2)*erfinv(2*p - 1)
-
 def Normal(name, mean, std):
     r"""
     Create a continuous random variable with a Normal distribution.
@@ -2128,14 +2112,13 @@ def Normal(name, mean, std):
     Examples
     ========
 
-    >>> from sympy.stats import Normal, density, E, std, cdf, skewness, quantile
+    >>> from sympy.stats import Normal, density, E, std, cdf, skewness
     >>> from sympy import Symbol, simplify, pprint, factor, together, factor_terms
 
     >>> mu = Symbol("mu")
     >>> sigma = Symbol("sigma", positive=True)
     >>> z = Symbol("z")
     >>> y = Symbol("y")
-    >>> p = Symbol("p")
     >>> X = Normal("x", mu, sigma)
 
     >>> density(X)(z)
@@ -2150,9 +2133,6 @@ def Normal(name, mean, std):
     -------------------- + -
              2             2
 
-    >>> quantile(X)(p)
-    mu + sqrt(2)*sigma*erfinv(2*p - 1)
-
     >>> simplify(skewness(X))
     0
 
diff --git a/sympy/stats/drv.py b/sympy/stats/drv.py
index a7fa7ff5c3..1f7870420d 100644
--- a/sympy/stats/drv.py
+++ b/sympy/stats/drv.py
@@ -138,31 +138,6 @@ def moment_generating_function(self, t, **kwargs):
                 return mgf
         return self.compute_moment_generating_function(**kwargs)(t)
 
-    @cacheit
-    def compute_quantile(self, **kwargs):
-        """ Compute the Quantile from the PDF
-
-        Returns a Lambda
-        """
-        x = symbols('x', integer=True, finite=True, cls=Dummy)
-        p = symbols('p', real=True, finite=True, cls=Dummy)
-        left_bound = self.set.inf
-        pdf = self.pdf(x)
-        cdf = summation(pdf, (x, left_bound, x), **kwargs)
-        set = ((x, p <= cdf), )
-        return Lambda(p, Piecewise(*set))
-
-    def _quantile(self, x):
-        return None
-
-    def quantile(self, x, **kwargs):
-        """ Cumulative density function """
-        if not kwargs:
-            quantile = self._quantile(x)
-            if quantile is not None:
-                return quantile
-        return self.compute_quantile(**kwargs)(x)
-
     def expectation(self, expr, var, evaluate=True, **kwargs):
         """ Expectation of expression over distribution """
         # TODO: support discrete sets with non integer stepsizes
@@ -363,10 +338,3 @@ def compute_moment_generating_function(self, expr, **kwargs):
             return Lambda(t, self.distribution.moment_generating_function(t, **kwargs))
         else:
             raise NotImplementedError()
-
-    def compute_quantile(self, expr, **kwargs):
-        if expr == self.value:
-            p = symbols("p", real=True, finite=True, cls=Dummy)
-            return Lambda(p, self.distribution.quantile(p, **kwargs))
-        else:
-            raise NotImplementedError()
diff --git a/sympy/stats/drv_types.py b/sympy/stats/drv_types.py
index 948c382150..eb954ea325 100644
--- a/sympy/stats/drv_types.py
+++ b/sympy/stats/drv_types.py
@@ -1,4 +1,5 @@
 """
+Discrete Random Variables - Prebuilt variable types
 
 Contains
 ========
@@ -6,29 +7,21 @@
 Poisson
 Logarithmic
 NegativeBinomial
-Poisson
 YuleSimon
 Zeta
 """
 
-
 from __future__ import print_function, division
 
-from sympy import (factorial, exp, S, sympify, And, I, zeta, polylog, log, beta,
-    hyper, binomial, Piecewise, floor)
+from sympy import (factorial, exp, S, sympify, And, I, zeta, polylog, log, beta, hyper, binomial,
+                   Piecewise, floor)
 from sympy.stats import density
 from sympy.stats.drv import SingleDiscreteDistribution, SingleDiscretePSpace
 from sympy.stats.joint_rv import JointPSpace, CompoundDistribution
 from sympy.stats.rv import _value_check, RandomSymbol
 import random
 
-__all__ = ['Geometric',
-'Logarithmic',
-'NegativeBinomial',
-'Poisson',
-'YuleSimon',
-'Zeta'
-]
+__all__ = ['Geometric', 'Logarithmic', 'NegativeBinomial', 'Poisson', 'YuleSimon', 'Zeta']
 
 
 def rv(symbol, cls, *args):
diff --git a/sympy/stats/frv.py b/sympy/stats/frv.py
index 52154ada0a..740c91c5bb 100644
--- a/sympy/stats/frv.py
+++ b/sympy/stats/frv.py
@@ -11,8 +11,8 @@
 
 from itertools import product
 
-from sympy import (Basic, Symbol, symbols, cacheit, sympify, Mul,
-        And, Or, Tuple, Piecewise, Eq, Lambda, exp, I, Dummy, nan)
+from sympy import (Basic, Symbol, cacheit, sympify, Mul,
+        And, Or, Tuple, Piecewise, Eq, Lambda, exp, I, Dummy)
 from sympy.sets.sets import FiniteSet
 from sympy.stats.rv import (RandomDomain, ProductDomain, ConditionalDomain,
         PSpace, IndependentProductPSpace, SinglePSpace, random_symbols,
@@ -308,14 +308,6 @@ def compute_expectation(self, expr, rvs=None, **kwargs):
         return sum([expr.xreplace(dict(elem)) * self.prob_of(elem)
                 for elem in self.domain])
 
-    def compute_quantile(self, expr):
-        cdf = self.compute_cdf(expr)
-        p = symbols('p', real=True, finite=True, cls=Dummy)
-        set = ((nan, (p < 0) | (p > 1)),)
-        for key, value in cdf.items():
-            set = set + ((key, p <= value), )
-        return Lambda(p, Piecewise(*set))
-
     def probability(self, condition):
         cond_symbols = frozenset(rs.symbol for rs in random_symbols(condition))
         assert cond_symbols.issubset(self.symbols)
diff --git a/sympy/stats/frv_types.py b/sympy/stats/frv_types.py
index de624bfff7..f586ecd6d8 100644
--- a/sympy/stats/frv_types.py
+++ b/sympy/stats/frv_types.py
@@ -22,15 +22,8 @@
 from sympy.core.logic import fuzzy_not, fuzzy_and
 from sympy.stats.frv import (SingleFinitePSpace, SingleFiniteDistribution)
 
-__all__ = ['FiniteRV',
-'DiscreteUniform',
-'Die',
-'Bernoulli',
-'Coin',
-'Binomial',
-'Hypergeometric',
-'Rademacher'
-]
+__all__ = ['FiniteRV', 'DiscreteUniform', 'Die', 'Bernoulli', 'Coin',
+        'Binomial', 'Hypergeometric']
 
 def rv(name, cls, *args):
     density = cls(*args)
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 3b44457e93..413dbb04d2 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -3,7 +3,7 @@
 
 Defines abstract random variable type.
 Contains interfaces for probability space object (PSpace) as well as standard
-operators, P, E, sample, density, where, quantile
+operators, P, E, sample, density, where
 
 See Also
 ========
@@ -1009,51 +1009,6 @@ def sample_iter(expr, condition=None, numsamples=S.Infinity, **kwargs):
     except TypeError:
         return sample_iter_subs(expr, condition, numsamples, **kwargs)
 
-def quantile(expr, evaluate=True, **kwargs):
-    r"""
-    Return the :math:`p^{th}` order quantile of a probability distribution.
-
-    Quantile is defined as the value at which the probability of the random
-    variable is less than or equal to the given probability.
-
-    ..math::
-        Q(p) = inf{x \in (-\infty, \infty) such that p <= F(x)}
-
-    Examples
-    ========
-
-    >>> from sympy.stats import quantile, Die, Exponential
-    >>> from sympy import Symbol, pprint
-    >>> p = Symbol("p")
-
-    >>> l = Symbol("lambda", positive=True)
-    >>> X = Exponential("x", l)
-    >>> quantile(X)(p)
-    -log(1 - p)/lambda
-
-    >>> D = Die("d", 6)
-    >>> pprint(quantile(D)(p), use_unicode=False)
-    /nan  for Or(p > 1, p < 0)
-    |
-    | 1       for p <= 1/6
-    |
-    | 2       for p <= 1/3
-    |
-    < 3       for p <= 1/2
-    |
-    | 4       for p <= 2/3
-    |
-    | 5       for p <= 5/6
-    |
-    \ 6        for p <= 1
-
-    """
-    result = pspace(expr).compute_quantile(expr, **kwargs)
-
-    if evaluate and hasattr(result, 'doit'):
-        return result.doit()
-    else:
-        return result
 
 def sample_iter_lambdify(expr, condition=None, numsamples=S.Infinity, **kwargs):
     """
diff --git a/sympy/stats/rv_interface.py b/sympy/stats/rv_interface.py
index f8da368ce0..84e0525876 100644
--- a/sympy/stats/rv_interface.py
+++ b/sympy/stats/rv_interface.py
@@ -2,14 +2,13 @@
 
 from .rv import (probability, expectation, density, where, given, pspace, cdf,
         characteristic_function, sample, sample_iter, random_symbols, independent, dependent,
-        sampling_density, moment_generating_function, _value_check, quantile)
-from sympy import Piecewise, sqrt, solveset, Symbol, S, log, Eq, Lambda, exp
-from sympy.solvers.inequalities import reduce_inequalities
+        sampling_density, moment_generating_function)
+from sympy import sqrt, log, Piecewise, Symbol, Eq, Lambda, exp
 
 __all__ = ['P', 'E', 'H', 'density', 'where', 'given', 'sample', 'cdf', 'characteristic_function', 'pspace',
         'sample_iter', 'variance', 'std', 'skewness', 'covariance',
         'dependent', 'independent', 'random_symbols', 'correlation',
-        'moment', 'cmoment', 'sampling_density', 'moment_generating_function', 'quantile']
+        'moment', 'cmoment', 'sampling_density', 'moment_generating_function']
 
 
 
@@ -251,8 +250,6 @@ def skewness(X, condition=None, **kwargs):
     return smoment(X, 3, condition, **kwargs)
 
 
-
-
 P = probability
 E = expectation
 H = entropy


Test patch:
diff --git a/sympy/polys/tests/test_monomials.py b/sympy/polys/tests/test_monomials.py
index 0ce8b5913a..8b7d717ead 100644
--- a/sympy/polys/tests/test_monomials.py
+++ b/sympy/polys/tests/test_monomials.py
@@ -17,150 +17,56 @@
 
 
 def test_monomials():
+    assert itermonomials([], -1) == set()
+    assert itermonomials([], 0) == {S(1)}
+    assert itermonomials([], 1) == {S(1)}
+    assert itermonomials([], 2) == {S(1)}
+    assert itermonomials([], 3) == {S(1)}
 
-    # total_degree tests
-    assert set(itermonomials([], 0)) == {S(1)}
-    assert set(itermonomials([], 1)) == {S(1)}
-    assert set(itermonomials([], 2)) == {S(1)}
+    assert itermonomials([x], -1) == set()
+    assert itermonomials([x], 0) == {S(1)}
+    assert itermonomials([x], 1) == {S(1), x}
+    assert itermonomials([x], 2) == {S(1), x, x**2}
+    assert itermonomials([x], 3) == {S(1), x, x**2, x**3}
 
-    assert set(itermonomials([], 0, 0)) == {S(1)}
-    assert set(itermonomials([], 1, 0)) == {S(1)}
-    assert set(itermonomials([], 2, 0)) == {S(1)}
+    assert itermonomials([x, y], 0) == {S(1)}
+    assert itermonomials([x, y], 1) == {S(1), x, y}
+    assert itermonomials([x, y], 2) == {S(1), x, y, x**2, y**2, x*y}
+    assert itermonomials([x, y], 3) == \
+        {S(1), x, y, x**2, x**3, y**2, y**3, x*y, x*y**2, y*x**2}
 
-    raises(StopIteration, lambda: next(itermonomials([], 0, 1)))
-    raises(StopIteration, lambda: next(itermonomials([], 0, 2)))
-    raises(StopIteration, lambda: next(itermonomials([], 0, 3)))
-
-    assert set(itermonomials([], 0, 1)) == set()
-    assert set(itermonomials([], 0, 2)) == set()
-    assert set(itermonomials([], 0, 3)) == set()
-
-    raises(ValueError, lambda: set(itermonomials([], -1)))
-    raises(ValueError, lambda: set(itermonomials([x], -1)))
-    raises(ValueError, lambda: set(itermonomials([x, y], -1)))
-
-    assert set(itermonomials([x], 0)) == {S(1)}
-    assert set(itermonomials([x], 1)) == {S(1), x}
-    assert set(itermonomials([x], 2)) == {S(1), x, x**2}
-    assert set(itermonomials([x], 3)) == {S(1), x, x**2, x**3}
-
-    assert set(itermonomials([x, y], 0)) == {S(1)}
-    assert set(itermonomials([x, y], 1)) == {S(1), x, y}
-    assert set(itermonomials([x, y], 2)) == {S(1), x, y, x**2, y**2, x*y}
-    assert set(itermonomials([x, y], 3)) == \
-            {S(1), x, y, x**2, x**3, y**2, y**3, x*y, x*y**2, y*x**2}
 
     i, j, k = symbols('i j k', commutative=False)
-    assert set(itermonomials([i, j, k], 0)) == {S(1)}
-    assert set(itermonomials([i, j, k], 1)) == {S(1), i, j, k}
-    assert set(itermonomials([i, j, k], 2)) == \
-           {S(1), i, j, k, i**2, j**2, k**2, i*j, i*k, j*i, j*k, k*i, k*j}
+    assert itermonomials([i, j, k], 0) == {S(1)}
+    assert itermonomials([i, j, k], 1) == {S(1), i, j, k}
+    assert itermonomials([i, j, k], 2) == \
+            {S(1), i, j, k, i**2, j**2, k**2, i*j, i*k, j*i, j*k, k*i, k*j}
 
-    assert set(itermonomials([i, j, k], 3)) == \
+    assert itermonomials([i, j, k], 3) == \
             {S(1), i, j, k, i**2, j**2, k**2, i*j, i*k, j*i, j*k, k*i, k*j,
-                    i**3, j**3, k**3,
-                    i**2 * j, i**2 * k, j * i**2, k * i**2,
-                    j**2 * i, j**2 * k, i * j**2, k * j**2,
-                    k**2 * i, k**2 * j, i * k**2, j * k**2,
-                    i*j*i, i*k*i, j*i*j, j*k*j, k*i*k, k*j*k,
-                    i*j*k, i*k*j, j*i*k, j*k*i, k*i*j, k*j*i,
-            }
-
-    assert set(itermonomials([x, i, j], 0)) == {S(1)}
-    assert set(itermonomials([x, i, j], 1)) == {S(1), x, i, j}
-    assert set(itermonomials([x, i, j], 2)) == {S(1), x, i, j, x*i, x*j, i*j, j*i, x**2, i**2, j**2}
-    assert set(itermonomials([x, i, j], 3)) == \
+                        i**3, j**3, k**3,
+                        i**2 * j, i**2 * k, j * i**2, k * i**2,
+                        j**2 * i, j**2 * k, i * j**2, k * j**2,
+                        k**2 * i, k**2 * j, i * k**2, j * k**2,
+                        i*j*i, i*k*i, j*i*j, j*k*j, k*i*k, k*j*k,
+                        i*j*k, i*k*j, j*i*k, j*k*i, k*i*j, k*j*i,
+                    }
+
+    assert itermonomials([x, i, j], 0) == {S(1)}
+    assert itermonomials([x, i, j], 1) == {S(1), x, i, j}
+    assert itermonomials([x, i, j], 2) == {S(1), x, i, j, x*i, x*j, i*j, j*i, x**2, i**2, j**2}
+    assert itermonomials([x, i, j], 3) == \
             {S(1), x, i, j, x*i, x*j, i*j, j*i, x**2, i**2, j**2,
-                            x**3, i**3, j**3,
-                            x**2 * i, x**2 * j,
-                            x * i**2, j * i**2, i**2 * j, i*j*i,
-                            x * j**2, i * j**2, j**2 * i, j*i*j,
-                            x * i * j, x * j * i
-            }
-
-    # degree_list tests
-    assert set(itermonomials([], [])) == {S(1)}
-
-    raises(ValueError, lambda: set(itermonomials([], [0])))
-    raises(ValueError, lambda: set(itermonomials([], [1])))
-    raises(ValueError, lambda: set(itermonomials([], [2])))
-
-    raises(ValueError, lambda: set(itermonomials([x], [1], [])))
-    raises(ValueError, lambda: set(itermonomials([x], [1, 2], [])))
-    raises(ValueError, lambda: set(itermonomials([x], [1, 2, 3], [])))
-
-    raises(ValueError, lambda: set(itermonomials([x], [], [1])))
-    raises(ValueError, lambda: set(itermonomials([x], [], [1, 2])))
-    raises(ValueError, lambda: set(itermonomials([x], [], [1, 2, 3])))
-
-    raises(ValueError, lambda: set(itermonomials([x, y], [1, 2], [1, 2, 3])))
-    raises(ValueError, lambda: set(itermonomials([x, y, z], [1, 2, 3], [0, 1])))
-
-    raises(ValueError, lambda: set(itermonomials([x], [1], [-1])))
-    raises(ValueError, lambda: set(itermonomials([x, y], [1, 2], [1, -1])))
-
-    raises(ValueError, lambda: set(itermonomials([], [], 1)))
-    raises(ValueError, lambda: set(itermonomials([], [], 2)))
-    raises(ValueError, lambda: set(itermonomials([], [], 3)))
-
-    raises(ValueError, lambda: set(itermonomials([x, y], [0, 1], [1, 2])))
-    raises(ValueError, lambda: set(itermonomials([x, y, z], [0, 0, 3], [0, 1, 2])))
-
-    assert set(itermonomials([x], [0])) == {S(1)}
-    assert set(itermonomials([x], [1])) == {S(1), x}
-    assert set(itermonomials([x], [2])) == {S(1), x, x**2}
-    assert set(itermonomials([x], [3])) == {S(1), x, x**2, x**3}
-
-    assert set(itermonomials([x], [3], [1])) == {x, x**3, x**2}
-    assert set(itermonomials([x], [3], [2])) == {x**3, x**2}
-
-    assert set(itermonomials([x, y], [0, 0])) == {S(1)}
-    assert set(itermonomials([x, y], [0, 1])) == {S(1), y}
-    assert set(itermonomials([x, y], [0, 2])) == {S(1), y, y**2}
-    assert set(itermonomials([x, y], [0, 2], [0, 1])) == {y, y**2}
-    assert set(itermonomials([x, y], [0, 2], [0, 2])) == {y**2}
-
-    assert set(itermonomials([x, y], [1, 0])) == {S(1), x}
-    assert set(itermonomials([x, y], [1, 1])) == {S(1), x, y, x*y}
-    assert set(itermonomials([x, y], [1, 2])) == {S(1), x, y, x*y, y**2, x*y**2}
-    assert set(itermonomials([x, y], [1, 2], [1, 1])) == {x*y, x*y**2}
-    assert set(itermonomials([x, y], [1, 2], [1, 2])) == {x*y**2}
-
-    assert set(itermonomials([x, y], [2, 0])) == {S(1), x, x**2}
-    assert set(itermonomials([x, y], [2, 1])) == {S(1), x, y, x*y, x**2, x**2*y}
-    assert set(itermonomials([x, y], [2, 2])) == \
-            {S(1), y**2, x*y**2, x, x*y, x**2, x**2*y**2, y, x**2*y}
-
-    i, j, k = symbols('i j k', commutative=False)
-    assert set(itermonomials([i, j, k], [0, 0, 0])) == {S(1)}
-    assert set(itermonomials([i, j, k], [0, 0, 1])) == {1, k}
-    assert set(itermonomials([i, j, k], [0, 1, 0])) == {1, j}
-    assert set(itermonomials([i, j, k], [1, 0, 0])) == {i, 1}
-    assert set(itermonomials([i, j, k], [0, 0, 2])) == {k**2, 1, k}
-    assert set(itermonomials([i, j, k], [0, 2, 0])) == {1, j, j**2}
-    assert set(itermonomials([i, j, k], [2, 0, 0])) == {i, 1, i**2}
-    assert set(itermonomials([i, j, k], [1, 1, 1])) == {1, k, j, j*k, i*k, i, i*j, i*j*k}
-    assert set(itermonomials([i, j, k], [2, 2, 2])) == \
-            {1, k, i**2*k**2, j*k, j**2, i, i*k, j*k**2, i*j**2*k**2,
-                    i**2*j, i**2*j**2, k**2, j**2*k, i*j**2*k,
-                    j**2*k**2, i*j, i**2*k, i**2*j**2*k, j, i**2*j*k,
-                    i*j**2, i*k**2, i*j*k, i**2*j**2*k**2, i*j*k**2, i**2, i**2*j*k**2
-            }
-
-    assert set(itermonomials([x, j, k], [0, 0, 0])) == {S(1)}
-    assert set(itermonomials([x, j, k], [0, 0, 1])) == {1, k}
-    assert set(itermonomials([x, j, k], [0, 1, 0])) == {1, j}
-    assert set(itermonomials([x, j, k], [1, 0, 0])) == {x, 1}
-    assert set(itermonomials([x, j, k], [0, 0, 2])) == {k**2, 1, k}
-    assert set(itermonomials([x, j, k], [0, 2, 0])) == {1, j, j**2}
-    assert set(itermonomials([x, j, k], [2, 0, 0])) == {x, 1, x**2}
-    assert set(itermonomials([x, j, k], [1, 1, 1])) == {1, k, j, j*k, x*k, x, x*j, x*j*k}
-    assert set(itermonomials([x, j, k], [2, 2, 2])) == \
-            {1, k, x**2*k**2, j*k, j**2, x, x*k, j*k**2, x*j**2*k**2,
-                    x**2*j, x**2*j**2, k**2, j**2*k, x*j**2*k,
-                    j**2*k**2, x*j, x**2*k, x**2*j**2*k, j, x**2*j*k,
-                    x*j**2, x*k**2, x*j*k, x**2*j**2*k**2, x*j*k**2, x**2, x**2*j*k**2
-            }
+                        x**3, i**3, j**3,
+                        x**2 * i, x**2 * j,
+                        x * i**2, j * i**2, i**2 * j, i*j*i,
+                        x * j**2, i * j**2, j**2 * i, j*i*j,
+                        x * i * j, x * j * i,
+                    }
+    assert itermonomials([x, y], 3, 0) == {1, x, x**2, x**2*y, x**3, x*y, x*y**2, y, y**2, y**3}
+    assert itermonomials([x, y], 3, 1) == {x, x**2, x**2*y, x**3, x*y, x*y**2, y, y**2, y**3}
+    assert itermonomials([x, y], 3, 2) == {x**2, x**2*y, x**3, x*y**2, y**2, y**3, x*y}
+    assert itermonomials([x, y], 3, 3) == {x**3, y**3, x**2*y, x*y**2}
 
 def test_monomial_count():
     assert monomial_count(2, 2) == 6
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index 9821e19872..1b0187d3ae 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -1,11 +1,9 @@
 from sympy import (Symbol, Abs, exp, S, N, pi, simplify, Interval, erf, erfc, Ne,
                    Eq, log, lowergamma, uppergamma, Sum, symbols, sqrt, And, gamma, beta,
-                   Piecewise, Integral, sin, cos, tan, atan, besseli, factorial, binomial,
+                   Piecewise, Integral, sin, cos, atan, besseli, factorial, binomial,
                    floor, expand_func, Rational, I, re, im, lambdify, hyper, diff, Or, Mul)
 from sympy.core.compatibility import range
 from sympy.external import import_module
-from sympy.functions.special.error_functions import erfinv
-from sympy.sets.sets import Intersection, FiniteSet
 from sympy.stats import (P, E, where, density, variance, covariance, skewness,
                          given, pspace, cdf, characteristic_function, ContinuousRV, sample,
                          Arcsin, Benini, Beta, BetaPrime, Cauchy,
@@ -17,7 +15,7 @@
                          QuadraticU, RaisedCosine, Rayleigh, ShiftedGompertz,
                          StudentT, Trapezoidal, Triangular, Uniform, UniformSum,
                          VonMises, Weibull, WignerSemicircle, correlation,
-                         moment, cmoment, smoment, quantile)
+                         moment, cmoment, smoment)
 from sympy.stats.crv_types import NormalDistribution
 from sympy.stats.joint_rv import JointPSpace
 from sympy.utilities.pytest import raises, XFAIL, slow, skip
@@ -42,7 +40,7 @@ def test_single_normal():
             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)
-    assert quantile(Y)(x) == Intersection(S.Reals, FiniteSet(sqrt(2)*sigma*(sqrt(2)*mu/(2*sigma) + erfinv(2*x - 1))))
+
     assert E(X, Eq(X, mu)) == mu
 
 
@@ -74,7 +72,6 @@ def test_ContinuousDomain():
 @slow
 def test_multiple_normal():
     X, Y = Normal('x', 0, 1), Normal('y', 0, 1)
-    p = Symbol("p", positive=True)
 
     assert E(X + Y) == 0
     assert variance(X + Y) == 2
@@ -93,7 +90,7 @@ def test_multiple_normal():
     assert smoment(X + Y, 3) == skewness(X + Y)
     assert E(X, Eq(X + Y, 0)) == 0
     assert variance(X, Eq(X + Y, 0)) == S.Half
-    assert quantile(X)(p) == sqrt(2)*erfinv(2*p - S.One)
+
 
 def test_symbolic():
     mu1, mu2 = symbols('mu1 mu2', real=True, finite=True)
@@ -250,12 +247,11 @@ def test_betaprime():
 def test_cauchy():
     x0 = Symbol("x0")
     gamma = Symbol("gamma", positive=True)
-    p = Symbol("p", positive=True)
 
     X = Cauchy('x', x0, gamma)
     assert density(X)(x) == 1/(pi*gamma*(1 + (x - x0)**2/gamma**2))
+    assert cdf(X)(x) == atan((x - x0)/gamma)/pi + S.Half
     assert diff(cdf(X)(x), x) == density(X)(x)
-    assert quantile(X)(p) == gamma*tan(pi*(p - S.Half)) + x0
 
     gamma = Symbol("gamma", positive=False)
     raises(ValueError, lambda: Cauchy('x', x0, gamma))
@@ -347,7 +343,6 @@ def test_erlang():
 def test_exponential():
     rate = Symbol('lambda', positive=True, real=True, finite=True)
     X = Exponential('x', rate)
-    p = Symbol("p", positive=True, real=True,finite=True)
 
     assert E(X) == 1/rate
     assert variance(X) == 1/rate**2
@@ -358,7 +353,6 @@ def test_exponential():
     assert P(X > 0) == S(1)
     assert P(X > 1) == exp(-rate)
     assert P(X > 10) == exp(-10*rate)
-    assert quantile(X)(p) == -log(1-p)/rate
 
     assert where(X <= 1).set == Interval(0, 1)
 
@@ -480,12 +474,10 @@ def test_laplace():
 def test_logistic():
     mu = Symbol("mu", real=True)
     s = Symbol("s", positive=True)
-    p = Symbol("p", positive=True)
 
     X = Logistic('x', mu, s)
     assert density(X)(x) == exp((-x + mu)/s)/(s*(exp((-x + mu)/s) + 1)**2)
     assert cdf(X)(x) == 1/(exp((mu - x)/s) + 1)
-    assert quantile(X)(p) == mu - s*log(-S(1) + 1/p)
 
 
 def test_lognormal():
diff --git a/sympy/stats/tests/test_finite_rv.py b/sympy/stats/tests/test_finite_rv.py
index 36d9594ab2..0e0b6fc205 100644
--- a/sympy/stats/tests/test_finite_rv.py
+++ b/sympy/stats/tests/test_finite_rv.py
@@ -1,13 +1,12 @@
-from sympy import (FiniteSet, S, Symbol, sqrt, nan,
+from sympy import (FiniteSet, S, Symbol, sqrt,
         symbols, simplify, Eq, cos, And, Tuple, Or, Dict, sympify, binomial,
-        cancel, exp, I, Piecewise)
+        cancel, exp, I)
 from sympy.core.compatibility import range
 from sympy.matrices import Matrix
 from sympy.stats import (DiscreteUniform, Die, Bernoulli, Coin, Binomial,
     Hypergeometric, Rademacher, P, E, variance, covariance, skewness, sample,
     density, where, FiniteRV, pspace, cdf,
-    correlation, moment, cmoment, smoment, characteristic_function,
-    moment_generating_function, quantile)
+    correlation, moment, cmoment, smoment, characteristic_function, moment_generating_function)
 from sympy.stats.frv_types import DieDistribution
 from sympy.utilities.pytest import raises
 
@@ -49,7 +48,7 @@ def test_discreteuniform():
 def test_dice():
     # TODO: Make iid method!
     X, Y, Z = Die('X', 6), Die('Y', 6), Die('Z', 6)
-    a, b, t, p = symbols('a b t p')
+    a, b, t = symbols('a b t')
 
     assert E(X) == 3 + S.Half
     assert variance(X) == S(35)/12
@@ -77,9 +76,6 @@ def test_dice():
     assert E(X + Y, Eq(X, Y)) == E(2*X)
     assert moment(X, 0) == 1
     assert moment(5*X, 2) == 25*moment(X, 2)
-    assert quantile(X)(p) == Piecewise((nan, (p > S.One) | (p < S(0))),\
-        (S.One, p <= S(1)/6), (S(2), p <= S(1)/3), (S(3), p <= S.Half),\
-        (S(4), p <= S(2)/3), (S(5), p <= S(5)/6), (S(6), p <= S.One))
 
     assert P(X > 3, X > 3) == S.One
     assert P(X > Y, Eq(Y, 6)) == S.Zero
@@ -162,13 +158,11 @@ def test_bernoulli():
     assert moment_generating_function(X)(t) == p * exp(a * t) + (-p + 1) * exp(b * t)
 
     X = Bernoulli('B', p, 1, 0)
-    z = Symbol("z")
 
     assert E(X) == p
     assert simplify(variance(X)) == p*(1 - p)
     assert E(a*X + b) == a*E(X) + b
     assert simplify(variance(a*X + b)) == simplify(a**2 * variance(X))
-    assert quantile(X)(z) == Piecewise((nan, (z > 1) | (z < 0)), (0, z <= 1 - p), (1, z <= 1))
 
     raises(ValueError, lambda: Bernoulli('B', 1.5))
     raises(ValueError, lambda: Bernoulli('B', -0.5))
@@ -216,17 +210,6 @@ def test_binomial_numeric():
             for k in range(n + 1):
                 assert P(Eq(X, k)) == binomial(n, k)*p**k*(1 - p)**(n - k)
 
-def test_binomial_quantile():
-    X = Binomial('X', 50, S.Half)
-    assert quantile(X)(0.95) == S(31)
-
-    X = Binomial('X', 5, S(1)/2)
-    p = Symbol("p", positive=True)
-    assert quantile(X)(p) == Piecewise((nan, p > S(1)), (S(0), p <= S(1)/32),\
-        (S(1), p <= S(3)/16), (S(2), p <= S(1)/2), (S(3), p <= S(13)/16),\
-        (S(4), p <= S(31)/32), (S(5), p <= S(1)))
-
-
 
 def test_binomial_symbolic():
     n = 2  # Because we're using for loops, can't do symbolic n
@@ -276,20 +259,16 @@ def test_rademacher():
 
 def test_FiniteRV():
     F = FiniteRV('F', {1: S.Half, 2: S.One/4, 3: S.One/4})
-    p = Symbol("p", positive=True)
 
     assert dict(density(F).items()) == {S(1): S.Half, S(2): S.One/4, S(3): S.One/4}
     assert P(F >= 2) == S.Half
-    assert quantile(F)(p) == Piecewise((nan, p > S.One), (S.One, p <= S.Half),\
-        (S(2), p <= S(3)/4),(S(3), True))
 
     assert pspace(F).domain.as_boolean() == Or(
         *[Eq(F.symbol, i) for i in [1, 2, 3]])
 
     raises(ValueError, lambda: FiniteRV('F', {1: S.Half, 2: S.Half, 3: S.Half}))
     raises(ValueError, lambda: FiniteRV('F', {1: S.Half, 2: S(-1)/2, 3: S.One}))
-    raises(ValueError, lambda: FiniteRV('F', {1: S.One, 2: S(3)/2, 3: S.Zero,\
-        4: S(-1)/2, 5: S(-3)/4, 6: S(-1)/4}))
+    raises(ValueError, lambda: FiniteRV('F', {1: S.One, 2: S(3)/2, 3: S.Zero, 4: S(-1)/2, 5: S(-3)/4, 6: S(-1)/4}))
 
 def test_density_call():
     from sympy.abc import p