Generalized Mandelbrot & Julia Stable sets.

R.H.Beresford. 7th. November 2004.

The Mandelbrot set ... has an extraordinarily elaborate structure, but it is not of any human design." (Roger Penrose [5]).

Summary.

The procedures GenMJArray and boundary[] are set up to create and display generalized Mandelbrot and Julia plots and the boundaries of their stable sets (these are often fractal). Single plots and arrays of plots with varying parameters can be created. Many different compiled functions are defined for use with GenMJArray, mainly as iterations of the generalized Mandelbrot-Julia iteration (Subscript[G, n + 1]) = (Subscript[G, n]^p) + k (Subscript[G, n]) + C (note that (Subscript[G, 0]) is replaced by (Subscript[G, n])); the result is the number of iterations required to reach instability or a prescribed iteration limit. When G is complex and p=2  the Mandelbrot plot is given by C=0 and k=1; Julia plots are given by k=0; "Mandelbar" plots can be obtained by suitable choice of k values. Functions using multiplication in many different small algebras (from [1]) are provided; fractal boundaries have been found in some cases. Investigation of these functions is rudimentary. A case-study approach has to be employed in many cases, as the array plots can only cope with four degrees of freedom, and fractal boundaries do not give good 3D plots. 
Many examples in the literature are not stable sets but are artefacts resulting from the use of too low an iteration limit; they should be called "plots" or "pseudo-sets" rather than "sets". The Mandelbrot set boundary is known [2] to be connected, so "blobs" that appear to be disconnected are actually linked to the main area by fine lines that rarely coincide with mesh points. Some "Mandelbar" iterations are only stable for a set of branched lines - i.e. a skeleton. These "skeleton sets" may be a new discovery.

1. Introduction. Revisions.

Stable sets comprise those starting values (Subscript[A, 0]) for which iterated values of (Subscript[A, n + 1]) = F[(Subscript[A, n])] remain within a bound. Mandelbrot & Julia sets ((Subscript[M, n + 1]) = (Subscript[M, n]^2) + (Subscript[M, 0 ]) or  (Subscript[J, n + 1]) = (Subscript[J, n]^2) + C, bounded in the complex field) are the best known examples; their boundaries are fractal and their appearance depends on the "Scale of Scrutiny" - similar (but not duplicate) patterns appear when the distance between grid points is reduced. Coloured plots of the region outside the set are often created [2], [3] by assigning colours according to the number of iterations required to exceed some unstable length (usually 2); these are beautiful but mathematically irrelevant artefacts that should be described as Mandelbrot or Julia plots rather than sets.

GenMJArray.m creates, and this notebook demonstrates, many stable sets. An essential feature of any generalization is that the function F is non-linear, so that chaotic and/or cyclic behaviour can arise. Generalizations include changing the function and changing from complex to other algebras, which need not be over a field. The stability criterion can be altered to distinguish between stable chaotic behaviour and stable cyclic behaviour with different cycle lengths. The plots can be of boundary points, or of (Subscript[A, n + 1]) vs (Subscript[A, n]), or of the length of cycles. Third order (Subscript[A, n + 2]) = F[(Subscript[A, n + 1]), (Subscript[A, n])] and higher order difference equations can be added.

Systems with two degrees of freedom can be plotted as {x,y}; 3D plots are only comprehensible for non-fractal boundaries; others need to be shown as multiple plots with some variables fixed. Complex Julia sets (((x + iy))^2) + 
        v + iw have four degrees of freedom. The following procedures use Dickau's presentation [4] with plots arranged in arrays corresponding to a grid of {v,w} parameter values; two principal variables {x,y} are explored in each plot. Systems with more than four degrees of freedom can be handled as case-studies by fixing some parameters.

	In GenMJArray, an iterative function is defined by the first parameter func (function descriptions and demonstrations are in Sections 3 & 4). Then come  vi, wi, which are iteration-lists for the coordinates v, w of the array of plots (though they can be single numbers to give a single plot); be warned that the w parameter plots lower values at the top. xi, yi are iteration-lists for x & y;  p is interpreted in different ways in different functions - as the power or as a set of extra parameters or coordinates; it is the maximum permitted number of iterations. (The plots can be suppressed by negating it). A final parameter can be Hue, to give a coloured plot, HueR to give a reverse-coloured plot, HueB  or HueW to give a coloured plot with black or white low points (i.e. the points with most iterations, usually the stable set), and BW to show only the low points. Another parameter should set the number of contours to be shown by ListDensityPlot, but this is not permitted by Mathematica. Hue4 &  Hue10 approximate to this, but repeat the same few colours.

As preliminary explorations often reveal symmetries, the ranges of variables may be restricted to eliminate the repeated or reflected regions. Smaller areas can be chosen for inspection at increased resolution. Dickau [4] stores the images as full arrays of iteration lengths at full resolution, showing them by ListDensityPlot[]; this is impractically voluminous for multidimensional functions. A procedure boundary[] can be used to identify points near the set boundary, giving a much more compact list.

	High resolution and high it values greatly extend the computation time, so rapid initial exploration with low values of these parameters can be used to check that the other parameters are correctly specified.

The first two examples are stable iterations of a real set with the minimal non-linearity, (Subscript[R, n + 1]) =Abs[(Subscript[R, n])]+ k(Subscript[R, n - 1])±1. There are cyclic zones where the starting point is revisited after c iterations, and chaotic zones within which every point is eventually visited.

	Dickau's Julia set example is reproduced (with some symmetric image rows suppressed) as the default case, Ex. 3. Many of his plots are outside the stable region, and should perhaps be described as Julia plots rather than Julia sets. A number of Julia & Mandelbrot examples follow.

The first generalization is to combine the Mandelbrot and Julia sets, introducing two new parameters p & k to give (Subscript[G, n + 1]) = (Subscript[G, n]^p) + k (Subscript[G, n]) + C. This reduces to the Julia and Mandelbrot sets when p=2 & either k=0 or C=0. The factorised case (Subscript[G, n + 1]) = ((Subscript[G, n]) + a) ((Subscript[G, n]) + b) reduces to Julia when a=-b, to Mandelbrot when b=0. (The Mandelbrot set is written (Subscript[M, n + 1]) = (Subscript[M, n]^2) + (Subscript[M, n]) rather than the equivalent (Subscript[M, n + 1]) = (Subscript[M, n]^2) + (Subscript[M, 0]) to allow this generalization.) The "Mandelbar" iteration [2] includes a complex conjugate term; a generalization is developed by modifying k. Some of the "typical" Mandelbar plots are misleading; the true stable set is a "skeleton" of branching lines, and is not the bounded area found by using a low iteration limit.

The next two generalizations use Quaternions {a,i,j,k} in place of complex variables, for the Mandelbrot, Julia, and extended Mandelbrot-Julia equation. Explicit expressions for (Q^2)& (Q^3)are used as the quaternion power calculation would be unacceptably slow. {i,j,k} symmetry makes these symmetrical about the {a} axis. [6] gives 3D quaternion pictures, but does not use Mathematica.

	Several other algebras (taken from [1]) are similarly investigated. Examples show regions with fractal boundaries.
	
	The plots are not included in this notebook because of their length. Even with compiled functions, plot generation is often slow; timings (3GH WindowsXP Mathematica 5) are supplied as warning in particularly slow examples.

Revisions.

Beta version, 21/9/04.
Version 1, 3/10/04. Hue4 added, used for Ex 14,15, 17/10/04
Version 1a, 7/11/04. BW added, revised Ex 12, 12a, 12b

2. Programs, Functions, Usage. (Initialise this.)

Off[
      General::"spell", General::"spell1", Part::"partd", 
      Table::"iterb"]

    tf = TraditionalForm;

Programs

Clear[cJulia, cMandelbrot, cMandelbar, cMJfactor, cMJpk, cQM2, cQJ2, 
      cQM3, cQJ3, cQMJ2, cQMJ3, cC3MJ2, cRabs, cC2p, cC2MJ2, cC2MJ3, cC3J2, 
      cC3MJ2, cC3MJ3, cDavMJ2, cKMJ2, cKMJ3, cC3KcMJ2, cC3KcMJ3, cC3KcM, 
      cC3KcJ, cC4MJ2, cC4MJ3, cCL2MJ2, cCL2MJ3]

Clear[GenMJArray, GenMJA];

    GenMJArray[func_:  cJulia, vi_: {-1.5, 3/4, 1/4}, 
          wi_: {4/3, -1/3, -1/3}, xi_: {-1.1, 1.1, 0.06}, 
          yi_: {-1.1, 1.1, 0.06}, p_:  2. , it_:  35, 
          colf_:  Automatic  ] := 
        If[NumericQ[vi], F0 = Hold[func]; P0 = p; Vi = vi; Wi = wi; 
          X0 = xi[[1]]; Y0 = yi[[1]]; Xd = xi[[3]]; 
          Yd = yi[[3]]; 
          itlim = -it - 
              1  ; 
txy = 
            Table[func[x, y, vi, wi, p, Abs[it]], {y, yi[[1]], 
                yi[[2]], yi[[3]]}, {x, xi[[1]], xi[[2]], 
                xi[[3]]}]; 
If[it > 1, 
            ListDensityPlot[txy, Mesh -> False, Frame -> False, 
              Axes -> False, AspectRatio -> Automatic, 
              ColorFunction -> 
                colf  ]], 
          Show[GraphicsArray[
              Table[ListDensityPlot[
                  Table[func[x, y, v, w, p, Abs[it]], {y, yi[[1]], 
                      yi[[2]], yi[[3]]}, {x, xi[[1]], 
                      xi[[2]], xi[[3]]}], Mesh -> False, 
                  Frame -> False, Axes -> False, 
                  AspectRatio -> 1, DisplayFunction -> Identity, 
                  ColorFunction -> 
                    colf  ], {w, 
                  wi[[1]], wi[[2]], wi[[3]]}, {v, vi[[1]], 
                  vi[[2]], vi[[3]]}]]]];

    

      GenMJA[func_, vi_, 
          wi_, xi_, yi_, p_, it_] = 
Table[
          func[x, y, vi, wi, p, it], {y, yi[[1]], yi[[2]], 
            yi[[3]]}, {x, xi[[1]], xi[[2]], 
            xi[[3]]}];

Clear[boundary]; txy =. ; X0 =. ; Xd =. ; Y0 =. ; 
    Yd =. ;
boundary[t_:  txy, a_:  4]:=
Module[{b = Min[t], m = {Length[t[[1]]]}, p, pind, pindtemp, 
              q, teml, temp, xtest},ptb = {};pindx = {};
(*Work through each line of the list of iteration lengths; p is the stable point list*)
Do[temp = {}; pind = {}; 
                p = Position[t[[y]], b]; 
                If[p ≠ {}, 
                  If[Last[p] == m, 
                    p = Most[p]]; 
xtest = -999; 
                  teml = {}; 
                  pindtemp = {}; 
  
Do[
                    If[t[[y - 1, x = p[[q, 1]]]] ≠ b || 
                        t[[y, x + 1]] ≠ b || 
                        t[[y + 1, x]] ≠ b || 
                        t[[y, x - 1]] ≠ b || 
                        If[a > 4, 
                          t[[y + 1, x + 1]] ≠ b || 
                            t[[y - 1, x + 1]] ≠ b || 
                            t[[y + 1, x - 1]] ≠ b || 
                            t[[y - 1, x - 1]] ≠ b || 
                            If[a > 
                                8, 
txy[[y - 2, 
                                      x - 1]] ≠ 
                                  b || 
txy[[y - 2, 
                                      x]] ≠ 
                                  b || 
txy[[y - 2, 
                                      x + 1]] ≠ 
                                  b || 
txy[[y + 2, 
                                      x - 1]] ≠ 
                                  b || 
txy[[y + 2, 
                                      x]] ≠ 
                                  b || 
txy[[y + 2, 
                                      x + 1]] ≠ 
                                  b || 
txy[[y - 1, 
                                      x - 2]] ≠ 
                                  b || 
txy[[y - 1, 
                                      x + 2]] ≠ 
                                  b || 
txy[[y, 
                                      x - 2]] ≠ 
                                  b || 
txy[[y, 
                                      x + 2]] ≠ 
                                  b || 
txy[[y + 1, 
                                      x - 2]] ≠ 
                                  b || 
txy[[y + 1, 
                                      x + 2]] ≠ b]], 
                      If[x == xtest + 1,   AppendTo[
                          pindtemp[[-1]], {x, 
                            y}]; 
AppendTo[
                          teml[[-1]], {X0 + Xd*x, 
                            Y0 + Yd*y - Yd}],   AppendTo[
                          pindtemp, {{x, y}}]; 
AppendTo[
                          teml, {{X0 + Xd*x, 
                              Y0 + Yd*y - Yd}}];]; 
xtest = 
                        x], {q, Length[p]}]; 
If[
                    teml ≠ {}, AppendTo[temp, teml]; 
                    AppendTo[pind, pindtemp]]]; 
If[
                  Length[Flatten[temp]] > 0, 
                  AppendTo[ptb, temp]; 
AppendTo[pindx, 
                    pind]], {y, 2, Length[t] - 1}];
ptb = Flatten[ptb, 2];
Flatten[pindx, 3]];

BW = 
        Hue[0, 0, If[# == 0, 0, 1]] &;
HueR = Hue[1 - #] &;
HueB = 
        Hue[#, 1, If[# == 0, 0, 1]] &;
HueW = 
        Hue[#, If[# == 0, 0, 1], 1] &;
Hue4 = 
        Hue[4  #, 1, If[# == 0, 0, 1]] &;
Hue10 = Hue[10  #, 1, If[# == 0, 0, 1]] &;

Functions.

Some functions (cRabs, cJulia, cMandelbrot etc) are specific cases; other functions (distinguished by the name ending with MJ2 or MJ3) implement the Mandelbrot-Julia generalization for (Subscript[G, n + 1]) = (Subscript[G, n]^2) + k (Subscript[G, n]) + C and (Subscript[G, n + 1]) = 
        (Subscript[G, n]^3) + k (Subscript[G, n]) + C, . "Pure" Mandelbrot (etc) cases can be investigated by restricting the parameters. All the functions are compiled; this usually necessitates a 2x2 array of parameters p. Algebras with 4 elements, such as quaternions, set up (Subscript[G, 0]) = {x,  y,  v,  w} and C={(Subscript[p, 1, 1]),(Subscript[p, 1, 2]),(Subscript[p, 1, 3]),(Subscript[p, 1, 4])}; k is not an element of the algebra, but is set up as individual multipliers {(Subscript[p, 2, 1]),(Subscript[p, 2, 2]),(Subscript[p, 2, 3]),(Subscript[p, 2, 4])} so that conjugates can be used. Thus k={2,-2,-2,-2} would give the quaternion conjugate and a multiplier of 2.

\(\(\(cJulia = 
        Compile[{x, y, v, w, p, {it, _Integer}}, 
          Module[{A = x + I\ y, 
              C = v + I\ w}, \(-Length[
                FixedPointList[#1\^p + C &, A, it, 
                  SameTest -> \((Abs[#2] > 2 &)\)]]\)]];\)\(
\)
    \)\)
cMandelbrot = 
        Compile[{x, y, v, w, p, {it, _Integer}}, 
          Module[{A = x + I y}, -Length[
                FixedPointList[#(1^2) + A &, A, it, 
                  SameTest -> (Abs[#2] > 
                          2 &)]]]];

    
\(\(\(cMandelbar = 
        Compile[{x, y, v, w, p, {it, _Integer}}, 
          Module[{A = x + I\ y, 
              k = w + I\ p}, \(-Length[
                FixedPointList[#\^v + k\ Conjugate[#] &, A, it, 
                  SameTest -> \((Abs[#2] > 
                          2 &)\)]]\)]];\)\(
\)
    \)\)
cMJ2factor = 
        Compile[{x, y, v, w, {p0, _Real, 2}, {it, _Integer}}, 
          Module[{A = x + I y, C = v + I w, 
              kr = p0[[1, 1]], ki = p0[[1, 2]], ab, apb}, 
            ab = C (kr + I ki);  
            apb = C + kr + I ki; -Length[
                FixedPointList[(# +  apb) (# + ab) &, A, it, 
                  SameTest -> (Abs[#2] > 
                          2 &)]]]];

    
cMJpk = 
        Compile[{x, y, v, w, {p0, _Real, 2}, {it, _Integer}}, 
          Module[{A = x + I y, C = v + I w, 
              p = p0[[1, 1]], 
              k = p0[[1, 2]]}, -Length[
                FixedPointList[#^p + # k + C &, A, it, 
                  SameTest -> (Abs[#2] > 
                          2 &)]]]];

    
cRAbs = 
        Compile[{a0, b0, d0, e0, p, {it, _Integer}}, 
          Module[{x = a0, y = b0, a = d0, i = 0, z, w}, 
            z = {x, y}; 
While[i++; 
              w = p Abs[y] - x + a; x = y; y = Chop[w]; 
              i < it && {x, y} ≠ z]; -i/6]];

\(\(\(cQM2 = 
        Compile[{a0, b0, c0, d0, 
            k, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, 
              q = {a0, b0, c0, d0}}, 
While[\(i++\); 
              q = {k\ a + a^2 - b^2 - c^2 - d^2, k\ b + 2*a*b, k\ c + 2*a*c, 
                  k\ d + 2*a*d}; 
a = q[\([1]\)]; 
              b = q[\([2]\)]; c = q[\([3]\)]; d = q[\([4]\)]; 
              a^2 + b^2 + c^2 + d^2 < 2 && i < it]; \(-i\)]];\)\(
\)
    \)\)
cQM3 = 
        Compile[{a0, b0, c0, d0, 
            k, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = q = {k a + a*(a^2 - 3*b^2 - 3*c^2 - 2*c*d - d^2), 
                    k b + 3*a^2*b + 2*a*c*(c - d) - 
                      b*(b^2 + c^2 + d^2), 
                    k c + 3*a^2*c + 2*a*b*(-c + d) - 
                      c*(b^2 + c^2 + d^2), 
                    k d + a^2*(c + 2*d) - 
                      c*(b^2 + c^2 + d^2)}; 
a = 
                q[[1]]; b = q[[2]]; c = q[[3]]; d = q[[4]]; 
              a^2 + b^2 + c^2 + d^2 < 2 && 
                i < it]; -i]];

    
cQJ2 = 
        Compile[{a0, b0, c0, 
            d0, {p0, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p0[[1, 1]], 
              cb = p0[[1, 2]], cc = p0[[1, 3]], cd = p0[[1, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + a^2 - b^2 - c^2 - d^2, cb + 2*a*b, cc + 2*a*c, 
                  cd + 2*a*d}; 
a = q[[1]]; 
              b = q[[2]]; c = q[[3]]; d = q[[4]]; 
              a^2 + b^2 + c^2 + d^2 < 2 && 
                i < it]; -i]];

    
cQJ3 = 
        Compile[{a0, b0, c0, 
            d0, {p0, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p0[[1, 1]], 
              cb = p0[[1, 2]], cc = p0[[1, 3]], cd = p0[[1, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + a*(a^2 - 3*b^2 - 3*c^2 - 2*c*d - d^2), 
                  cb + 3*a^2*b + 2*a*c*(c - d) - b*(b^2 + c^2 + d^2), 
                  cc + 3*a^2*c + 2*a*b*(-c + d) - 
                    c*(b^2 + c^2 + d^2), 
                  cd + a^2*(c + 2*d) - 
                    c*(b^2 + c^2 + d^2)}; 
a = 
                q[[1]]; b = q[[2]]; c = q[[3]]; d = q[[4]]; 
              a^2 + b^2 + c^2 + d^2 < 2 && 
                i < it]; -i]];

    
\(\(\(cQMJ2 = 
        Compile[{a0, b0, c0, 
            d0, {p0, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p0[\([1, 1]\)], 
              cb = p0[\([1, 2]\)], cc = p0[\([1, 3]\)], cd = p0[\([1, 4]\)], 
              ka = p0[\([2, 1]\)], kb = p0[\([2, 2]\)], kc = p0[\([2, 3]\)], 
              kd = p0[\([2, 4]\)], 
              q = {a0, b0, c0, d0}}, 
While[\(i++\); 
              q = {ka\ a + ca + a^2 - b^2 - c^2 - d^2, kb\ b + cb + 2*a*b, 
                  kc\ c + cc + 2*a*c, 
                  kd\ d + cd + 2*a*d}; 
a = q[\([1]\)]; 
              b = q[\([2]\)]; c = q[\([3]\)]; d = q[\([4]\)]; 
              a^2 + b^2 + c^2 + d^2 < 2 && i < it]; \(-i\)]];\)\(
\)
    \)\)
cQMJ3 = 
        Compile[{a0, b0, c0, 
            d0, {p0, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p0[[1, 1]], 
              cb = p0[[1, 2]], cc = p0[[1, 3]], cd = p0[[1, 4]], 
              ka = p0[[2, 1]], kb = p0[[2, 2]], kc = p0[[2, 3]], 
              kd = p0[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ka a + ca + a*(a^2 - 3*b^2 - 3*c^2 - 2*c*d - d^2), 
                  kb b + cb + 3*a^2*b + 2*a*c*(c - d) - 
                    b*(b^2 + c^2 + d^2), 
                  kc c + cc + 3*a^2*c + 2*a*b*(-c + d) - 
                    c*(b^2 + c^2 + d^2), 
                  kd d + cd + a^2*(c + 2*d) - 
                    c*(b^2 + c^2 + d^2)}; 
a = 
                q[[1]]; b = q[[2]]; c = q[[3]]; d = q[[4]]; 
              a^2 + b^2 + c^2 + d^2 < 2 && i < it]; -i]];

cC2p = 
        Compile[{a0, b0, d, e, p, {it, _Integer}}, 
          Module[{a = a0, b = b0, i = 0, x, y}, 
            While[i++; x = Sign[a - b]*(Abs[a - b]^p/2) + d; 
              y = Sign[a + b]*(Abs[a + b]^p/2) + 
                  e; 
i < it && Abs[a - b] < 2 && 
                Abs[a + b] < 2, a = x; 
              b = y;]; -i]];

    
cC2MJ2 = 
        Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, i = 0, v, w, ca = p[[1, 1]], 
              cb = p[[1, 2]], ka = p[[2, 1]], kb = p[[2, 2]]}, 
            While[i++; v = (a^2) + (b^2) + ka a + ca; 
              w = 2 a b + kb b + cb; 
i < it && 
                Abs[a + b] < 2 && Abs[a - b] < 2, a = v; 
              b = w]; -i]];

    
cC2MJ3 = 
        Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, i = 0, v, w, ca = p[[1, 1]], 
              cb = p[[1, 2]], ka = p[[2, 1]], kb = p[[2, 2]]}, 
            While[i++; v = (a^3) + 3 a (b^2) + ka a + ca; 
              w = 3 (a^2) b + (b^3) + kb b + cb; 
i < it && 
                Abs[a + b] < 2 && Abs[a - b] < 2, a = v; 
              b = w]; -i]];

cC3J2=Compile[{a0, b0, c0, d, {p, _Real, 2}, {it, _Integer}},Module[{a = a0, b = b0, c = c0, e = p[[1, 1]], 
                      f = p[[1, 2]], i = 0, x, y},While[i++;x = (a^2) + 2 b c + d ;y = 2 a b + (c^2) + e;c = (b^2) + 2 a c + f;
a = x;b = y;
i < it&& (((a - b))^2)+(((b - c))^2)+(((c - a))^2)<4&&Abs[a + b + c] < 6];-i]];

  
cC3MJ2=Compile[{a0, b0, c0, k, {p, _Real, 2}, {it, _Integer}},Module[{a = a0, b = b0, c = c0, ca = p[[1, 1]], 
                      cb = p[[1, 2]], cc = p[[1, 3]], 
                      da = p[[2, 1]], db = p[[2, 2]], 
                      dc = p[[2, 3]], i = 0, x, y},While[i++;x = (a^2) + 2 b c + a k da  + ca;y = 2 a b + (c^2) + b k db + cb;c = (b^2) + 2 a c + c k dc + cc;
a = x;b = y;
i < it&& (((a - b))^2)+(((b - c))^2)+(((c - a))^2)<4&&Abs[a + b + c] < 6];-i]];

  
cC3MJ3=Compile[{a0, b0, c0, k, {p, _Real, 2}, {it, _Integer}},Module[{a = a0, b = b0, c = c0, ca = p[[1, 1]], 
                      cb = p[[1, 2]], cc = p[[1, 3]], 
                      da = p[[2, 1]], db = p[[2, 2]], 
                      dc = p[[2, 3]], i = 0, x, y},While[i++;x = 
                            (a^3) + (b^3) + 6 a b c + (c^3) + a k da  + 
                              ca;y = 
                            3 ((a^2) b + (b^2) c + a (c^2)) + b k db + 
                              cb;c = 
                            3 (a (b^2) + (a^2) c + b (c^2)) + c k dc + 
                              cc;
a = x;b = y;
i < it&& (((a - b))^2)+(((b - c))^2)+(((c - a))^2)<4&&Abs[a + b + c] < 6];-i]];

cC3MbJ2=Compile[{a0, b0, c0, k, {p, _Real, 2}, {it, _Integer}},Module[{a = a0, b = b0, c = c0, ca = p[[1, 1]], 
                    cb = p[[1, 2]], cc = p[[1, 3]], 
                    da = p[[2, 1]], db = p[[2, 2]], 
                    dc = p[[2, 3]], i = 0, x, y},While[i++;x = (a^2) + 2 b c + a k da  + ca;y = 2 a b + (c^2) + c k db + cb;c = (b^2) + 2 a c + b k dc + cc;
a = x;b = y;
i < it&& (((a - b))^2)+(((b - c))^2)+(((c - a))^2)<4&&Abs[a + b + c] < 6];-i]];

cDavMJ2 = 
        Compile[{a0, b0, c0, 
            d0, {p, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p[[1, 1]], cb = p[[1, 2]], 
              cc = p[[1, 3]], cd = p[[1, 4]], ka = p[[2, 1]], 
              kb = p[[2, 2]], kc = p[[2, 3]], kd = p[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + ka a + a^2 - b^2 - c^2 + d^2, 
                  cb + kb b + 2*a*b - 2  c*d, cc + kc c + 2*a*c - 2  b*d, 
                  cd + kd d + 2 b c + 2 a d}; 
a = 
                q[[1]]; b = q[[2]]; c = q[[3]]; 
              d = q[[4]]; (((b + c))^2) + (((a - d))^2) + (((b - c))^2) + (((a + d))^2) < 4 && i < it]; -i]];

\(\(\( (*Arbitrary\ Stability\ \(Limit!\)\ *) \)\(
\)\(\(cKMJ2 = 
        Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, c = c0, d = d0, ca = p[\([1, 1]\)], 
              cb = p[\([1, 2]\)], cc = p[\([1, 3]\)], cd = p[\([1, 4]\)], 
              ka = p[\([2, 1]\)], kb = p[\([2, 2]\)], kc = p[\([2, 3]\)], 
              kd = p[\([2, 4]\)], i = 0, v, w, x, y}, 
            While[\(i++\); v = ca + ka\ a + a\^2 + b\^2 + c\^2 + d\^2; 
              w = cb + kb\ b + 2\ a\ b + 2\ c\ d; 
              x = cc + kc\ c + 2\ a\ c + 2\ b\ d; 
              y = cd + kd\ d + 2\ b\ c + 2\ a\ d; 
i < it && 
                Abs[a] < 2 && Abs[b] < 2 && Abs[c] < 2 && Abs[d] < 2, a = v; 
              b = w; c = x; 
              d = y]; \(-i\)]];\)

 (*Arbitrary\ Stability\ \(Limit!\)\ *) 

    \(cKMJ3 = 
        Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, c = c0, d = d0, ca = p[\([1, 1]\)], 
              cb = p[\([1, 2]\)], cc = p[\([1, 3]\)], cd = p[\([1, 4]\)], 
              ka = p[\([2, 1]\)], kb = p[\([2, 2]\)], kc = p[\([2, 3]\)], 
              kd = p[\([2, 4]\)], i = 0, v, w, x, y}, 
            While[\(i++\); 
              v = ca + ka\ a + a\^3 + 3\ a\ b\^2 + 3\ a\ c\^2 + 6\ b\ c\ d + 
                  3\ a\ d\^2; 
              w = cb + kb\ b + 3\ a\^2\ b + b\^3 + 3\ b\ c\^2 + 6\ a\ c\ d + 
                  3\ b\ d\^2; 
              x = cc + kc\ c + 3\ a\^2\ c + 3\ b\^2\ c + c\^3 + 6\ a\ b\ d + 
                  3\ c\ d\^2; 
              y = cd + kd\ d + 6\ a\ b\ c + 3\ a\^2\ d + 3\ b\^2\ d + 
                  3\ c\^2\ d + d\^3; 
i < it && Abs[a] < 2 && 
                Abs[b] < 2 && Abs[c] < 2 && Abs[d] < 2, a = v; b = w; c = x; 
              d = y]; \(-i\)]];\)\)\)\)

cC3KcMJ2 = 
        Compile[{a0, b0, c0, 
            d0, {p, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p[[1, 1]], cb = p[[1, 2]], 
              cc = p[[1, 3]], cd = p[[1, 4]], ka = p[[2, 1]], 
              kb = p[[2, 2]], kc = p[[2, 3]], kd = p[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ka a + ca + a^2 - b^2 + c^2 - d^2, 
                  kb b + cb + 2*a*b + 2  c*d, kc c + cc + 2*a*c - 2  b*d, 
                  kd d++ cd + 2 b c + 2 a d}; 
a = 
                q[[1]]; b = q[[2]]; c = q[[3]]; 
              d = q[[4]]; (((a - c))^2) + (((b - d))^2) + (((a + c))^2) + (((b + d))^2) < 4 && i < it]; -i]];

cC2MJ2 = 
        Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, i = 0, v, w, ca = p[[1, 1]], 
              cb = p[[1, 2]], ka = p[[2, 1]], kb = p[[2, 2]]}, 
            While[i++; v = (a^2) + (b^2) + ka a + ca; 
              w = 2 a b + kb b + cb; 
i < it && 
                Abs[x] < 50 && Abs[a + b] < 2 && Abs[a - b] < 2, a = v; 
              b = w]; -i]];

\(\(\(cC4MJ2 = 
        Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, c = c0, d = d0, i = 0, v, w, x, y, 
              ca = p[\([1, 1]\)], cb = p[\([1, 2]\)], cc = p[\([1, 3]\)], 
              cd = p[\([1, 4]\)], ka = p[\([2, 1]\)], kb = p[\([2, 2]\)], 
              kc = p[\([2, 3]\)], kd = p[\([2, 4]\)]}, 
            While[\(i++\); v = k\ a + ca + a\^2 + c\^2 + 2\ b\ d; 
              w = k\ b + cb + 2\ a\ b + 2\ c\ d; 
              x = kc\ c + cc + b\^2 + 2\ a\ c + d\^2; 
              y = kd\ d + cd + 2\ b\ c + 2\ a\ dkd\ d + 
                  cd; 
i < it && Abs[x] < 50 && 
                v\^2 + w\^2 + x\^2 + y\^2 < 2, a = v; b = w; c = x; 
              d = y]; \(-i\)]];\)\(
\)
    \)\)
cC4MJ3 = 
        Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, c = c0, d = d0, i = 0, v, w, x, y, 
              ca = p[[1, 1]], cb = p[[1, 2]], cc = p[[1, 3]], 
              cd = p[[1, 4]], ka = p[[2, 1]], kb = p[[2, 2]], 
              kc = p[[2, 3]], kd = p[[2, 4]]}, 
            While[i++; 
              v = ka a + ca + (a^3) + 3 (b^2) c + 3 a (c^2) + 6 a b d + 
                  3 c (d^2); 
              w = kb b + cb + 3 (a^2) b + 3 b (c^2) + 3 (b^2) d + 
                  6 a c d + (d^3); 
              x = kc c + cc + 3 a (b^2) + 3 (a^2) c + (c^3) + 6 b c d + 
                  3 a (d^2); 
              y = kd d + cd + (b^3) + 6 a b c + 3 (a^2) d + 3 (c^2) d + 
                  3 b (d^2); 
i < it && Abs[x] < 50 && 
                (v^2) + (w^2) + (x^2) + (y^2) < 2, a = v; b = w; c = x; 
              d = y]; -i]];

cCL2MJ2 = 
        Compile[{a0, b0, c0, 
            d0, {p, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p[[1, 1]], cb = p[[1, 2]], 
              cc = p[[1, 3]], cd = p[[1, 4]], ka = p[[2, 1]], 
              kb = p[[2, 2]], kc = p[[2, 3]], kd = p[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + ka a + a^2 + b^2 - c^2 + d^2, cb + kb b + 2*a*b, 
                  cc + kc c + 2*a*c, 
                  cd + kd d + 2 a d}; 
a = q[[1]]; 
              b = q[[2]]; c = q[[3]]; d = q[[4]]; 
              Abs[a] < 2 && Abs[b] < 2 && Abs[c] < 2 && Abs[d] < 2 && 
                Abs[(a^2) - (b^2) + (c^2) - (d^2)] < 2 && 
                i < it]; -i]];

    
cCL2MJ3 = 
        Compile[{a0, b0, c0, 
            d0, {p, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p[[1, 1]], cb = p[[1, 2]], 
              cc = p[[1, 3]], cd = p[[1, 4]], ka = p[[2, 1]], 
              kb = p[[2, 2]], kc = p[[2, 3]], kd = p[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + ka a + (a^3) + 3 a (b^2) - 3 a (c^2) + 3 a (d^2), 
                  cb + kb b + 3 (a^2) b + (b^3) - b (c^2) + b (d^2), 
                  cc + kc c + 3 (a^2) c + (b^2) c - (c^3) + c (d^2), 
                  cd + kd d + 3 (a^2) d + (b^2) d - (c^2) d + 
                    (d^3)}; 
a = q[[1]]; b = q[[2]]; 
              c = q[[3]]; d = q[[4]]; 
              Abs[a] < 2 && Abs[b] < 2 && Abs[c] < 2 && Abs[d] < 2 && 
                Abs[(a^2) - (b^2) + (c^2) - (d^2)] < 2 && 
                i < it]; -i]];

\(\(\( (*No\ fractals\ found . \ 
          Not\ \(\(initialized\)\(.\)\)*) \)\(
\)\(\(cKrMJ2 = Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
          Module[{a = a0, b = b0, c = c0, d = d0, i = 0, v, w, x, y, 
              ca = p[\([1, 1]\)], cb = p[\([1, 2]\)], cc = p[\([1, 3]\)], 
              cd = p[\([1, 4]\)], ka = p[\([2, 1]\)], kb = p[\([2, 2]\)], 
              kc = p[\([2, 3]\)], kd = p[\([2, 4]\)]}, 
            While[\(i++\); v = a\^2 + b\^2 + c\^2 + d\^2 + ka\ a + ca; 
              w = 2\ a\ b - 2\ c\ d + kb\ b + cb; 
              x = 2\ a\ c - 2\ b\ d + kc\ c + cc; 
              y = \(-2\)\ b\ c + 2\ a\ d + kd\ d + cd; 
i < 
                  it && Abs[x] < 50 && Abs[a - b - c - d] < 2 && 
                Abs[a + b + c - d] < 2 && Abs[a + b - c + d] < 2 && 
                Abs[a - b + c + d] < 2, a = v; b = w; c = x; 
              d = y]; \(-i\)]];\)


    Clear[
      cKrMJ3];  (*No\ fractals\ found . \ 
          Not\ \(\(initialized\)\(.\)\)*) 
cKrMJ3 = 
      Compile[{a0, b0, c0, d0, {p, _Real, 2}, {it, _Integer}}, 
        Module[{a = a0, b = b0, c = c0, d = d0, i = 0, v, w, x, y, 
            ca = p[\([1, 1]\)], cb = p[\([1, 2]\)], cc = p[\([1, 3]\)], 
            cd = p[\([1, 4]\)], ka = p[\([2, 1]\)], kb = p[\([2, 2]\)], 
            kc = p[\([2, 3]\)], kd = p[\([2, 4]\)]}, 
          While[\(i++\); 
            v = a\^3 + 3\ a\ b\^2 + 3\ a\ c\^2 - 6\ b\ c\ d + 3\ a\ d\^2 + 
                ka\ a + ca; 
            w = 3\ a\^2\ b + b\^3 + 3\ b\ c\^2 - 6\ a\ c\ d + 3\ b\ d\^2 + 
                kb\ b + cb; 
            x = 3\ a\^2\ c + 3\ b\^2\ c + c\^3 - 6\ a\ b\ d + 3\ c\ d\^2 + 
                kc\ c + cc; 
            y = \(-6\)\ a\ b\ c + 3\ a\^2\ d + 3\ b\^2\ d + 3\ c\^2\ d + 
                d\^3 + kd\ d + cd; 
i < it && Abs[x] < 50 && 
              Abs[a - b - c - d] < 2 && Abs[a + b + c - d] < 2 && 
              Abs[a + b - c + d] < 2 && Abs[a - b + c + d] < 2, a = v; b = w; 
            c = x; d = y]; \(-i\)]];\)\)\)

cC4C2rMJ2 = 
        Compile[{a0, b0, c0, 
            d0, {p, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p[[1, 1]], cb = p[[1, 2]], 
              cc = p[[1, 3]], cd = p[[1, 4]], ka = p[[2, 1]], 
              kb = p[[2, 2]], kc = p[[2, 3]], kd = p[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + ka a + (a^2) + (c^2) - 2 b d, 
                  cb + kb b + 2 a b - 2 c d, 
                  cc + kc c + (b^2) + 2 a c + (d^2), 
                  cd + kd d - 2 b c + 2 a d}; 
a = 
                q[[1]]; b = q[[2]]; c = q[[3]]; 
              d = q[[4]]; (((a - c))^2) + (((b + d))^2) < 2 && 
                Abs[a + b + c - d] < 2 && Abs[a - b + c + d] < 2 && 
                i < it]; -i]];

cP4MJ2 = 
        Compile[{a0, b0, c0, 
            d0, {p, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p[[1, 1]], cb = p[[1, 2]], 
              cc = p[[1, 3]], cd = p[[1, 4]], ka = p[[2, 1]], 
              kb = p[[2, 2]], kc = p[[2, 3]], kd = p[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + ka a + (a^2) + (b^2) + (c^2) + (d^2), 
                  cb + kb b + 2 a b, cc + kc c + 2 a c, 
                  cd + kd d + 2 a d}; 
a = q[[1]]; 
              b = q[[2]]; c = q[[3]]; d = q[[4]]; 
              Abs[(a^2) - (b^2) - (c^2) - (d^2)] < 2 && (a^2) < 8 && (b^2) < 8 && 
                (c^2) < 8 && (d^2) < 8 && i < it]; -i]];

cP4MJ3 = 
        Compile[{a0, b0, c0, 
            d0, {p, _Real, 
              2}, {it, _Integer}}, 
Module[{a = a0, b = b0, 
              c = c0, d = d0, i = 0, ca = p[[1, 1]], cb = p[[1, 2]], 
              cc = p[[1, 3]], cd = p[[1, 4]], ka = p[[2, 1]], 
              kb = p[[2, 2]], kc = p[[2, 3]], kd = p[[2, 4]], 
              q = {a0, b0, c0, d0}}, 
While[i++; 
              q = {ca + ka a + (a^3) + 3 a (b^2) + 3 a (c^2) + 3 a (d^2), 
                  cb + kb b + 3 (a^2) b + (b^3) + b (c^2) + b (d^2), 
                  cc + kc c + 3 (a^2) c + (b^2) c + (c^3) + c (d^2), 
                  cd + kd d + 3 (a^2) d + (b^2) d + (c^2) d + 
                    (d^3)}; 
a = q[[1]]; b = q[[2]]; 
              c = q[[3]]; d = q[[4]]; 
              Abs[(a^2) - (b^2) - (c^2) - (d^2)] < 2 && (a^2) < 8 && (b^2) < 8 && 
                (c^2) < 8 && (d^2) < 8 && i < it]; -i]];

Usage.

boundary::usage="boundary[t_:txy,a_:4] works through t to identify and output near-boundary points as a list of lattice indices. Also,a list of cell coordinates (calculated using global values X0,Xd Y0,Yd, set up when creating a txy list) is stored as ptb. Points that have passed the stability test are black. (They may be far from the true stable set boundary if it is low - see examples 13 to 15.)A near-boundary point is black and has at least one non-black near neighbour. If a≤4, only the edge-neighbours are tested, giving the smallest list; if a>4 the corner neighbours are also tested; if a>8 another 12 neighbours (within a (Sqrt[5]) circle) are tested: \!\(\*GridBox[{
{\(\\ \), \(x\), \(x\), \(x\), \(\\ \)},
{\(x\), \(8\), \(4\), \(8\), \(x\)},
{\(x\), \(4\), \(0\), \(4\), \(x\)},
{\(x\), \(8\), \(4\), \(8\), \(x\)},
{\(\\ \), \(x\), \(x\), \(x\), \(\\ \)}
}]\).
ptb and a>4 are not used in this notebook.";

cJulia::usage = "See GenMJArray::usage";

cMandelbrot::usage="The compiled function called by GenMJArray to plot Mandelbrot sets (Subscript[A, n + 1])=(Subscript[A, n]^2)+(Subscript[A, n]). The arguments are the same as for cJulia, but vi, wi & p are not used and array plots are irrelevant.";

BW::usage = "A colorfunction that shows the points that have \
reached the iteration in black, others in white.";

cQMJ2::usage="The compiled function called by GenMJArray to plot Quaternion Mandelbrot-Julia (square) sets \!\(\*
StyleBox[\(A\_\(n\\  + 1\)\),
FontSlant->"Plain"]\)=(Subscript[A, n]^2)+k (Subscript[A, n])+ C, where A={a,b,c,d} and C={ca,cb,cc,cd}are quaternions and k is a set of four reals that allow the (Subscript[A, n]) term to be a conjugate. Iterators vi & wi define c & d, xi & yi define a & b, so array plots are slices of the 4D function. p is {{ca,cb,cc,cd},{ka,kb,kc,kd}}, the Julia quaternion constant C and the extended Mandelbrot factors k";

\(cQMJ3::usage = \*""The compiled function called by GenMJArray to plot the generalised Quaternion-cube sets (Subscript[A, n + 1])=(Subscript[A, n]^3)+k (Subscript[A, n])+C. See cQMJ2"";\)

cC4C2rMJ2::usage = "The compiled function called by GenMJArray to \
plot sets for the C4C2r algebra, for which {a,b,c,d},{e,f,g,h}=
{a e-d f+c g-b h, b e+a f-d g-c h, c e+b f+a g+d h, d e-c f-b g+a \
h}";

GenMJArray::usage="GenMJArray[func_:cJulia,vi_:{-1.5,3/4,1/4},wi_:{-4/3,1/3,1/3},xi_:{-1.1,1.1,0.06},yi_:{-1.1,1.1,0.06},p_:2.,it_:35,colf_:Automatic,conum_:10] creates Mandelbrot, Julia and similar stable iterative sets. The default parameters generate rows of Dickau's [1] array of purported Julia sets (only the central plots are stable sets). The lower rows are central symmetries of the upper rows  so all but one are omitted. The first parameter is a function name  - cJulia is a compiled function that iterates \!\(\*
StyleBox[\(A\_\(\(n\\  + \\ 1\)\(=\)\)\),
FontSlant->"Plain"]\)(Subscript[A, n]^p)+C for up to it iterations or until the length of the complex vector A exceeds 2 (i.e. is outside the circle of stability). vi & wi are iterators for v & w,giving C=v+Iw (A single ListDensityplot (with contours->conum) is produced if vi & wi are single numbers rather than 3-member lists; global values F0=Hold[func];P0=p;Vi=vi;Wi=wi;X0=xi[[1]];Y0=yi[[1]];Xd=xi[[3]];Yd=yi[[3]];itlim=-it-1  are also created in this case).
xi & yi are iterators for x & y giving starting values (Subscript[A, 0])=x+Iy for the plots. p=2 gives the standard Julia set; other powers (excluding 0 & 1) can be used  but their circles of stability may exceed 2. Other functions interpret these arguments differently - see the function usages. Plots are suppressed if conum<2.";

HueR::usage = "The colorFunction to give a reversed-colour \
plot.";
HueB::usage = "The colorFunction to give a colour plot with black \
minimal points.";
HueW::usage = "The colorFunction to give a colour plot with white \
minimal points."; 
    Hue4::usage = "The colorFunction to give a colour plot with black \
minimal points whilst subdividing the coloured areas into 4 repeated colour \
zones."; 
    Hue10::usage = "The colorFunction to give a colour plot with black \
minimal points whilst subdividing the coloured areas into 10 repeated colour \
zones.";

3. Real & Complex Algebra Demonstrations.

RAbs

The simplest non-linear iteration is RAbs (Abs for Real numbers). It gives cyclic/chaotic behaviour by iterating w=p Abs[y]-x+a; x=y; y=Chop[w]. When p=-1 there are cyclic regions (grey) resembling a rabbit, embedded in a chaotic zone and surrounded by more cyclic zones. When p=+1, the central cyclic zones resemble a spaceman. Other p values scale the figures.

GenMJArray[cRAbs, 1, 
        0, {4.1, -5, - .211}, {4.1, -5, - .211}, -1, 
        180]//Timing

GenMJArray[cRAbs, 1, 
        0, {8, -5.1, - .211}, {8, -5.1, - .211}, 1, 
        180]//Timing

Dickau Example, Julia Sets.

Ex. 3. is Dickau's Julia example [4], iterating (Subscript[J, n + 1]) = (Subscript[J, n]^2) + Cin the complex plane. The lower rows are reflections of the upper rows, so all but one are omitted. The default case is:-
GenMJArray[func_:cJulia,vi_:{-1.5,3/4,1/4},wi_:{4/3,-1/3,-1/3},
  xi_:{-1.1,1.1,0.06},yi_:{-1.1,1.1,0.06},p_:2,it_:35]

GenMJArray[];

Many of these plots are strictly outside the Julia Set, for they reach an unstable point in a finite number of iterations. Ex. 4 repeats the fourth row, sixth column with it=200; -Min[txy] is the largest number of iterations for any point. As it is 188, rather than 201, no points in this plot are within the set. The black "hubs" arise because the number of contours for ListDensityPlot cannot be changed from 10.

(*Ex .4 .  Single high resolution plot with it=200 shows no points in the Julia set for k=-2.  *)
GenMJArray[cJulia, 0, 
          2/3, {-1.2, 1.2, 0.01}, {-1.15, 1.15, 0.01}, 2, 200, 
          HueB];

      -Min[txy]

The following single plots are analysed by boundary[] to create a compact representation of the boundary of the "pseudo-set" (those points that pass the test for it iterations). The plots have the {x,y} ranges adjusted to leave a narrow band outside the set limits because boundary[] accesses one point beyond the extreme pseudo-set points.
Two named fractal Julia sets [3] (demonstrating the use of HueB & HueR), and their boundary plots:-

\(\(\( (*Ex  .5 . \ Douady' s\ Rabbit*) \)\(GenMJArray[
      cJulia, \(- .123\),  .745, {\(-1.3\), 1.3, 0.01}, {\(-1.1\), 1.1, 
        0.01}, 2, 50, HueB] (*Douady' s\ Rabbit*) 

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*Ex  .6 . \ San\ Marco*) \)\(\(GenMJArray[cJulia, \(- .75\), 
        0, {\(-1.7\), 1.7, 0.01}, {\(- .9\),  .9, 0.01}, 2, 50, HueR];\)

    \(ListPlot[boundary[txy, 4], PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

In Ex.7, the lower half is omitted as it is a reflection of the upper half through the centre-point. Ex.8 magnifies the right-hand tip of this plot and demonstrates the effect of HueW. Only the boundary point plots are retained.

\(\(\( (*\ 
      Ex  .7 . \ Single\ Julia\ Half - 
        plot\ with\ boundary\ \(\(plot\)\(.\)\)\ *) \)\(\(GenMJArray[
        cJulia,  .3, \(- .5\), {\(- .65\), 1.1, 0.01}, {0, 1.2, 0.01}, 
        2, \(-35\)];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*\ 
      Ex  .8 . \ Enlarged\ right\ hand\ extreme\ of\ Ex . \ 
          7*) \)\(\(GenMJArray[
        cJulia,  .3, \(- .5\), { .95, 1.05, 0.0005}, { .5,  .65, 0.0005}, 
        2, \(-35\)];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

Mandelbrot Set.

There are two, equivalent, expressions for the Mandelbrot set, (Subscript[M, n + 1]) = (Subscript[M, n]^2) + (Subscript[M, 0]) &  (Subscript[M, n + 1]) = (Subscript[M, n]^2) + (Subscript[M, n]). The second is preferred as it allows a generalization (next section.) The Julia sets have 4 parameters {v,w,x,y}; I included an extra parameter p so that \(TraditionalForm\`\((x + iy)\)\_\(n + 1\) = \((x + iy)\)\_n\%p + v + iw\) can be explored. The Mandelbrot set  has only two parameters {x,y}, so v, w & p are not used; they are included as dummy parameters in the general procedure.
In Ex.9, the y=0 line is included by careful choice of the initial y value and step length; otherwise the extension to x=-2 could be missed. The x step is chosen to find the blob at {-1.77, ±.016}, but misses the smaller one at {-1.477, ±.005}; these are shown as Ex. 10 & Ex. 11. As the boundary is continuous [2] any isolated points must be connected to the central area by unplottably fine filaments.  Examples 12 to 15 (together with Ex.9) reproduce Fig.1.2 in Penrose[5]. Many iterations are required to approximate the true boundary set at high resolution, and the mesh points rarely coincide with these filaments.

\(\(\( (*Ex  .9 . \ MandelBrot . \ Penrose\ 2004\ Fig . \ 1.2  
        a\ *) \)\(\(GenMJArray[cMandelbrot, 0, 
        0, {\(-2.256\),  .65,  .007325}, {\(-1.125\), 1.125, 1/128. }, 
        0, \(-499\)];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*\ 
      Ex  .10 . \ Mandelbrot\ Set . \ Upper\ half\ of\ mini - 
        blob\ at\ \(\({\(-1.477\), 0}\)\(.\)\)\ *) \)\(\(GenMJArray[
        cMandelbrot, 0, 
        0, {\(-1.483\), \(-1.473\),  .0001}, {\(- .0001\),  .0046,  .0001}, 
        2, \(-999\)];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*\ 
      Ex  .11 . \ Mandelbrot\ Set . \ Upper\ half\ of\ larger\ mini - 
        blob\ at\ \(\({\(-1.77\), 0}\)\(.\)\)\ *) \)\(\(GenMJArray[
        cMandelbrot, 0, 
        0, {\(-1.78\), \(-1.745\),  .0005}, {\(- .0005\),  .017,  .0005}, 
        2, \(-999\)];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

The other Mandelbrot sets illustrated in [5] are reproduced (Slow - many minutes):-

GenMJArray[cMandelbrot, 0, 
        0, {- .6, - .4,  .0005}, { .5,  .7,  .0005}, 0, 1999, 
        HueB];

Ex.12a and 12b demonstrate the trade-off between low and high iteration limits, by magnifying the "blob" near the top left corner of the above figure. A limit of 100 picks up some of the filaments, at the cost of filling in the valleys.

GenMJArray[cMandelbrot, 0, 
        0, {- .6005, - .595,  .00001}, { .6605,  .6655,  .00001}, 0, 
        100, BW];

GenMJArray[cMandelbrot, 0, 
        0, {- .6005, - .595,  .00001}, { .6605,  .6655,  .00001}, 0, 
        1000, BW];

GenMJArray[
        cMandelbrot, 0, 
        0, {- .797, - .7865,  .00002}, { .14,  .152,  .00002}, 2, 
        500, HueW];

GenMJArray[cMandelbrot, 
        0, 0, {- .797, - .7865,  .00001}, { .14,  .152,  .00001}, 2, 
        2000, BW];

\(\(\( (*Ex  .15 . \ Penrose\ 2004\ Fig  .1  .2  d, \ 
      100  iteration\ limit\ shows\ "hubs"\ that\ are\ regions\ of\ very\ slow\ divergence\ and\ appear\ black\ when\ a\ low\ iteration\ limit\ is\ employed . \ These\ shrink\ when\ the\ iteration\ limit\ is\ increased, \ 
      but\ the\ connecting\ filaments\ are\ \(\(lost\)\(.\)\)\ *) \)\(\(GenMJArray[cMandelbrot, 0, 
        0, {\(- .7491\), \(- .7478\),  .000005}, { .148,  .1495,  .000005}, 
        2, 99, HueB];\)

    \(GenMJArray[cMandelbrot, 0, 
        0, {\(- .7491\), \(- .7478\),  .000005}, { .148,  .1495,  .000005}, 
        2, 9999, BW];\)\)\)\)

Mandelbar, Skeletons.

The name "Mandelbar Set" has been applied to the stable initial points of \(TraditionalForm\`M\_\(n + 1\) = \(M\&_\)\_n\%p + M\_0\). The more general iteration \(TraditionalForm\`M\_\(n + 1\) = M\_n\%p + k  M\&_\_n\) gives the same points when k=1. cMandelbar is set up as \(TraditionalForm\`M\_\(n + 1\) = 
        M\_n\%v + \((w + i\ p)\) M\&_\_n\) so that the array shows different powers ({1/2, 1, 3/2, 2, 5/2, 3, 7/2, 4} in Ex.16) horizontally and multipliers (1,1/2 0(boring) -1/2,-1} vertically.
Although p can be used to make k a complex number, the results are messy and none is shown.

GenMJArray[
        cMandelbar, {3/2, 4, 1/2}, {-1, 1, 1/2}, {-2.1, 2.1, 
          0.05}, {-1.7, 1.7, 0.05}, 0, 99];

The bottom row of Ex.16 (which extends the Mandelbar examples in [2]) shows artefacts resulting from a low iteration limit and/or insufficient contours. The stable set is not an area with a boundary, but a set of branching lines (a "skeleton" in image-analysis jargon). This is demonstrated for the base case in Ex. 17. As it is increased, the curves at the junctions shrink. Sloping lines may only seen as rows of dots because they rarely coincide with grid points.

\(\(\( (*Ex  .17 . \ 
          Basic\ Mandelbar\ is\ a\ set\ of\ \(lines!\)\ *) \)\(\(GenMJArray[
        cMandelbar, 2, 1, {\(-2. \), 1.4, 0.01}, {\(-1.7\), 1.7, 0.01}, 0, 
        19];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*Ex  .18 . \ 
          Mandelbar\ w1th\ negative\ k\ gives\ areas*) \)\(\(GenMJArray[
        cMandelbar, 3/2, \(-1\), {\(-2.5\), 3, 0.005}, {\(-2\), 2, 0.005}, 0, 
        49, HueB];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*Ex  .19 . \ 
          Cubic\ Mandelbar\ is\ a\ set\ of\ \(\(lines\)\(.\)\)\ *) \)\(\(GenMJArray[cMandelbar, 3, 
        1, {\(-2.1\), 2.1, 0.02}, {\(-1.7\), 1.7, 0.02}, 0, 49, HueB];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

Mandelbrot-Julia, MJ.

I generalize the Mandelbrot iteration by making the last term k (Subscript[M, n]) (in place of (Subscript[M, 0])), including p as a power, and adding the Julia constant, to give a generalization, G, (Subscript[G, n + 1])=(Subscript[G, n]^p)+k (Subscript[G, n]) +C. The Mandelbrot and Julia set relationships are simplifications with p=2 and either C =0 or k =0. (For some other algebras, I use k to give the conjugate \(TraditionalForm\`G\&_\_n\).). 
In the compiled version,  p has to be a real matrix (e.g.{{2.,2.}}) to provide compilable parameters,  p & k.
The factorized p=2 version (Subscript[G, n + 1])=(Subscript[G, n]^2)+(a+b) (Subscript[G, n]) +ab is provided as cMJ2factor.
I have not explored these functions in detail

GenMJArray[
        cMJpk, {-3/4, 1/2. , 1/4}, {-3/4, 0, 1/4}, {-2.1,  .1, 
          0.05}, {-1.5, 1.5, 0.05}, {{2. , 2. }}, 35];

\(\(\( (*Ex  .21\ MandelbrotJulia\ Set\ single\ Fractal\ \(\(image\)\(.\)\)\ *) \)\(\(GenMJArray[
        cMJpk,  .3,  .5, {\(-2.1\),  .1, 0.005}, {\(-1.15\), 1.15, 
          0.005}, {{2. , 2. }}, 35];\)

    \(ListPlot[
boundary[], 
        PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

GenMJArray[
        cMJ2factor, {-1/4, 1/4. , 1/8}, {-1/2, 1/8, 1/8}, {-2.6, 
          1.3, 0.05}, {-1.5, 2, 0.05}, {{2.1, - .9}}, 
        35];

GenMJArray[cMJ2factor, 
        1/4, -1/4, {-2.1, - .4, 0.005}, {- .1, 2.1, 
          0.005}, {{2. , -1. }}, 99];

GenMJArray[
        cMJpk, {0, 3/4, 1/4}, {-1/2, 1/2, 1/4}, {- .9,  .9, 
          0.05}, {-2.2, 2.2, 0.05}, {{3. , 2. }}, 35];

\(\(\( (*Ex  .25 . \ 
          Fractal\ Cubic\ MandelbrotJulia\ Set\ with\ Fatou\ \(\(dust\)\(.\)\)*) \)\(\(GenMJArray[cMJpk, 
        0.1, \(- .1\), {\(- .6\),  .6, 0.005}, {\(-2\), 2, 
          0.005}, {{3. , 2}}, 35];\)

    \(ListPlot[
boundary[], 
        PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*Ex  .26\ Rough\ boundary, \ 
      Cubic\ MandelbrotJulia\ \(\(Set\)\(.\)\)*) \)\(\(GenMJArray[
        cMJpk, \(-7\)/8, \(- .3\), {\(-2\), 1.1, 0.01}, {\(-1.1\), 1.1, 
          0.01}, {{2. , 1. }}, \(-35\)];\)

    \(ListPlot[
boundary[], 
        PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

4. Other Algebras.

GenMJArray handles any iteration (Subscript[A, n + 1]) = F[(Subscript[A, n])] for which a func has been defined, with A being general number - a set of m coefficients that will be mapped onto {x,y..} etc. - it is NOT restricted to complex numbers A=x+i y with (Subscript[A, n]^p)evaluated in the complex plane.Whilst arbitrary functions are fun, my interest is in  algebras. Several different algebras are demonstrated, with appropriate multiplication rules and stability limits obtained from [1]. The power is restricted to 2 or 3 for speed; the general cases of the Mandelbrot-Julia function iterate (Subscript[A, n + 1])= (Subscript[A, n]^2)+k (Subscript[A, n])+ C or (Subscript[A, n + 1])= (Subscript[A, n]^3)+k (Subscript[A, n])+ C. Quaternions are typical of 4-element algebras, with {v,w,x,y} being a quaternion; their multiple plots show a 4D surface as {x,y} slices at a range of {v,w} values. Algebras with more than 4 elements, conjugates, and functions with extra parameters, are handled by packing items (on a case-by-case basis) into p.

Quaternions.

The explicit quaternion square and cube generalized Mandelbrot & Julia sets are supplied, as cQpMJ = (Subscript[A, n + 1])= (Subscript[A, n]^p)+k (Subscript[A, n])+C. Each point has 4 coordinates {a,b,c,d} and extra parameters p[[1, 1 to 5]] = {(Subscript[C, a]),(Subscript[C, b]),(Subscript[C, c]),(Subscript[C, d]), k}. The 4D surface is shown as 2D {a,b} slices on a quadrant of the array of {c,d} values. {b,c,d} symmetries mean that their negative values need not be shown. I have only looked at a few examples. A 3D representation is in[2], but this does not use Mathematica.

GenMJArray[
        cQMJ2, {0, 4/5, 1/5}, {-4/5, 0, 1/5}, {- .52,  .35, 
          0.02}, {- .02, 1.1, 
          0.02}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. }}, 35] // 
      Timing

GenMJArray[
        cQMJ3, {0, 4/5, 1/5}, {-4/5, 0, 1/5}, {- .8,  .7, 0.05}, {0, 
          1.45, 0.05}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. }}, 35] // 
      Timing

\(\(\( (*Ex  .29, \ 
      Quaternion\ Cube\ Mandelbrot\ Fractal\ with\ k = 
        1\ \((Quadrant\ only . \ \(SLOW!\))\)*) \)\(GenMJArray[cQMJ3, 0, 
        0, {\(- .01\),  .75, 0.005}, {\(- .01\), 1.45, 
          0.005}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. }}, \(-35\)] // 
      Timing

    \(ListPlot[
boundary[], 
        PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

GenMJArray[
        cQMJ2, {0, 4/5, 1/5}, {-4/5, 0, 1/5}, {-1.2, 1, 
          0.05}, {-1.1, 1.1, 
          0.05}, {{0. ,  .3, - .1,  .432}, {0. , 0. , 0. , 0. }}, 35] // 
      Timing

GenMJArray[
        cQMJ2, {0, 4/5, 1/5}, {-4/5, 0, 1/5}, {-1. , 1.2, 
          0.05}, {-1.1, 1.1, 
          0.05}, {{- .27,  .03, - .5,  .32}, {0. , 0. , 0. , 0. }}, 
        35] // Timing

GenMJArray[
        cQMJ2, {0, 4/5, 1/5}, {-4/5, 0, 1/5}, {-1.2, 1, 
          0.05}, {-1.1, 1.1, 
          0.05}, {{0. , 0. , 0,  .0}, { .5,  .5,  .5,  .5}}, 35] // 
      Timing

GenMJArray[
        cQMJ2, {0, 1, 1/5}, {-3/5, 0, 1/5}, {-2, 0, 0.05}, {-1, 
          1, 0.05}, {{0. , 0. , 0,  .0}, {2. , 2. , 2. , 2. }}, 35] // 
      Timing

GenMJArray[
        cQMJ2, {0, 3/5, 1/5}, {-3/5, 0, 1/5}, {-1.5, 0, 
          0.05}, {-1, 1, 
          0.05}, {{- .27,  .03, - .15,  .32}, {2. , 2. , 2. , 2. }}, 
        35];

Ex.35 shows that a 2D slice of the cQMJ2 4D surface can have detached areas and so a 3D slice would be difficult to display.

\(\(\( (*Ex  .35\ Quaternion\ Square\ Julia\ Set\ expanded\ \(\(slice\)\(.\)\)*) \)\(GenMJArray[
        cQMJ2,  .57,  .56, {\(- .6\),  .3, 0.01}, {\(- .75\),  .7, 
          0.01}, {{\(- .3\),  .35,  .37,  .22}, {0. , 0. , 0. , 
            0. }}, \(-50\)] // Timing

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

The cQMJ2 and cQMJ3 functions allow the use of four individual k values as multipliers for the individual linear term elements, so that the conjugate ("Mandelbar") function can be shown.

GenMJArray[
        cQMJ2, {0, 3/4, 1/4}, {-3/4, 0, 1/4}, {-1. , 1.5, 
          0.05}, { .05, 
          1.75,  .05}, {{0. , 0. , 0. , 
            0. }, {1. , -1. , -1. , -1. }}, 35];

\(\(\( (*Ex  .37\ Quaternion\ conjugate\ M2\ skeleton\ *) \)\(\(GenMJArray[cQMJ2, 1/8, 
        1/4, {0, 1.1, 0.01}, {0, 1.7, 
          0.01}, {{0. , 0. , 0. , 0. }, {1. , \(-1. \), \(-1. \), \(-1. \)}}, 
        35, HueW];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .02], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*Ex  .38\ Quaternion\ conjugate\ M3\ skeleton\ *) \)\(\(GenMJArray[cQMJ3, 1/8, 
        1/4, {0, 1.05, 0.01}, {0, 1.05, 
          0.01}, {{0. , 0. , 0. , 0. }, {1. , \(-1. \), \(-1. \), \(-1. \)}}, 
        35, HueW];\)

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .02], 
        AspectRatio -> Automatic];\)\)\)\)

Terplex

This is the algebra of the Cyclic 3-element group, C3, with{a,b,c}*{e,f,g}= {a e+c f+b g,b e+a f+c g,c e+b f+a g}. The stability limit is (((a - b))^2)+(((b - c))^2)+(((c - a))^2)>4, a+b+c>2. Four versions, (Subscript[A, n])=(Subscript[A, n]^2)+C (cC3J2) , (Subscript[A, n])=(Subscript[A, n]^2)+k A+C (cC3MJ2), (Subscript[A, n])=(Subscript[A, n]^2)+k Abar+C (cC3MbJ2)and (Subscript[A, n])=(Subscript[A, n]^2)+k A+C (cC3MJ3), are programmed.

Terplex Julia.

GenMJArray[
        cC3J2, {- .2,  .6,  .2}, {- .2,  .6,  .2}, {- .9,  .9, 
          0.1}, {- .9,  .9, 0.1}, {{- .1, - .2}}, 
        99];

\(\(\( (*Ex  .40\ Julia\ Squared\ Terplex {a, b,  .2}\^2 + { .2, \(- .08\), \(- .1\)}\ V . \ SLOW\ *) \)\(GenMJArray[
        cC3J2,  .2,  .2, {\(- .85\),  .85, 0.01}, {\(-1.1\), 1.05, 
          0.01}, {{\(- .08\), \(- .1\)}}, \(-999\)] // Timing

    \(ListPlot[
boundary[], 
        PlotStyle -> PointSize[ .01], 
        AspectRatio -> Automatic];\)\)\)\)

Look at the RH hook in greater detail:-

GenMJArray[
        cC3J2,  .2,  .2, { .57,  .82, 0.002}, {- .08,  .06, 
          0.002}, {{- .08, - .1}}, 999] // Timing

and even greater detail:-

GenMJArray[
        cC3J2,  .2,  .2, { .57,  .64, 0.0005}, {- .07,  .04, 
          0.0005}, {{- .08, - .1}}, 999] // Timing

GenMJArray[
        cC3J2, {-1/5, 4/5, 1/5}, {-3/5, 1/5, 1/5}, {-1, 1, 
          0.05}, {- .9,  .9, 0.05}, {{- .3,  .2}}, 49] // 
      Timing

GenMJArray[
        cC3J2,  .6,  .2, {- .9,  .6, 0.005}, {- .9,  .4, 
          0.005}, {{- .3,  .2}}, 99] // Timing

(*Ex  .45 Bottom left corner from previous plot.*) GenMJArray[
        cC3J2,  .6,  .2, {- .075, - .005, 
          0.0005}, {- .49, - .43, 0.0005}, {{- .3,  .2}}, 99] // 
      Timing

Terplex MandelbrotJulia.

The general MandelbrotJulia function array shows {a,b} plots, with {c,k} varying in the array. p contains 3 C values.

GenMJArray[
        cC3MJ2, {-3/5, 4/5, 1/5}, {-3/2, 1, 1/2}, {-1.1, 
          2,  .05}, {-1.1, 
          1.1,  .05}, {{0. , 0. , 0. }, {1. , 1. , 1. }}, 35] // 
      Timing

GenMJArray[cC3MJ2, 
        1/5, -1, {-1.1, 2,  .01}, {-1.1, 
          1.1,  .01}, {{0. , 0. , 0. }, {1. , 1. , 1. }}, 35] // 
      Timing

GenMJArray[cC3MJ2, 
        1/5, -1, {- .8, - .45,  .001}, {0,  .4,  .001}, {{0. , 
            0. , 0. }, {1. , 1. , 1. }}, 99] // Timing

GenMJArray[
        cC3MJ2, {-2/5, 4/5, 1/5}, {0, 3/2, 
          1/2}, {-1.8,  .75,  .05}, {-1, 
          1,  .05}, {{ .1, - .12,  .27}, {1. , 1. , 1. }}, 35] // 
      Timing

Conjugate Terplex MandelbrotJulia.

The Terplex conjugate of {a,b,c} is {a,c,b}, and so b k db changes to c k db, c k dc changes to b k dc for the "bar" case.

GenMJArray[
        cC3MbJ2, {-3/5, 4/5, 1/5}, {-3/2, 1, 1/2}, {-1.1, 
          2,  .1}, {-1.1, 1.5,  .1}, {{0. , 0. , 0. }, {1. , 1. , 1. }}, 
        35] // Timing

GenMJArray[cC3MbJ2, 
        1/5, 1/2, {-1.2, 
          1,  .01}, {-1.1,  .8,  .01}, {{0. , 0. , 0. }, {1. , 1.5, 
            1. }}, 35] // Timing

GenMJArray[cC3MbJ2, 
        1/5, 1/2, {-1. , - .55,  .001}, { .35,  .8,  .001}, {{0. , 
            0. , 0. }, {1. , 1.5, 1. }}, 35] // Timing

    ListPlot[
boundary[], 
        PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];

K.

This is the algebra of the Klein group with {a,b,c,d}*{e,f,g,h} = {a e+b f+c g+d h,b e+a f+d g+c h,c e+d f+a g+b h,d e+c f+b g+a h}. The plots are not fractal.

GenMJArray[
        cKMJ2, {0, 3/4, 1/4}, {-3/4, 0, 1/4}, {-1.5,  .5, 
          0.05}, {-1, 1, 
          0.05}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. } }, 
        49];

GenMJArray[
        cKMJ3, {0, 3/4, 1/4}, {-3/4, 0, 1/4}, {-1, 1, 0.05}, {-1, 
          1, 0.05}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. } }, 
        49];

GenMJArray[
        cKMJ2, {0, 1, 1/4}, {-3/4, 0, 1/4}, {-1, 1, 0.05}, {-1, 
          1, 0.05}, {{0.1, 0.2, 0.3, -0.321}, {0. , 0. , 0. , 0. } }, 
        49];

Other 4 element algebras.

Davenport

The Daveport product is:-

\(genTimes[{a, b, c, d}, {e, f, g, h}, "Dav"] = 
{a\ e - b\ f - 
          c\ g + d\ h, b\ e + a\ f - d\ g - c\ h, c\ e - d\ f + a\ g - b\ h, 
        d\ e + c\ f + b\ g + a\ h}\)

GenMJArray[
        cDavMJ2, {0, 3/4, 1/8}, {-3/4, 0, 1/8}, {-1.1, 1, 
          0.1}, {-1, 1, 
          0.1}, {{- .3,  .35, - .1,  .22}, {1. , 1. , 1. , 1. }}, 
        15];

GenMJArray[cDavMJ2, 
        1/4, -1/8, {-1.2,  .8, 0.01}, {- .9,  .9, 
          0.01}, {{- .3,  .35, - .1,  .22}, {1. , 1. , 1. , 
            1. }}, -35] // Timing

    ListPlot[
boundary[], 
        PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];

GenMJArray[
        cDavMJ2, {0, 3/4, 1/8}, {-5/8, 0, 1/8}, {-1.3,  .3, 
          0.05}, {-1.2, 1.2, 
          0.05}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. }}, 15] // 
      Timing

\(\(\( (*Ex  .59, \ 
      Davenport\ Mandelbrot\ Quadrant\ with\ k = 1*) \)\(GenMJArray[cDavMJ2, 
        1/4, 1/8, {\(-1.3\),  .3, 0.01}, {\(- .9\),  .9, 
          0.01}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. }}, \(-35\)] // 
      Timing

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .005], 
        AspectRatio -> Automatic];\)\)\)\)

\(\(\( (*Ex  .60, \ 
      Davenport\ Mandelbrot\ tip\ with\ k = 1*) \)\(GenMJArray[cDavMJ2, 1/4, 
        1/8, { .13,  .23, 0.001}, { .55,  .7, 
          0.001}, {{0. , 0. , 0. , 0. }, {1. , 1. , 1. , 1. }}, \(-35\)] // 
      Timing

    \(ListPlot[boundary[], PlotStyle -> PointSize[ .001], 
        AspectRatio -> Automatic];\)\)\)\)

Clifford(2)

As the conserved size is (a^2) - (b^2) + (c^2) - (d^2), the variables can increase without bound if this is used as the stability criterion. They are arbitrarily limited to ±2.

GenMJArray[
        cCL2MJ2, {0, 3/2, 1/8}, {-3/4, 0, 1/8}, {-1.5, 1.5, 
          0.1}, {-1.5, 1.5, 
          0.1}, {{- .3,  .35, - .1,  .22}, {0. , 0. , 0,  .0}}, 35] // 
      Timing

GenMJArray[cCL2MJ2, 
        11/8, -1/8, {-2, 2, 0.01}, {-2.1, 2.1, 
          0.01}, {{- .3,  .35, - .1,  .22}, {0. , 0. , 0,  .0}}, 35] // 
      Timing

C4C2r

GenMJArray[
        cC4C2rMJ2, {0, 5/8, 1/8}, {-5/8, 0, 1/8}, {-1.1, 1, 
          0.05}, {-1, 1, 
          0.05}, {{- .3,  .35, - .1,  .22}, {0. , 0. , 0. , 0. }}, 
        35] // Timing

GenMJArray[cC4C2rMJ2, 
        1/8, -1/4, {- .4, - .15, 0.001}, {0,  .3, 
          0.001}, {{ .15,  .3, - .1,  .22}, {0. , 0. , 0. , 0. }}, 35] // 
      Timing

References.

[1] http://library.wolfram.com/infocenter/MathSource/4894/

[2] http://mathworld.wolfram.com/MandelbrotSet.html

[3] http://mathworld.wolfram.com/JuliaSet.html

[4] \(\(\(\(http\)\(:\)\) // \(\(mathforum . org/advanced\)/
                    robertd\)/julias . html\)\+_\)

[5] Roger Penrose, The Road To Reality, Jonathon Cape 2004.

[6] \(\(\(\(http\)\(:\)\) // 
                www . physcip . uni - \(stuttgart . de/phy11733\)/
                    index_e . html\)\+_\)