Excellent, Dave! Thanks for that info. With the correct formula for testing root validity, the standard chaos method works, but I was right about it needing reseeding. It also has very poor coverage by default. I was able to get fair results by boosting the runs of same and alternating choices, but less satisfactory than for the quadratic Julia.

Peter's version of MIIM using contractivity is more practical than what I read about using pixel hit counts.

I used Peter Liepa's source code to write the following short Mathematica program. It seems to be working correctly:

imax = 100; z = Table[0, {imax}]; dz = Table[1, {imax}]; roots = Table[1, {imax}];

power = 2; zc = -0.74543 + 0.11301I; x1 = -1.5; y1 = -1.5; x2 = 1.5; y2 = 1.5; n = 275; image = Table[0, {n}, {n}]; i = 2;

While[i > 1, z[[ i ]] = If[roots[[ i ]] == 1, 1, -1]Sqrt[z[[i - 1]] - zc]; dz[[ i ]] = power Abs[z[[ i ]]]^(power - 1)dz[[i - 1]]; {x, y} = Round[n {(Re[z[[ i ]]] - x1)/(x2 - x1), (Im[z[[ i ]]] - y1)/(y2 - y1)}]; If[1 <= x <= n && 1 <= y <= n, image[[y, x]] = 1]; If[i < imax && dz[[ i ]] < 1000.0, i++; roots[[ i ]] =1, While[i > 1 && roots[[ i ]] == power, roots[[ i ]] =1; i--]; roots[[ i ]]++]];

ListDensityPlot[image, Mesh -> False, Frame -> False]

I tried to modify the code to create a Glynn Julia set but unfortunately I don't think it looks right (see image below):

CPow[z_, n_, k_] := Module[{theta = n(Arg[z] + 2Pi k)}, Abs[z]^n(Cos[theta] + I Sin[theta])];

imax = 1000; z = Table[0, {imax}]; dz = Table[1, {imax}]; roots = Table[1, {imax}];

power = 1.5; zc = -0.2; nroot = 3; n = 275; image = Table[0, {n}, {n}];

{x1, y1} = {-1.5, -1.5}; {x2, y2} = {1.5, 1.5}; i = 2; count = 0;

While[i > 1 && count < 500000, z[[ i ]] = CPow[z[[ i - 1]] - zc, 1/power, roots[[ i ]] - 1]; dz[[ i ]] = Abs[power]Abs[z[[ i ]]]^(power - 1)dz[[ i - 1]]; If[i < imax && dz[[ i ]] < 1.0*^10, i++; roots[[ i ]] = 1, {x, y} = Round[n {(Re[z[[ i ]]] - x1)/(x2 - x1), (Im[z[[ i ]]] - y1)/(y2 - y1)}]; If[1 <= x <= n && 1 <= y <= n, image[[y, x]]++]; While[i > 1 && roots[[ i ]] == nroot, roots[[ i ]] = 1; i--]; roots[[ i ]]++]; count++];

ListDensityPlot[Map[If[ # == 0, -1, Log[ # ]] &, image, {2}], Mesh -> False, Frame -> False, PlotRange -> {-1, 8}]

Do you have any thoughts about what I am doing wrong here?