\( \DeclareMathOperator{\abs}{abs} \newcommand{\ensuremath}[1]{\mbox{$#1$}} \)

Raising a complex base to a rational exponent,
numerically.
Underlying assumption:
Such 'roots' number (q), if the exponent was (p/q),
and if (p,q) are integers with a GCD of (1).
Therefore, the best that can be accomplished is,
that an arbitrary 'first' exponent be computed, and
that, in order to find the additional roots, that
first exponent be rotated in the complex plane.
An expression is later to be raised to (1/3), that is already
familiar, and that goes as follows:

(%i1) expr : ((sqrt(5003)*%i)/6+709/54);
\[\tag{expr}\label{expr}\frac{\sqrt{5003}\%{}i}{6}+\frac{709}{54}\]

Approach:
Compute logarithm, resulting in one complex number,
Multiply by (p/q),
Compute antilogarithm.
Catch:
The logarithmS of a complex number comprise
an infinite series, differing from one solution to
the next by (2*%pi*%i).

(%i2) ComplexLog(cc, n) :=
   float(log(cabs(cc))) +
   (carg(cc) + 2 * %pi * n) * %i$

The following should be a simple example,
in which the polar angle of the argument,
in the complex plane, is 45⁰:

(%i3) ComplexLog(sqrt(0.5) * (1 + %i), 0);
\[\tag{\%{}o3}\label{o3} \frac{\%{}i\ensuremath{\pi} }{4}+2.220446049250312{{10}^{-16}}\]
(%i4) ComplexLog(sqrt(0.5) * (1 + %i), 1);
\[\tag{\%{}o4}\label{o4} \frac{9\%{}i\ensuremath{\pi} }{4}+2.220446049250312{{10}^{-16}}\]
(%i5) ComplexExp(cc) :=
   float(exp(realpart(cc)) * cos(imagpart(cc))) +
   %i * float(exp(realpart(cc)) * sin(imagpart(cc)))$
(%i6) ComplexExp(%o3);
\[\tag{\%{}o6}\label{o6} 0.7071067811865476\%{}i+0.7071067811865476\]
(%i7) ComplexExp(%o4);
\[\tag{\%{}o7}\label{o7} 0.7071067811865476\%{}i+0.7071067811865476\]
(%i8) ComplexPow(cc, r, n) :=
   ComplexExp(ComplexLog(cc, n) * r)$

In this case, the polar angle is being changed,
from 45⁰ to 30⁰, by exponentiation to 2/3:

(%i9) ComplexPow(sqrt(0.5) * (1 + %i), 2/3, 0);
\[\tag{\%{}o9}\label{o9} 0.5000000000000001\%{}i+0.8660254037844388\]
(%i10) ComplexPow(sqrt(0.5) * (1 + %i), 2/3, 1);
\[\tag{\%{}o10}\label{o10} -1.0\%{}i\]
(%i11) ComplexPow(sqrt(0.5) * (1 + %i), 2/3, 2);
\[\tag{\%{}o11}\label{o11} 0.5000000000000001\%{}i-0.8660254037844388\]

Finally, the case from a previous exercise:

(%i12) H : [0.0, 0.0, 0.0]$
(%i13) for i : 1 thru 3 do (
   H[i] : ComplexPow(expr, 1/3, i-1)
)$
(%i14) H;
\[\tag{\%{}o14}\label{o14} [0.6286416626037741\%{}i+2.526378324364056,1.873586977167749\%{}i-1.807608811874181,-2.502228639771522\%{}i-0.7187695124898743]\]

Validation:

(%i15) Roots: [0.0, 0.0, 0.0]$
(%i16) for i : 1 thru 3 do (
   Roots[i] : rectform(H[i] + 61 / (9 * H[i]) + 8/3)
)$
(%i17) Roots;
\[\tag{\%{}o17}\label{o17} [7.71942331539478-3.330669073875469{{10}^{-16}}\%{}i,-1.110223024625156{{10}^{-15}}\%{}i-0.9485509570816979,1.332267629550187{{10}^{-15}}\%{}i+1.229127641686917]\]
(%i18) ProductSeries(list) := block (
   product : 1,
   for elem in list do (
       product : product * -elem
   ),
   rectform(product)
)$
(%i19) ProductSeries(Roots);
\[\tag{\%{}o19}\label{o19} 1.990086891200246{{10}^{-14}}\%{}i+8.999999999999977\]

Again, the exact answer, corresponding to the constant term,
was +9 .

(%i20) Poly(x) := rectform(x^3-8*x^2+x+9);
\[\tag{\%{}o20}\label{o20} \operatorname{Poly}(x):=\operatorname{rectform}\left( {{x}^{3}}-8{{x}^{2}}+x+9\right) \]
(%i21) poly_values : map(Poly, Roots);
\[\tag{poly\_ values}\label{poly\_ values}[-1.873758345831649{{10}^{-14}}\%{}i-2.48689957516035{{10}^{-14}},2.48689957516035{{10}^{-14}}-2.095663872197948{{10}^{-14}}\%{}i,1.77635683940025{{10}^{-14}}-1.882997489707991{{10}^{-14}}\%{}i]\]

In general, the phenomenon was happening, that
values close to 10^-14 or 10^-16 were being generated,
that could be seen as equivalent to zero, but
that were not so, just due to arithmetic round-off errors.
The last list above, is essentially equivalent to a list
of 3 zeros, each with an imaginary and a real component.

Created with wxMaxima.