I need `solve()`

to return only positive real solutions to a symbolic complex equation, but it fails to find a useful solution:

```
sage: var('x, y', domain='positive')
sage: z = 1/(i*x + 1/(i*y + 1))
sage: Equation = z == 2
sage: solve(Equation, x, y)
([x == (-I*y + 1)/(2*y - 2*I)], [1])
```

I want solutions in the form `[x == foo, y == bar]`

. Splitting Equation into real and imaginary components does the trick:

```
sage: Equation = [z.real() == 2, z.imag() == 0]
sage: solve(Equation, x, y)
[[x == (-1/2), y == -1], [x == (1/2), y == 1]]
```

But why is that reformulation necessary to produce useful results? Isn't it equivalent to the original Equation?

Also, why does `solve()`

ignore my variables' `domain='positive'`

clause? Solve also ignores positive x and y assumptions:

```
sage: assume(x>0)
sage: assume(y>0)
sage: solve(Equation, x, y)
[[x == (-1/2), y == -1], [x == (1/2), y == 1]]
```

And it bugs out when positive x & y clauses are added to Equation:

```
sage: Equation.append(x>0, y>0)
sage: solve([Equation, x>0], x, y)
[[0 < x, [1/((y^2 + 1)*(1/(y^2 + 1)^2 + (x - y/(y^2 + 1))^2)) == 2, -(x - y/(y^2 + 1))/(1/(y^2 + 1)^2 + (x - y/(y^2 + 1))^2) == 0, x > 0, y > 0]]]
```