# Maximizing efficiency of a power-tower (tetration) calculation

I'm working on a power tower $$\displaystyle z^{z^{.^{.^{.}}}}$$ and wish to maximize the efficiency of my algorithm to compute a cycle map. This is a color-coded diagram of the cycle period of each point in a region of the complex plane mapped by the power tower. For example, if I start with $$z=-0.727+0.108i$$, the tetration enters a 6-cycle period:

{-0.727 + 0.108 I, -0.540034 - 0.726564 I,
1.83707 - 10.2354 I, -8.29029*10^12 + 8.0833*10^12 I, 0. + 0. I,
1. + 0. I, -0.727 + 0.108 I}


This point is given a value of $$6$$ which is assigned a color by the function Colorize. If the tetration grows without bounds, it's flagged before it reaches the max machine number and assigned a value of -1 or black by Colorize.

And although NestList can implement the tetration, I'm not able to implement the necessary checking of underflow and overflow in a compiled routine using NestList so I coded it directly (below).

Overflow is flagged when $$Re[z_{i+1}\log(z_i)]>500$$ and given a color code of -1 (Black). (Approximate) Underflow, when $$Re[z_{i+1}\log(z_i)]<-500$$. However underflow yields $$z_{i}^0=1, z_{i+1}=z^1$$ and the cycle repeats. The cycle can be quite long and may not appear until a sufficiently large number of steps are generated so I have to generate a relatively long list then delete the duplicates.

My problem is that's it's quite slow. Processing {{x,-1,0,0.001},{y,-1,0,0.001}} took 20 minutes on a dual-core 2.2 GHz machine. I was wondering if someone could look at my code and suggest ways of improving it? I don't have a C-compiler on my machine so can't compile to C.

getTetrationList =
Compile[{{z, _Complex}, {max, _Integer}},
Module[{znext, i, highFlag, myTable = Table[1 + I, {max}]},
znext = z^z;
myTable[] = z;
myTable[] = znext;
i = 3;
highFlag = False;
While[i <= max && ! highFlag,
If[Re[znext Log[z]] > 500,
(* flag this point as overflow *)
myTable[] = -10;
highFlag = True;
,
If[Re[znext Log[z]] < -500,
(* this is underflow, set to 0 to continue the cycle *)
myTable[[i]] = 0;
znext = 0;
,
znext = z^znext;
myTable[[i]] = znext;
];
];
i++;
];
myTable
], RuntimeAttributes -> {Listable}, Parallelization -> True]

tolerance = 0.001;
maxIterations = 500;
AbsoluteTiming[pointSummary = ParallelTable[
tetList = getTetrationList[x + I y, maxIterations];
(* check if point is overflow, if so, set to -1 for black color *)

If[tetList[] == -10,
theVal = -1;
,
theVal =
Length@DeleteDuplicates[tetList[[maxIterations - 100 ;;]],
Abs[#1 - #2] < tolerance &];
];
theVal, {y, -1, 0, 0.001}, {x, -1, 0, 0.001}];] • (1) It seems that underflow should also give an early termination. (2) You can create down values for the znext values, and terminate once a repeat is found (since it denotes the presence of a cycle). This would entail a change of how the cycle length is found. – Daniel Lichtblau Jan 7 at 15:01
• Thanks for that Daniel. I agree terminating after underflow would improve speed. Down values unfortunately are new to me but found some threads here about them. Also, seems some cycles aren't stabilizing after 500 iterations but need more so would take still longer. – Dominic Jan 7 at 16:07
• Here there are some suggestions for implementation of power towers – vi pa Jan 8 at 13:19
• For a demonstration on tetration see here – vi pa Jan 8 at 13:32
• Thanks guys. Very much appreciate any and all suggestions. I'm making some improvements. My goal is not so much to draw a nice picture but to study it's algebraic geometry and so an efficient and robust (error-free) algorithm to generate the data is my goal at present. – Dominic Jan 8 at 14:06

Sorry to hear you're stuck with the slow interpreted DeleteDuplicates. This new function is equivalent to DeleteDuplicates[list, Abs[#1-#2]<tolerance&], but compilation made it 12 times faster on my machine, and simply replacing the relevant part of your original code made it 4 times faster overall.

compiledDeleteDuplicates =
Compile[{{list, _Complex, 1}, {tolerance, _Real}},
Module[{unique = 0.*list, n = 0, i = 0, j = 0},
Do[
j = 1;
While[j <= n && Abs[unique[[j]] - list[[i]]] >= tolerance, j++];
If[j > n,
n++;
unique[[n]] = list[[i]]
];
, {i, 1, Length[list]}];
unique[[1 ;; n]]
]]

rnd = RandomComplex[{-1 - I, 1 + I}, 1000];
tolerance = 0.001;
AbsoluteTiming[result1 = DeleteDuplicates[rnd, Abs[#1 - #2] < tolerance &];]
AbsoluteTiming[result2 = compiledDeleteDuplicates[rnd, tolerance];]
result1 == result2
(*{0.782959, Null}
{0.06176, Null}
True*)

• Outstanding aooiii! Using the compiled version dropped execution time for the {{-1,0},{-1,0}} with $\Delta x=\Delta y=0.001$ run to 4 minutes. – Dominic Jan 11 at 11:13

I will update this with the plot (if it ever finishes!). Here is some code that I believe computes the value:

thePoints = ParallelTable[
{x, y, NestWhile[
{#[] + 1, Quiet@Chop[zStart^#[]]} &,
{1, zStart = x + I*y},
(#[] != 1) && (#[] < 500) &]},
{y, -1, 0, 0.001}, {x, -1, 0, 0.001}];


This doesn't take that long to run. Now, armed with the points, we can perform the limits by looking at any that have the cycle count equal to the limit, which is 500 in this case. Upper clip are the values that are large and lower the values that are small. Then, we use the third argument (the cycle count) to get a color. I'm new to doing this, so would appreciate any advice on a better method.

I'm using this:

ListPlot[Style[Take[#, 2], ColorData["NeonColors"][1/#[[-1]]]] & /@ Select[Flatten[thePoints, 1], #[[-1]] != 500 &]]


Here I didn't do the clipping and instead threw out the limit values (will come back to that in an edit). As I say, the plotting is taking a very long time (been going for over an hour and still not done).

EDIT 1: It appears that the reason the graph took so long is that each point was individually colored. By grouping the points by their color, the time drops to something reasonable, still about 2-4 minutes. Here is one such graph: To get this takes a few new steps. We flatten (because of the 2D table), then, if we are at the limit of the cycles, clip either up or down.

thePoints2 = Flatten[thePoints, 1];
thePointsClipped =
If[#[[3, 1]] == 500,
ReplacePart[#, {3, 1} ->
If[Abs[#[[3, 2]]] > 1000, 0, Infinity]], #] & /@ thePoints2;


Now, we get the pieces we want:

thePointsClipped2 = {#[], #[], #[[3, 1]]} & /@ thePointsClipped;
thePointsClipped3 = GatherBy[thePointsClipped2, Last];
finiteColors =
Select[Table[{i, ColorData[i][]}, {i, 1, 114}], #[[2, 2]] !=
Infinity &];
colorMappingRule =
MapIndexed[#2[] -> #1 &,
Flatten[MapIndexed[{#1, #2[]} &, #] & /@ (Table[#[], {#[[2,
2]]}] & /@ Take[finiteColors, 52]), 1]];
pointsToPlot =
MapIndexed[
Style[#1[[;; , {1, 2}]],
If[#1[[1, -1]] == 0, Black,
If[#1[[1, -1]] == Infinity, ColorData["Pastel"][.2],
ColorData @@ Lookup[colorMappingRule, First[#2]]]]] &,
thePointsClipped3];
ListPlot[pointsToPlot, AspectRatio -> 1, ImageSize -> Large]


There are other interesting plots by omitting the points that clipped either high or low using

ListPlot[Drop[pointsToPlot,2],AspectRatio -> 1, ImageSize -> Large]


EDIT 2: I implemented the code from OP and did similar work with the points to generate a graph. Here is what I see: This seems much closer to the graph shown in the original pointing.

• You don't define zstart in the first block of code. I assume a table has to be generated for some zstart in a grid of the complex plane for this? – Dominic Jan 10 at 14:05
• Nevermind about comment above, I understand now after looking at code. – Dominic Jan 10 at 14:10

(Not an actual answer, but a couple more suggestions):

1. (RuntimeAttributes -> {Listable}, Parallelization -> True) don't do anything if you call the compiled function as defined. They only make difference by parallelizing over arrays with more dimensions than in the definition. Try calling it like this for a speedup:
array = Table[x + I y, {y, -1, -0.001, 0.001}, {x, -1, -0.001, 0.001}];
getTetrationList[array, maxIterations]
1. DeleteDuplicates[list, test] is very slow compared to compiled functions. Not only it makes lots of Mathematica function calls, the number of these calls grows quadratically with the length of the list. After all, the main reason to use Compile is to get rid of millions of function calls. This is much faster - no extra calls and O(n*log(n)) complexity:
countUnique[x_, tol_] := Length[Union[Round[x/tol]]]

I admit it's not an ideal replacement, because e.g. countUnique[{0.49,0.51},1], countUnique[{0.499,0.501},1], countUnique[{0.4999,0.5001},1] and so on all return 2, but I think it's still applicable, and if not, maybe you could write a better compiled function for the task.

That's outstanding Mark! Give me some time to go over it. I've been trying to improve the speed of my code using the other suggestions above without success. It still takes about 20 minutes to generate the plot below with $${\Delta x,\Delta y}=0.001$$ but it's accurate. For example, looking at the scales in the plot, if I iterate $$-0.48-0.74i$$, then based on the color map below, it should have and 8-cycle.

Ran into problems trying to increase speed by using (or a variation of):

countUnique[x_, tol_] := Length[Union[Round[x/tol]]]


So I went back to DeleteDuplicates. This is the code I've come up with so far

(* have to run Colorize first *)
Colorize[{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}]
colorizeColors[n_] :=
Part[RGBColor[N@#1/256, N@#2/256, N@#3/256] & @@
ImageColorOperationsDumphashcolor[#1] & /@ Range[1, 101], n]
colorLegend =
Graphics[Table[{colorizeColors[x], Rectangle[{x, 0}, {x + 1, 1}],
White, Text[x, {x + 0.5, -.5}]}, {x, -1, 35}],
Background -> Black, ImageSize -> 400];

colorDataRules =
Flatten@Table[ToRules[n == colorizeColors[n]], {n, 1, 101}];
PrependTo[colorDataRules, -1 -> Black];

(* code to generate tetration cycle and manually check for overflow \
and underflow and return an array of iterations with array[] equal \
to the number of data points in the array *)

getTetrationList2 =
Compile[{{z, _Complex}, {max, _Integer}},
Module[{znext, zprev, i, highFlag, myTable = Table[1 + I, {max}]},
(* myTable[] stores the size of the data buffer in myTable *)
znext = z^z;
myTable[] = z;
myTable[] = znext;
i = 4;
highFlag = False;
While[i <= max - 1 && ! highFlag,
If[Re[znext Log[z]] > 500,
(* flag this point as overflow *)
myTable[] = -10;
highFlag = True;
,
If[Re[znext Log[z]] < -500,
(* this is underflow, set myTable[i]=0, myTable[i+1]=
1 and stop iteration *)
myTable[[i]] = 0;
myTable[[i + 1]] = 1;
myTable[] = i + 1;
highFlag = True;
,
(* neither over or underflow, continue iteration *)
zprev = znext;
znext = z^zprev;
If[Abs[zprev - znext] < 10^-5,
(* flag this as 1-cycle *)
highFlag = True;
myTable[] = -15;
,
myTable[[i]] = znext;
];
];
];
i++;
];
If[! highFlag,
myTable[] = max - 1;
];
myTable
], RuntimeAttributes -> {Listable}]

(* routine to acquire the iteration data *)
ymax = 0.03;
ymin = -1.0014;
xmax = 0.03;
xmin = -1.0014;
delta = 0.001;
totalX = Abs[xmax - xmin];
totalY = Abs[ymax - ymin];
array = Table[
x + I y, {y, ymax, ymin, -delta}, {x, xmin, xmax, delta}];
maxIterations = 500;

stats1 = AbsoluteTiming[
pointSummary = getTetrationList2[array, maxIterations];]

(* routine to acquire the cycle data *)
tolerance = 0.001;
rowLength = Length[pointSummary];
resultTable = Table[0, {i, 1, rowLength}, {j, 1, rowLength}];
AbsoluteTiming[For[i = 1, i <= rowLength, i++,
For[j = 1, j <= rowLength, j++,
totalLength = Re[pointSummary[[i, j, 1]]];
(* check if point is overflow, if so, set to -1 for black color *)

If[totalLength == -10,
theVal = -1;
,
If[totalLength == -15,
theVal = 1;
,
(* not under or overflow -- get length *)
If[totalLength > maxIterations,
Print[totalLength];
Abort[];
];
If[totalLength >= 200,
theVal =
Length@DeleteDuplicates[
pointSummary[[i, j, totalLength - 100 ;; totalLength]],
Abs[#1 - #2] < tolerance &];
,
theVal =
Length@DeleteDuplicates[
pointSummary[[i, j, 2 ;; totalLength]],
Abs[#1 - #2] < tolerance &];
];
If[theVal > 100,
theVal = 100;
];
];
];
resultTable[[i, j]] = theVal;
];
];]

myValues = pointSummary[[All, All, 2]];

(* use ArrayPlot to plot the data but the data ranges are now in \
terms of matrix size so for deltax=deltay=0.001 and region \
{{0,1},{0,1} the ranges will go from 1 to 1000.  Will need to \
manually re-scale axes for plot *)

myPlotA = ArrayPlot[resultTable, ColorRules -> colorDataRules];
myTetrationPlotLegend = Show[colorLegend];

myXVal[n_] := xmin + n/802;
totalTicks = Abs[xmin - xmax]/delta
ytickMax = Abs[ymin - ymax] /delta;
ydif = Abs[ymin - ymax]/4;
xdif = Abs[xmin - xmax]/4;
yDelta = Abs[ymin - ymax];
xDelta = Abs[xmin - xmax];
MtoRx[n_] := xmin + n/totalTicks (xDelta);
MtoRy[n_] := ymin + n/totalTicks (yDelta);

xtickMax = Abs[xmin - xmax]/ delta;
yaxis = Line@{{xtickMax/2, 0}, {xtickMax/2, ytickMax}};
xaxis = Line@{{0, ytickMax/2}, {xtickMax, ytickMax/2}};

myTetrationPlot =
Show[{myPlotA, Graphics@{xaxis, yaxis}}, Axes -> True,
FrameTicks -> {{{{0, ymin}, {ytickMax/4, ymin + ydif}, {ytickMax/2,
ymin + 2 ydif}, {3/4 ytickMax, ymin + 3 ydif}, {ytickMax,
ymax}}, {{0, ymin}, {ytickMax/4, ymin + ydif}, {ytickMax/2,
ymin + 2 ydif}, {3/4 ytickMax, ymin + 3 ydif}, {ytickMax,
ymax}}}, {{{0, xmin}, {xtickMax/4, xmin + xdif}, {ytickMax/2,
xmin + 2 xdif}, {3/4 ytickMax, xmin + 3 xdif}, {xtickMax,
xmax}}, {{0, xmin}, {xtickMax/4, xmin + xdif}, {ytickMax/2,
xmin + 2 xdif}, {3/4 ytickMax, xmin + 3 xdif}, {xtickMax,
xmax}}}}]
myTetrationPlotLegend = Show[colorLegend]


Just an update to the code base on using the compiled version of DeleteDuplicates written by a00iii above.

Here are the stats for range {{-1,0},{-1,0}} with $$\Delta x=\Delta y=0.001$$:

Time to acquire the tetration list: 2.5 minutes.

Time to identify the cycles using compiledDeleteDuplicates: 1.8 minutes.

That's excellent. I had planned to place a 100pt bounty for anyone who could get it below 10 minutes on my machine but can't seem to do so. If anyone knows how, I would be happy to award that to a00iii for his excellent solution to this problem.

Here's the updated code and a plot of the results:

(* routine to acquire the cycle data *)
tolerance = 0.001;
rowLength = Length[pointSummary];
resultTable = Table[0, {i, 1, rowLength}, {j, 1, rowLength}];
AbsoluteTiming[For[i = 1, i <= rowLength, i++,
For[j = 1, j <= rowLength, j++,
totalLength = Re[pointSummary[[i, j, 1]]];
(* check if point is overflow, if so, set to -1 for black color *)

If[totalLength == -10,
theVal = -1;
,
If[totalLength == -15,
theVal = 1;
,
(* not under or overflow -- get length *)
If[totalLength > maxIterations,
Print[totalLength];
Abort[];
];
If[totalLength >= 200,

theVal =
Length@compiledDeleteDuplicates[
pointSummary[[i, j, totalLength - 100 ;; totalLength]],
tolerance];
,

theVal =
Length@compiledDeleteDuplicates[
pointSummary[[i, j, 2 ;; totalLength]], tolerance];
];
If[theVal > 100,
theVal = 100;
];
];
];
resultTable[[i, j]] = theVal;
];
];]

myValues = pointSummary[[All, All, 2]]; 