Sympy yields `TypeError: unsupported operand type(s) for *: 'interval' and 'complex'` for complex rational expression

I have this sympy script:

from sympy import *
from sympy.plotting import plot_implicit

x, y = symbols('x y', real=True)

alpha = sqrt(2) / 2
expr = 1 + ((1-alpha) * x + y*I) / (1 - alpha * (x + y*I))**2

expr = Eq(abs(expr), 1)
p1 = plot_implicit(expr)

Which is trying to solve R(z) = 1 where:

As far as I can tell the expression I am giving it is correct. Why am I getting incompatible types? Note that if I get rid of the denominator I do get a plot.

Full error:

Traceback (most recent call last):
  File "/home/makogan/Documents/University/Math521/A3/A5.py", line 49, in <module>
    p1 = plot_implicit(expr)
  File "/home/makogan/.local/lib/python3.10/site-packages/sympy/plotting/plot_implicit.py", line 430, in plot_implicit
    p.show()
  File "/home/makogan/.local/lib/python3.10/site-packages/sympy/plotting/plot.py", line 240, in show
    self._backend.show()
  File "/home/makogan/.local/lib/python3.10/site-packages/sympy/plotting/plot.py", line 1533, in show
    self.process_series()
  File "/home/makogan/.local/lib/python3.10/site-packages/sympy/plotting/plot.py", line 1530, in process_series
    self._process_series(series, ax, parent)
  File "/home/makogan/.local/lib/python3.10/site-packages/sympy/plotting/plot.py", line 1398, in _process_series
    points = s.get_raster()
  File "/home/makogan/.local/lib/python3.10/site-packages/sympy/plotting/plot_implicit.py", line 87, in get_raster
    func(xinterval, yinterval)
  File "/home/makogan/.local/lib/python3.10/site-packages/sympy/plotting/experimental_lambdify.py", line 272, in __call__
    return self.lambda_func(*args, **kwargs)
  File "<string>", line 1, in <lambda>
TypeError: unsupported operand type(s) for *: 'interval' and 'complex'

Solution

That's caused by the fact that sympy's plotting module is using an old and obsolete evaluation module. If you try this new plotting module you'll get the answer:

from sympy import *
from spb import *
var("x:z")

x, y = symbols('x y', real=True)

alpha = sqrt(2) / 2
expr = 1 + ((1 - alpha) * z) / (1 - alpha * z)**2
expr = expr.subs(z, x + I*y)

expr = Eq(abs(expr), 1)
p1 = plot_implicit(expr, x)