# Numerical instability of an extended tetration

For bases $$a\in(1,e^{1/e})$$, $${}^na=a^{({}^{n-1}a)}=a^{a^{a^{.^{.^{.^a}}}}}$$ converges to a value denoted as $${}^\infty a$$. By observing the convergence rate of this sequence, we can derive the limit:

$$\lim_{n\to\infty}\frac{{}^\infty a-{}^{n+x}a}{{}^\infty a-{}^na}=[\ln({}^\infty a)]^x$$

By supposing we seek a continuous version of tetration that satisfies this, and rearranging so that $${}^xa$$ is solved for, we derive:

$${}^xa=\lim_{n\to\infty}\log_a^{\circ n}({}^\infty a-({}^\infty a-{}^na)[\ln({}^\infty a)]^x)\tag1$$

where $$\log^{\circ n}$$ is the logarithm applied $$n$$ times. As an example, with $$n=10$$, I obtained the following plot: which looks really nice. Then looking at $$n=15$$, I get this: which raises concern. For $$a$$ close to $$1$$ and $$n=10$$, I get It would seem to work well for small $$n$$ and large $$a$$, but then for larger $$n$$ or smaller $$a$$, it becomes unstable. As far as I can tell, this issue is due to the amount of numerical precision required while evaluating $$(1)$$, especially when the base is closer to $$1$$.

So the first question is whether this is due to numerical precision, or if it's simply because $$(1)$$ does not converge.

If it's the former, then is there any way to circumvent this without brute forcing with more precision? And how should I pick the values of $$n$$ for a given base $$a$$ (and $$x$$)?

If it's the latter, then does it converge anywhere?

Code for computing $$(1)$$, showing the following for $$a=\sqrt2$$ and $$x=1.5$$:

n   1.4142135623730^^n
--------------------------
0   1.42291711861386
1   1.4657586018199498
2   1.4910645646490854
3   1.5069501895748705
4   1.5172760309843982
5   1.5241342747726574
6   1.528753204049527
7   1.5318927292918296
8   1.5340399138955585
9   1.5355145848360043
10  1.5365302824374432
...
45  1.538805432574356
46  1.5388054445894592
47  1.5388054519338499
48  1.5388054652284342
49  1.5388054823911386
50  1.538805506512146
...
90  1.7233534923554696
91  1.755592017472159
92  2.0000000000000004
93  2.000000000000001
94  2.000000000000001


Showing the apparent value of $$^{1.5}\sqrt2\simeq1.5388$$ followed by divergence.

• Wait... I hope I didn't miss my 7 o'clock question hat! xd – Simply Beautiful Art Dec 26 '19 at 0:09
• look at this paper arxiv.org/abs/1908.05559 and also the 1981 article in the Monthly quoted in the references – Math101 Dec 26 '19 at 0:32
• Giving a quick read, the mentioned paper is just an elaborate discussion on the convergence of ${}^na$ as $n\to\infty$. The question at hand discusses the convergence of $$\log_a^{\circ n}({}^\infty a-({}^\infty a-{}^na)[\ln({}^\infty a)]^x)$$ as $n\to\infty$, for $a\in(1,e^{1/e})$. – Simply Beautiful Art Dec 26 '19 at 1:03
• @SimplyBeautifulArt It does look like a round-off error. You calculate ${}^\infty a$ by bisection with $n=100$. Then around $k=100$ you get instability, recovering basically $2={}^\infty(\sqrt{2})$. I guess it is because ${}^\infty a-{}^n a$ becomes zero. You may try Newton method instead of bisection to see if there will be any difference. – A.Γ. Dec 26 '19 at 1:39
• +1 for the way to extend tetration that actually seems somewhat natural and makes sense. – URL Dec 27 '19 at 8:21

Conclusion of what follows below is that the formula/implementation is subject to catastrophic cancellation in the part $${}^\infty a-{}^na$$.

I have implemented your code verbatim in MATLAB. The result with double precision is identical (subject to $$\pm 1$$ shift somehow in the counter $$k$$) for $$k=25,50,100$$.

>> test_double(25)

inf_a - n_a = 4.5947e-05

res =

1.538798995072006

>> test_double(50)

inf_a - n_a = 4.8182e-09

res =

1.538805596316793

>> test_double(100)

inf_a - n_a = 1.3323e-15  % Close to machine epsilon for 64 bits (!!!)

res =

2.000000000000001


Now increasing the precision to quadruple (digits = 32)

>> test_vpa(25,32)

inf_a - n_a = 0.00004594710217926553685144416637951

res =

1.538798995052591060169673385115

>> test_vpa(50,32)

inf_a - n_a = 0.0000000048181840417120357729198170741061

res =

1.5388054167659345033199452184372

>> test_vpa(100,32)

inf_a - n_a = 0.00000000000000005298542667702117495721500428964

res =

1.5388054174394144484626385875352


One can even go for more $$k$$, but for $$k=200$$ it becomes complex (maybe time to improve the bisection(?))

>> test_vpa(150,32)

inf_a - n_a = 0.00000000000000000000000058267893924374831695038060802142

res =

1.5388052522455767035302332127693

>> test_vpa(200,32)

inf_a - n_a = -0.0000000000000000000000000000002180085870192506738614539008334

res =

4.5484210435493524026130653949673 + 0.27008458814289299475281495540935i


Just for fun: $$k=300$$, 50 digits precision, bisection has 200 iterations:

>> test_vpa(300,50)

inf_a - n_a =

ans =

0.00000000000000000000000000000000000000000000000077490532526797719918187015369561788340262562842819

res =

1.5388054183673573750604847502205024901223348381524

• Thanks! Good to know it's most likely instability and not lack of convergence. – Simply Beautiful Art Dec 26 '19 at 16:26

Another example suggesting the problem is numerical instability of the double-precision based computations.

Using Pari/GP with internal precision of 200 decimal digits I get the following list:

    n  a^^1.5 (based on n iterations)
---------------------------------------
1  1.4657586018199495028
2  1.4910645646490847958
3  1.5069501895748684802
4  1.5172760309843957556
.....
27  1.5388009658025622487
28  1.5388023318106759423
29  1.5388032786497405906
30  1.5388039349458743480
...
83  1.5388054174394108432
84  1.5388054174394125092
85  1.5388054174394136640
86  1.5388054174394144644
87  1.5388054174394150192
88  1.5388054174394154038
89  1.5388054174394156704
90  1.5388054174394158551
91  1.5388054174394159832
92  1.5388054174394160720
93  1.5388054174394161335
94  1.5388054174394161762
95  1.5388054174394162057
96  1.5388054174394162262
97  1.5388054174394162404
98  1.5388054174394162503
99  1.5388054174394162571
100  1.5388054174394162618
....
200  1.538805417439416272500504006016378087864
201  1.538805417439416272500504006016378484370
202  1.538805417439416272500504006016378759207
203  1.538805417439416272500504006016378949709
204  1.538805417439416272500504006016379081755
205  1.538805417439416272500504006016379173283
206  1.538805417439416272500504006016379236725
207  1.538805417439416272500504006016379280699
208  1.538805417439416272500504006016379311180
209  1.538805417439416272500504006016379332308


The Pari/GP-routine was

default(realprecision,200) \\ this is my standard internal precision request
\\ performs tetration to (i)nteger "heights" by base b
itet(x,h)=for(k=1,h,x=b^x);for(k=1,-h,x=log(x)/log(b));return(x)

b=sqrt(2)
\\ one could find the fixpoint (or "b^^infinity") by simple binary search
b_inf = solve(t=1,exp(1), b^t - t )
\\ but more efficient and more general is the following formula:
t0=b_inf=exp(-LambertW(-log(b)))  \\ fixpoint: compute b^^infinity

\\ tetration of an initial value z0 with base b
\\           to fractional heights x
\\               with n iterations,
\\                      where n should go to infinity
{tetx(z0,x,max_n=10)=my(w);
w= t0-(t0-itet(z0,max_n))*log(t0)^x;
itet(w,-max_n)}

\\ compute a list of estimates based on n iterations, n=1..100
list = Mat(vectorv(100,n, [n, tetx(1,1.5,n)]))

\\ compute additional approximations, n=200.. 209
list1 = Mat(vectorv(10,n, [199+n, tetx(1, 1.5, 199+n )]))


We compare consecutive terms. Let $$q=\ln({}^\infty a)$$. Note that:

\begin{align}\log_a({}^\infty a-({}^\infty a-{}^{n+1}a)q^x)&=\log_a({}^\infty a)+\log_a\left(\frac{{}^\infty a-{}^{n+1}a}{{}^\infty a}q^x\right)\\&={}^\infty a+\frac{{}^\infty a-{}^{n+1}a}{\ln({}^\infty a)}q^x+\mathcal O(q^{2n})\\&={}^\infty a+({}^\infty a-{}^na)q^x+\mathcal O(q^{2n})\end{align}

Now we need to log both sides $$n$$ times. Note that

$$\log_a(\alpha+\epsilon)=\log_a(\alpha)+\frac{\epsilon}{\alpha\ln(a)}+\mathcal O(\epsilon^2)$$

We can also verify through a similar process that the limit is monotone increasing, and hence $$\alpha\ge\alpha'$$, which implies the error between consecutive terms to be $$\mathcal O(q^2/\alpha'\ln(a))^n\simeq\mathcal O(q({}^\infty a/\alpha'))^n$$. It is known that $$q<1$$, so all that remains to be seen is that $$\alpha'$$ is sufficiently close to $${}^\infty a$$. From this, a rough estimate of $$\alpha'$$ can be deduced as $${}^xa$$, with the error given from above, and hence for all $$a$$ and sufficiently large $$x$$, the limit exists, and furthermore, converges monotonically and uniformly to an analytic function.

By disregarding the smallest $$N$$ values that $$\alpha$$ takes, we can get a better bound of $$\alpha'\simeq{}^{x+N}a$$, and hence we can get the error's ratio to go below $$1$$ eventually, and hence this holds everywhere.

The above works for any $$x\in\mathbb C$$ since $$[\ln({}^\infty a)]^{x+N}\to0$$ as $$N\to\infty$$ even for complex $$x$$.

• This whole ansatz seem to me to reflect the Schroeder's method, reduced to the first nonzero term of the resulting power-series (theta-series?), where the following terms are neglected due to the limit process. If you're interested in this aspect I could put something more explicite in an additional answer. – Gottfried Helms Dec 29 '19 at 12:30
• Sure why not.${}$ – Simply Beautiful Art Dec 29 '19 at 13:17