# Can this power tower function be optimized to perform faster?

Take the following function defined in Mathematica:

Itr[x_, p_, n_] := x^Nest[Power[p, #] &, 1, n];


Evaluating this function even for small values of p results in very slow and processor-intensive evaluation. Is there an equivalent expression that will evaluate more efficiently and quickly?

For example, Itr[3, 1/2, 25] // N takes nearly 3 minutes to evaluate on a 2016 MacBook Pro.

Yes, use floating point numbers right from the start:

 Itr[3., 0.5, 50000000]


takes about one second on my machine.

You can also perform this in higher precision, but that will take longer; for example the following computation with 100-digit precision needs also about one second on my laptop:

 Itr[3.100, 0.5100, 500000];


Assuming that the nested powers converge quickly, one can also use FixedPoint in order to iterate as long as it is needed:

 f[x_, p_] := x^FixedPoint[Power[p, #] &, 1];


# Edit

Assuming that the sequence $$x_{n+1} = p^{x_n}$$ converges, towards an $$x \in \mathbb{R}$$ the limit point $$x$$ has to be a fixed point of

$$\varPhi_p(x) = p^x.$$

We can use Mathematica to solve this fixed point equation:

ClearAll[x,p];
F[p_] = Simplify[x /. Solve[Power[p, x] == x][[1]], {p > 0}]


-(ProductLog[-Log[p]]/Log[p])

So, we may use

f2[x_, p_] = x^F[p];


as even faster version. For example, we can compute 100000 limit points of Itr within a single second:

n = 100000;
x = RandomReal[{1, 40}, n];
p = RandomReal[{0.1, .9}, n];
data = f2[x, p]; // AbsoluteTiming // First


1.00969