The comment by Tom gave me an idea. This is no complete answer yet, so please expand if you can think of ways to proceed. Let's start with what he wrote:
$$\exists z\in\mathbb C:\quad
\lvert z-z_1\rvert=\lvert z-z_2\rvert=\lvert z-z_3\rvert=1$$
Let's square things, and express absolute values using conjugation.
$$\exists z\in\mathbb C\,\forall k\in\{1,2,3\}:\quad
1=(z-z_k)(\bar z-\bar z_k)=z\bar z-z\bar z_k-\bar zz_k+z_k\bar z_k$$
Now we can treat $a=-\bar z$ and $b=-z$ as two distinct variables and obtain
$$ab + az_k + b\bar z_k + \lvert z_k\rvert^2 - 1 = 0$$
Three conditions for two variables, which is still too much. So let's treat $c=ab$ as a third variable.
$$az_k + b\bar z_k + c = 1 - \lvert z_k\rvert^2$$
So now we have a $3\times 3$ system of equations, and once we have solved that, we can check whether $c=ab$ holds. I guess that if it does, then $a=\bar b$ will hold as well. To verify that guess: For reasons of symmetry $c$ is real. Therefore $az_k$ and $b\bar z_k$ must have opposite imaginary parts. Furthermore, $az_k\cdot b\bar z_k=c\lvert z_k\rvert^2\in\mathbb R$ so they must have opposite argument as well. If two numbers have opposite imaginary part and opposite argument, then they are conjugate to one another.
Unfortunately solving a system of three linear equations is still harder than what I'd have hoped for. One big win is that this approach is very symmetric, and can be done symbolically as well. Unfortunately, the result still looks evil:
$$
z_1^2\bar z_1^2z_2\bar z_2 - z_1\bar z_1^2z_2^2\bar z_2 - z_1^2\bar z_1z_2\bar z_2^2 + z_1\bar z_1z_2^2\bar z_2^2 - z_1^2\bar z_1^2\bar z_2z_3 + \bar z_1^2z_2^2\bar z_2z_3 + z_1^2\bar z_1\bar z_2^2z_3 - \bar z_1z_2^2\bar z_2^2z_3 + z_1\bar z_1^2\bar z_2z_3^2 - \bar z_1^2z_2\bar z_2z_3^2 - z_1\bar z_1\bar z_2^2z_3^2 + \bar z_1z_2\bar z_2^2z_3^2 - z_1^2\bar z_1^2z_2\bar z_3 + z_1\bar z_1^2z_2^2\bar z_3 + z_1^2z_2\bar z_2^2\bar z_3 - z_1z_2^2\bar z_2^2\bar z_3 + z_1^2\bar z_1^2z_3\bar z_3 - \bar z_1^2z_2^2z_3\bar z_3 - z_1^2\bar z_2^2z_3\bar z_3 + z_2^2\bar z_2^2z_3\bar z_3 - z_1\bar z_1^2z_3^2\bar z_3 + \bar z_1^2z_2z_3^2\bar z_3 + z_1\bar z_2^2z_3^2\bar z_3 - z_2\bar z_2^2z_3^2\bar z_3 + z_1^2\bar z_1z_2\bar z_3^2 - z_1\bar z_1z_2^2\bar z_3^2 - z_1^2z_2\bar z_2\bar z_3^2 + z_1z_2^2\bar z_2\bar z_3^2 - z_1^2\bar z_1z_3\bar z_3^2 + \bar z_1z_2^2z_3\bar z_3^2 + z_1^2\bar z_2z_3\bar z_3^2 - z_2^2\bar z_2z_3\bar z_3^2 + z_1\bar z_1z_3^2\bar z_3^2 - \bar z_1z_2z_3^2\bar z_3^2 - z_1\bar z_2z_3^2\bar z_3^2 + z_2\bar z_2z_3^2\bar z_3^2 + \bar z_1^2z_2^2 - 2z_1\bar z_1z_2\bar z_2 + z_1^2\bar z_2^2 - 2\bar z_1^2z_2z_3 + 2z_1\bar z_1\bar z_2z_3 + 2\bar z_1z_2\bar z_2z_3 - 2z_1\bar z_2^2z_3 + \bar z_1^2z_3^2 - 2\bar z_1\bar z_2z_3^2 + \bar z_2^2z_3^2 + 2z_1\bar z_1z_2\bar z_3 - 2\bar z_1z_2^2\bar z_3 - 2z_1^2\bar z_2\bar z_3 + 2z_1z_2\bar z_2\bar z_3 - 2z_1\bar z_1z_3\bar z_3 + 2\bar z_1z_2z_3\bar z_3 + 2z_1\bar z_2z_3\bar z_3 - 2z_2\bar z_2z_3\bar z_3 + z_1^2\bar z_3^2 - 2z_1z_2\bar z_3^2 + z_2^2\bar z_3^2=0
$$
Still, there are some applications where this formulation is better than the naive approach of splitting real and imaginary part. Namely when you are working in some subfield of $\mathbb C$ where the imaginary part is not in that field and therefore computing it comes at a high conversion cost, e.g. in a cyclotomic field. The above can stay in that field the whole time, without conversions.