I find a solution for the example in UPD using solveset and as_set:

from sympy import *
var('x', reals=True)
range = (x > 2) | (x < 0)
equation = Eq(x ** 2 - 1, 0)
print(solveset(equation, x, range.as_set()))
# Output: {-1}