How fast do iterated exponentiation converge? Iterated exponentiation is defined by 
$$x \mapsto x^{x^{x^{\cdot^{\cdot^{\cdot}}}}}$$
or more conveniently, we denote by $^rx$ the expression $\underbrace{x^{x^{\cdot^{\cdot^{\cdot^{x}}}}}}_{r \text{ times}}$. Let us define the function 
$$\begin{array}{cccc}
f_x:& \mathbb{N}& \to & \mathbb{R}\\
& r & \mapsto & ^rx.
\end{array} $$
Euler proved that $\lim_{r \to \infty} f_x(r)$ exists for real numbers $x \in [e^{-e}, e^{1/e}]$ ([1]) (there is a convergence result in $\mathbb{C}$, but I am only interested in real numbers here).
Now we consider the Lambert $W$ function, which is defined to be the function satisfying 
$$W(z)e^{W(z)} = z.$$
Then we know from [1] that when $f_x$ converges, it converges to 
$$\lim_{r \to \infty} f_x(r) = \frac{W(-\ln(x))}{-\ln(x)}.$$
My question is: Given $\epsilon > 0$, how can we determine $r_0 \in \mathbb{N}$ such that $\left|f_x(r_0) - \frac{W(-\ln(x))}{-\ln(x)} \right| < \epsilon$? Is it possible to solve this algebraically, or is it only possible to do it numerically?
 A: Once you have picked an $x$ to begin with, you're iterating the function
$$ y \mapsto x^y $$
to find a fixed point for it. The general theory of function iteration now applies.
When $x$ is sufficiently small, the function will generally have two fixpoints, of which the smaller one is attractive and the larger is repelling. The slope $\frac{d}{dy}x^y$ is less than one at the attractive fixpoint, meaning that once you've gotten reasonably close to it, the distance between two successive iterations will decrease by about the same factor for each step.
During the iteration, you can compute and monitor these differences, and once they appear to stabilize you can use the theory of infinite geometric series to compute either an estimate of the fixpoint, or an estimate of the ratio between ${}^{r}x - {}^{r-1}x$ and the actual distance from ${}^rx$ to the fixpoint. This will allow you to compute a pretty good estimate of the limiting relation between $r$ and the remaining error.

Of course, if all you're interested in is finding the fixpoint, then this also suggests a quicker way to proceed than actually iterating $y\mapsto x^y$ one step at the time. You could either use the geometric-series estimate of the fixpoint alluded to above, or use Newton-Raphson iteration to find a root of $x^y-y$.
