Repository: sympy/sympy

Issue:
[bug] subs option to evalf isn't working for sets
Evalf does not seem to work on solveset. 

Code to replicate:
```
from sympy.stats import Uniform, cdf
from sympy import Symbol, Eq, solveset, solve

x = Symbol("x")
a = Symbol("a")
b = Symbol("b")

D = Uniform("d", a, b)
F_x = cdf(D)(x)

u = Symbol("u")
expr = solveset(Eq(F_x, u), x)
expr.evalf(subs={u:0.6, a: 0, b: 2})
>> returns {a - u*(a - b)}
```

If instead I use solve, it works:

```
expr = solve(Eq(F_x, u), x)[0]
expr.evalf(subs={u:0.6, a: 0, b: 2})
>> returns 1.20000000000000
```

Discussion:

Patch:
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index 812e9de5f8..5901d06fa2 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -1473,7 +1473,10 @@ def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None,
             result = evalf(self, prec + 4, options)
         except NotImplementedError:
             # Fall back to the ordinary evalf
-            v = self._eval_evalf(prec)
+            if hasattr(self, 'subs') and subs is not None:  # issue 20291
+                v = self.subs(subs)._eval_evalf(prec)
+            else:
+                v = self._eval_evalf(prec)
             if v is None:
                 return self
             elif not v.is_number:
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index 72b2fa4a0a..4ea90f2f0d 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -36,7 +36,7 @@
 
 
 @sympify_method_args
-class Set(Basic):
+class Set(Basic, EvalfMixin):
     """
     The base class for any kind of set.
 
@@ -648,6 +648,9 @@ def _boundary(self):
     def _measure(self):
         raise NotImplementedError("(%s)._measure" % self)
 
+    def _eval_evalf(self, prec):
+        return self.func(*[arg.evalf(n=prec_to_dps(prec)) for arg in self.args])
+
     @sympify_return([('other', 'Set')], NotImplemented)
     def __add__(self, other):
         return self.union(other)
@@ -866,7 +869,7 @@ def __bool__(self):
         return all([bool(s) for s in self.sets])
 
 
-class Interval(Set, EvalfMixin):
+class Interval(Set):
     """
     Represents a real interval as a Set.
 
@@ -1133,7 +1136,7 @@ def _eval_Eq(self, other):
             return false
 
 
-class Union(Set, LatticeOp, EvalfMixin):
+class Union(Set, LatticeOp):
     """
     Represents a union of sets as a :class:`Set`.
 
@@ -1293,15 +1296,6 @@ def as_relational(self, symbol):
     def is_iterable(self):
         return all(arg.is_iterable for arg in self.args)
 
-    def _eval_evalf(self, prec):
-        try:
-            return Union(*(set._eval_evalf(prec) for set in self.args))
-        except (TypeError, ValueError, NotImplementedError):
-            import sys
-            raise (TypeError("Not all sets are evalf-able"),
-                   None,
-                   sys.exc_info()[2])
-
     def __iter__(self):
         return roundrobin(*(iter(arg) for arg in self.args))
 
@@ -1526,7 +1520,7 @@ def as_relational(self, symbol):
         return And(*[set.as_relational(symbol) for set in self.args])
 
 
-class Complement(Set, EvalfMixin):
+class Complement(Set):
     r"""Represents the set difference or relative complement of a set with
     another set.
 
@@ -1731,7 +1725,7 @@ def _boundary(self):
         return S.EmptySet
 
 
-class FiniteSet(Set, EvalfMixin):
+class FiniteSet(Set):
     """
     Represents a finite set of discrete numbers.
 


Test patch:
diff --git a/sympy/core/tests/test_evalf.py b/sympy/core/tests/test_evalf.py
index dc3a1663ca..8e3b655990 100644
--- a/sympy/core/tests/test_evalf.py
+++ b/sympy/core/tests/test_evalf.py
@@ -579,3 +579,16 @@ def test_issue_13425():
 
 def test_issue_17421():
     assert N(acos(-I + acosh(cosh(cosh(1) + I)))) == 1.0*I
+
+
+def test_issue_20291():
+    from sympy import FiniteSet, Complement, Intersection, Reals, EmptySet
+    a = Symbol('a')
+    b = Symbol('b')
+    A = FiniteSet(a, b)
+    assert A.evalf(subs={a: 1, b: 2}) == FiniteSet(1.0, 2.0)
+    B = FiniteSet(a-b, 1)
+    assert B.evalf(subs={a: 1, b: 2}) == FiniteSet(-1.0, 1.0)
+
+    sol = Complement(Intersection(FiniteSet(-b/2 - sqrt(b**2-4*pi)/2), Reals), FiniteSet(0))
+    assert sol.evalf(subs={b: 1}) == EmptySet