Computations with complex algebraic numbers?

edit

I am trying to perform some Newton-type computations: z-p(z)/p'(z) where p(z) is a polynomial over the field of algebraic numbers.

For example, this shows the sort of thing I'm trying to do:

p = (z^2+1)*(z^2+4)
q = p/(z-I)
pd = diff(p,z)
qd = diff(q,z)
sol = solve(qd=0,solution_dict='true')
a0 = sol[0][z]
b0 = a0-p.subs(z=a0)/pd.subs(z=a0)

However, this doesn't work because the expression z-I is not recognized as a factor of p. So I might try defining p as a polynomial over the field QQbar:

R.<z> = PolynomialRing(QQbar)
p = (z^2+1)*(z^2+4)
q = p.quo_rem(z-I)[0]
pd = p.differentiate(z)
qd = q.differentiate(z)

The trouble now is that I can only seem to solve qd over CC: that is, in numerical form. Both the commands

qd.roots()
qd.factor()

produce numerical roots; I can't use solve on qd to obtain closed-form expressions for the roots. I've also tried the first method above prefaced with z=QQbar['z'].0 but to no avail.

Is there some way of doing this entire computation over QQbar, that is, with exact complex numbers instead of numerical approximations?