What does dtype=complex do?

No idea. The original code was like this...

I guess it is supposed to be complex=True or something. It should be removed though because it doesn't do anything (ideally this would raise an error).

OK! Will do!

I think evaluate=False just needs to be used somewhere.

Problem is that the args are flipped in zoo - x**2.

I guess that the cause of the problem is evaluate=False being used somewhere...

I tried to rewrite the expressions some time ago, but couldn't figure out where to put the evaluate=False. May have another go, but indeed it would be better to only get one solution...

To recreate the output of solveset we need evaluate=False for the Add in zoo - x**2:

In [15]: sol = Union(FiniteSet(Integer(0)), Intersection(S.Reals, FiniteSet(Pow(Add(Mul(Integer(-1),
    ...:             Pow(x, Integer(2))), zoo, evaluate=False), Rational(1, 2)), Mul(Integer(-1),
    ...:             Pow(Add(Mul(Integer(-1), Pow(x, Integer(2))), zoo, evaluate=False), Rational(1, 2))))))
    ...: 
    ...: expr = -2*y*exp(-x**2 - y**2)*Abs(x)
    ...: 
    ...: print(solveset(expr, y, Reals) == sol)
True

It's the 2-arg Mul:

In [6]: e = zoo - x**2

In [7]: e.args
Out[7]: 
⎛       2⎞
⎝zoo, -x ⎠

In [8]: e2 = -1 * e

In [9]: e2
Out[9]: 
 2      
x  + zoo

In [10]: e2.args
Out[10]: 
⎛ 2     ⎞
⎝x , zoo⎠

In [11]: e2 == e2.func(*e2.args)
Out[11]: False

It can be fixed with

diff --git a/sympy/core/mul.py b/sympy/core/mul.py
index 9933eb601a..960b5baf25 100644
--- a/sympy/core/mul.py
+++ b/sympy/core/mul.py
@@ -281,7 +281,7 @@ def flatten(cls, seq):
                         ar = a*r
                         if ar:
                             bargs.insert(0, ar)
-                        bargs = [Add._from_args(bargs)]
+                        bargs = [Add(*bargs)]
                         rv = bargs, [], None
             if rv:
                 return rv

or maybe _unevaluated_Add(*args) -- I read the issue but am not seeing enough context to understand the problem. Add._from_args is much faster than Add and since you are making this change in core it makes me nervous.

I'm not sure this does fix things in a more general way. It removes zoo but loses the reason for removing zoo.

Well, the option is to check each and every inverse function that the input will not lead to zoo as you suggested. I cannot really see that keeps the reason better.

My point is that once the function class is identified it will in almost all cases call _invert_real again with an evaluated version of the inverse function. So rather than checking for every evaluation if the input may yield zoo, which is a bit of a hassle and cases will have to be added once identified, it seems more general to check if the input to (the next call of) _invert_real is zoo. The only possible way that this may not work is if there are some functions that map zoo back to reals in a "controlled" way (not just inverting again). It doesn't seem to be any of those cases in the tests so far.

I added a test when I wrote this code to illustrate the point:
assert solveset(atan(x**2 - y**2)-pi/2, y, S.Reals) == EmptySet() which currently returns Intersection({-sqrt(x**2 + zoo), sqrt(x**2 + zoo)}, Reals).

I think it would be better if invert_real only returned sets that are strictly contained in the domain of a function. There are other ways this can go wrong like if if the subset test is indeterminate then a zoo case is not removed.

I think it would be better if invert_real only returned sets that are strictly contained in the domain of a function

Since it has the domain argument, I agree. I think I fixed a problem where the domain was integers but it wasn't return only integers. Those are troublesome bugs because you assume that the routine can be counted on to respect the domain and that you won't have worry about this in the the calling routine.

I agree in general. The only other way I might think of is

g_ys = Intersection(g_ys, S.Reals)

which will also cover other cases (however, I believe that all the inverse functions return real-valued results for real-valued input, except for those that returns zoo, but I may be wrong). I'm a bit worried that this could explode in certain cases. (Intersection(ImageSet(..., Intersection(ImageSet(..., Intersection(...)), S.Reals), S.Reals)) Also, I am not sure how computationally efficient this is, but can give it a try.

Will also be interesting see if that deals with: solveset(I*pi-log(x), x, S.Reals), which, if I guess correctly, will have an imaginary intermediate value.

I tried to fix this at the intersection-level in #22045.

Isn't this always True? (Or rather is it not safe to assume that it is always True?)

It may be that you are correct.

The thing was that just adding the domain sometimes led to some bad looking results. Typically, the domain was not really required in some cases and still remained. I do not remember an exact example, but, say, Intersection(ConditionSet(... S.Reals), S.Reals) or something similar. This was one of many attempts to get rid of those that at least didn't add anything "stupid".

Maybe it would be better to just add the domain unless it is S.Complexes and use that as input for improving Intersection?

Ahh, forgot what I had already written. I'll try to find a better way of doing this.

This was the problem:
solveset(Abs(x) - n, x, S.Reals) currently returns ConditionSet(x, Contains(n, Interval(0, oo)), {-n, n}), by Intersecting with the Reals domain we get ConditionSet(x, Contains(_n, Interval(0, oo)), Intersection({-_n, _n}, Reals)).

I realize that I was trying to solve issues that I didn't fully understood. So I'll just revert this part.