Section A: Introduction to fixed points of iterated exponentials

$$ \newcommand{\bint}{\displaystyle{\int\hspace{-10.4pt}\Large\mathit{8}}} \newcommand{\res}{\displaystyle{\text{Res}}} \newcommand{Log}{\text{Log}} \newcommand{Arg}{\text{Arg}} \newcommand{pLog}{\text{pLog}} $$

1. Background:

This page reviews algorithms and Mathematica code used to numerically compute and index fixed points of $z^{z^w}$ described in the paper A systematic approach to computing and indexing the fixed points of an iterated exponential. Readers are encouraged to review the paper first.

$1$-cycle and $2$-cycle expressions for constant $z$ are: \begin{equation} \begin{aligned} T_1(w;z)&=w-z^w\equiv w-e^{w\Log(z)}, \\ T_2(w;z)&=w-z^{z^w}\equiv w-e^{z^w \Log(z)}\equiv w-e^{\Log(z)e^{w\Log(z)}}, \end{aligned} \label{equation:eqn0} \end{equation} where $\Log(z)=\ln|z|+i\theta$ with $\theta=\Arg(z)$ and $-\pi<\Arg(z)\leq \pi$. The notation $\log(z)$ is the multivalued logarithmic function base $e$. $\Log(z)$ is the principal value of $\log(z)$. Likewise, $\arg(z)$ is the multivalued argument function, and $\Arg(z)$ is the principal value of $\arg(z)$.

It's relatively easy to generate a plot showing some of the zeros of (1): Let $\Log(z)=a+bi$ and $w=x+iy$. In the case of $T_2$, first expand the $1$-cycle expressions \begin{equation} \begin{aligned} z^w&=e^{(a+bi)(x+iy)}=e^{(ax-by)+i(ay+bx)}\\ &=e^{(ax-by)}\bigg[\cos(ay+bx)+i\sin(ay+bx)\bigg]. \end{aligned} \label{equation:eqn21} \end{equation} Now let \begin{equation} \begin{aligned} e^{(ax-by)}\cos(ay+bx)&=c \\ e^{(ax-by)}\sin(ay+bx)&=d. \label{equation:eqn22} \end{aligned} \end{equation} Then $$ \begin{aligned} z^{z^w}&=e^{s(a+bi)}=e^{(a+bi)(c+di)}\\ &=\text{exp}\big\{(ac-bd)\left[\cos(ad+bc)+i\sin(ad+bc)\right]\big\}=x+iy, \end{aligned} $$ so that in order for $w=z^{z^w}$: $$ \begin{aligned} y&=e^{(ac-bd)}\sin(ad+bc)\\ x&=e^{(ac-bd)}\cos(ad+bc).\\ \label{equation:eqn22b} \end{aligned} $$ Plugging in $c$ and $d$ above, define \begin{equation} \begin{split} \text{imagF2B}(x,y)&=\text{exp}\bigg\{e^{ax-by}\big[a\cos(ay+bx)-b\sin(ay+bx)\big]\bigg\}\\ &\hspace{20pt}\times\sin\bigg\{e^{ax-by}\big[a\sin(ay+bx)+b\cos(ay+bx)\big]\bigg\} \end{split} \end{equation} \begin{equation} \begin{split} \text{realF2B}(x,y)&=\text{exp}\bigg\{e^{ax-by}\big[a\cos(ay+bx)-b\sin(ay+bx)\big]\bigg\}\\ &\hspace{20pt}\times\cos\bigg\{e^{ax-by}\big[a\sin(ay+bx)+b\cos(ay+bx)\big]\bigg\}. \end{split} \end{equation}

Next use Mathematica's $\texttt{ContourPlot}$ function to generate a plot of $\text{realF2B}(x,y)=x$ and $\text{imagF2B}(x,y)=y$:

Mathematica code
zSample = 2;
reim = Im;
a = Log[Abs[zSample]];
b = Arg[zSample];
realF1[x_, y_] := Exp[a x - b y] Cos[a y + b x];
imagF1[x_, y_] := Exp[a x - b y] Sin[a y + b x];
realF2[x_, y_] :=Exp[a realF1[x, y] - b imagF1[x, y]] Cos[a imagF1[x, y] + b realF1[x, y]];
imagF2[x_, y_] := Exp[a realF1[x, y] - b imagF1[x, y]] Sin[a imagF1[x, y] + b realF1[x, y]];
realContourP = ContourPlot[realF2[x, y] == x, {x, -2, 6}, {y, -15, 15},ContourStyle -> Red, PlotPoints -> 25];
imagContourP =ContourPlot[imagF2[x, y] == y, {x, -2, 6}, {y, -15, 15},ContourStyle -> Blue,PlotPoints -> 25];
Show[{realContourP, imagContourP}]


Figure 1: Contour plot for z=2
Figure 1 shows a small section of the contour plot with the real contours in red and the imaginary contours in blue. The zeros are segregated on foliated "branch bulbs" at the intersections of the red and blue contours. Branches are delineated by the purple and cyan line segments left of the bulbs with the branch number between the segments. Zeros are indexed as $\{n,m,p\}$ with $m$ being the branch number, and $n$ the root number or "leaf number" corresponding to a contour intersection. The third index, $p$, is zero except when multiple roots are assigned to an $\{n,m\}$ pair. In general, leaf numbers at the top of the branch bulb have positive $n$ and those at the bottom, negative $n$. Root numbers with $n=0$ are usually at the apex of a branch bulb and those near the origin will sometimes have multiple root assignments. See paper for more information about this indexing scheme.

A simple way of finding a zero would be to use the Mathematica "get coordinate" feature to choose a value of $w$ in the plot close to a root, then use $\texttt{FindRoot}$. Let's say in a Mathematica notebook with this plot, we place the cursor near one of the intersections and obtain $(0.8,1.7)$. Then simply code: $$ \texttt{FindRoot[w==zSample^(zSample^w),{w,0.8+1.7I}]} $$

Mathematica returns $0.824679 + 1.56743 I$ which is a very close approximation to this root. We can increase the accuracy by adding a working precision: $$ \texttt{FindRoot[w==zSample^(zSample^w),{w,0.8+1.7I},WorkingPrecsion->25]} $$ Which returns w= $0.8246785461420742223140646 + 1.567432123849647861058574 I$. But this of course is not a practical method of computing roots for several reasons:
• Iterating $T_2$ in its exponential form $z^{z^w}$ frequently produces underflow and overflow.
• Mathematica interprets $z^{z^w}$ as principal value. However, the underlying complex geometry is infinitely-valued due to the implicit complex $\log$ terms.
• As can be seen in the plot, the intersections grow increasingly closer together; the root density along the branch bulb increases as the bulb extends outward towards the right of the diagram. For example, on the same branch bulb, root one-trillion and one-trillion+1 differ only after the 11'th decimal place: $$ \begin{aligned} 1000000000000&: 43.0433996410655766643648805958518803348+2.2661800709127327385371991312567672208 I \\ 1000000000001 &: 43.0433996410670193594057688254300272329+2.2661800709127327385371999877483600194 I \end{aligned} $$
And notice one particularly interesting phenomenon of iterated exponentials: There are about one-trillion roots with positive $n$ on this branch lobe in the interval $(0,43.05)$. How many are in the interval $(0,100)$? Here are two just beyond this: $$ \small \begin{aligned} 1000000000000000000000000000000 &:\\ &\hspace{-100pt}102.83810534903808684840035995956853+ 2.2661800709135969048138414728572695209249858117859292708894581544i\\ 1000000000000000000000000000001 &:\\ &\hspace{-100pt}102.83810534903808684840035995956998+ 2.2661800709135969048138414728572695209249858117859292708894592150i. \end{aligned} $$ If a zero is removed from the index numbers, the real part of the root drop below $100$. Let's define a function for this: $\Phi(m,\rho)$ gives the number of roots of $T_2(w,2)$ on branch m with real part less than $\rho$. For example $\Phi(0,10)=226$.

Two "bracket" roots are associated with $\Phi(0,\rho)$: $\{\phi_{\rho},0\}$ and $\{\phi_{\rho}+1,0\}$ with $\phi_{\rho}>0$ so that $\text{Re}\{\phi_{\rho},0\}\lt \rho$ and $\text{Re}\{\phi_{\rho}+1,0\}\geq\rho$.

2. A challenge for the interested reader:

Compute $\Phi(0,10^2)$. Here are the first and last digits for several others: $$ \begin{aligned} \Phi(0,10^3)&=2364129446. . .1647386392 \\ \Phi(0,10^4)&=4401819481. . .1892131900 \\ \Phi(0,10^5)&=2204154263. . .7983989378\\ \Phi(0,10^6)&=2184437229. . .6948769810 \\ \Phi(0,10^7)&=1996711872. . .0967813102\\ \Phi(0,10^8)&=8129684803. . .6940652584\\ \Phi(0,10^9)&=1017786728. . .7141566978\\ \end{aligned} $$

3. Newton Iteration:

Experimentally, a fixed-point iteration of the exponential form of $T_2$ was found to be problematic: the exponential terms often cause overflow and underflow. But if the iterated exponential is converted to its logarithmic form, this effect was virtually eliminated in all test cases studied thus far: \begin{equation} \begin{aligned} w&=e^{\Log(z)e^{w\Log(z)}}\\ (\Log(w)+2n\pi i)&=\Log(z)e^{w\Log(z)};\quad n\in\mathbb{Z} \\ \frac{1}{\Log(z)}\bigg[\Log(w)+2n\pi i\bigg]&=e^{w\Log(z)}\\ w&=\frac{1}{\Log(z)}\bigg\{\Log\bigg[\frac{1}{\Log(z)}\left(\Log(w)+2n\pi i\right)\bigg]+2m\pi i\bigg\};\quad (n,m)\in\mathbb{Z}^2\\ w&=\pLog(w,n,m;z) \end{aligned} \label{equation:eqn1} \end{equation} so that $T_2$ in it's logarithmic form is \begin{equation} T_2(w,n,m;z)=w-\pLog(w,n,m;z);\quad (n,m)\in\mathbb{Z}^2. \end{equation} We will work with the hypothesis that all the roots of $T_2$ can be enumerated by the set of triplets $\{n,m,p\}$ with $(n,m)$ enumerated by $\mathbb{Z}^2$ discussed in the reference and that these zeros are organized primarily in the form of foliated branch leaves in a contour plot including degenerate leaves near the origin. The plot above shows three branches and their corresponding (red and blue) leaf lobes. And that the branching extends to infinity and the number of leaves on each branch is infinite as well.

A Newton iteration of $T_2$ takes the form \begin{equation} w_{k+1}=w_{k}-\frac{w_{k}-\pLog\left(w_{k},n,m\right)}{1-\frac{1}{w_{k}\Log(z)\left(\Log(w_{k})+2n\pi i\right)}}, \label{eqn:eqn500} \end{equation} where \begin{equation} \frac{d}{dw}\pLog(w,n,m)=\frac{1}{w\Log(z)\left(\Log(w)+2n\pi i\right)}. \end{equation}

The immediate question regarding (8) is, what are the basins of attraction? Unfortunately, determining the basins of attraction directly for a fixed-point iteration of a complex function is usually not possible: One usually iterates a region of the complex plane identifying what sub-regions converge or diverge. The reference paper cited above includes several diagrams of basin plots showing convergence regions for roots near the origin. All test cases, including those not included in the paper showed the root basins for all but the roots near the origin encompassed a large area surrounding multiple branches. And using the techniques described in the paper, we can always compute analytically, an iteration seed close to the branch head (bulb) with arbitrary precision. And since sequential roots $n,n+1,n+2,\cdots$ on the same branch are on sequential leaf bulbs which grow closer together, for a sequential set of roots, we can always update the seed with the previously-calculated root to obtain seeds which grow increasingly close to the respective root.

The task now becomes, designing a robust algorithm to implement (8) using arbitrary precision arithmetic. Before we do this, we review the computations to compute the branching metrics. This is covered in Section B: Computing the branching parameters