For any complex value of c, the following polynomial has 6 complex root values of p:
$$1+c+2 c^2+c^3+p+2 c p+c^2 p+p^2+3 c p^2+3 c^2 p^2+p^3+2 c p^3+p^4+3 c p^4+p^5+p^6$$
For a general 6th order polynomial we cannot derive an exact formula for the roots. However in this particular case, we also happen to know that the 6 roots are in two cyclic groups where
$$p_1 = p_3^2+c$$ $$p_2 = p_1^2+c$$ $$p_3 = p_2^2+c$$
and
$$p_4 = p_6^2+c$$ $$p_5 = p_4^2+c$$ $$p_6 = p_5^2+c$$
Is there any way we can use that information to derive exact formulas for $p_1$ and $p_4$ in terms of c? If not, can we use our knowledge of the cyclicity of the roots to simplify or improve numerical solving? And is there some interesting relationship between $p_1$ and $p_4$?
Numerical calculation of roots without using the cyclicity information is not that slow for low-order polynomials, but I need a solution I can generalize to much higher order. Polynomials of arbitrarily high order of the sort I'm working with can be generated with the following Mathematica code. The first two lines just define $q_{0}\equiv p$ and $q_{n>0}\equiv q_{n-1}^{2}+c$, and the rest is simple enough to read as pseudo-code:
q[0, p_, c_] := p;
q[n_, p_, c_] := q[n - 1, p, c]^2 + c;
maxCyclicity = 5;
Do[
w[i] = q[i, p, c] - p;
Do[
If[Mod[i, j] == 0,
w[i] = PolynomialQuotient[w[i], w[j], p]],
{j, 1, i - 1}],
{i, 1, maxCyclicity}];
For any $m$ the polynomial $w_m(p)$ has $n*m$ roots in $n$ cycles of $m$ each, where $n$ is an integer. I used $m=3$ in the simple example above. Is it possible to somehow solve directly for just $n$ roots, one from each $m$-cycle? And is there some explicable relationship between those $n$ roots?