Repository: sympy/sympy

Issue:
solveset result is too complicated when AppliedUndef
I'm not sure why the 3rd result below is so complicated. `solve` doesn't make it so and if `f(0)` is replaced with `y` it is not complicated. 
```
>>> solveset(x**2 + y + 1, x)
{-sqrt(-y - 1), sqrt(-y - 1)}
>>> solve(x**2 + f(0) + 1, x)
[-sqrt(-f(0) - 1), sqrt(-f(0) - 1)]
>>> solveset(x**2 + f(0) + 1, x)
{-I*((-re(f(0)) - 1)**2 + im(f(0))**2)**(1/4)*sin(atan2(-im(f(0)), -re(f(0)) - 1
)/2) - ((-re(f(0)) - 1)**2 + im(f(0))**2)**(1/4)*cos(atan2(-im(f(0)), -re(f(0))
- 1)/2), I*((-re(f(0)) - 1)**2 + im(f(0))**2)**(1/4)*sin(atan2(-im(f(0)), -re(f(
0)) - 1)/2) + ((-re(f(0)) - 1)**2 + im(f(0))**2)**(1/4)*cos(atan2(-im(f(0)), -re
(f(0)) - 1)/2)}
```
Perhaps a simplification decision is being based on the result not having any free symbols (rather than seeing if it is a number): 
```
>>> f(0).free_symbols  # don't simplify based on this
set()
>>> f(0).is_number # simplify into re and imag parts if *this* is true (maybe)
False
```

Discussion:

Patch:
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 59b60b2e47..4cf16df1bf 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -12,7 +12,7 @@
 from sympy.core.sympify import sympify
 from sympy.core import S, Pow, Dummy, pi, Expr, Wild, Mul, Equality
 from sympy.core.numbers import I, Number, Rational, oo
-from sympy.core.function import (Lambda, expand_complex)
+from sympy.core.function import (Lambda, expand_complex, AppliedUndef)
 from sympy.core.relational import Eq
 from sympy.core.symbol import Symbol
 from sympy.simplify.simplify import simplify, fraction, trigsimp
@@ -478,10 +478,10 @@ def _solve_as_poly(f, symbol, domain=S.Complexes):
     if result is not None:
         if isinstance(result, FiniteSet):
             # this is to simplify solutions like -sqrt(-I) to sqrt(2)/2
-            # - sqrt(2)*I/2. We are not expanding for solution with free
-            # variables because that makes the solution more complicated. For
-            # example expand_complex(a) returns re(a) + I*im(a)
-            if all([s.free_symbols == set() and not isinstance(s, RootOf)
+            # - sqrt(2)*I/2. We are not expanding for solution with symbols
+            # or undefined functions because that makes the solution more complicated.
+            # For example, expand_complex(a) returns re(a) + I*im(a)
+            if all([s.atoms(Symbol, AppliedUndef) == set() and not isinstance(s, RootOf)
                     for s in result]):
                 s = Dummy('s')
                 result = imageset(Lambda(s, expand_complex(s)), result)


Test patch:
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index a4527e63ad..121e283700 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -893,6 +893,7 @@ def test_solve_lambert():
 
 def test_solveset():
     x = Symbol('x')
+    f = Function('f')
     raises(ValueError, lambda: solveset(x + y))
     assert solveset(x, 1) == S.EmptySet
     assert solveset(x - 1, 1) == FiniteSet(x)
@@ -915,6 +916,8 @@ def test_solveset():
     assert solveset(exp(x) - 1, x) == imageset(Lambda(n, 2*I*pi*n), S.Integers)
     assert solveset(Eq(exp(x), 1), x) == imageset(Lambda(n, 2*I*pi*n),
                                                   S.Integers)
+    # issue 13825
+    assert solveset(x**2 + f(0) + 1, x) == {-sqrt(-f(0) - 1), sqrt(-f(0) - 1)}
 
 
 def test_conditionset():