# Calculating the nth super-root when n is greater than 2?

Tetration (literally "4th operator iteration") is iterated exponentiation, much like how exponents are iterated multiplying. For example, $$2$$^^$$3$$ is the same as $$2^{2^2}$$, which is $$2^4$$ or 16. This applies to any two positive real numbers, and possibly even all real numbers.

There are two inverses of tetration, iterated logarithms "super-logarithms" or "tetralogarithms" and iterated roots "super-roots". For instance, the value of $$\sqrt{256}_s$$ is $$4$$ because $$4$$^^$$2 = 4^4 = 256$$.

It is very difficult to calculate with or approximate either of these inverse operations to tetration. Super-logarithms use some complicated recursive algorithms (see this wiki article), but super-roots are even harder; the example I showed just turned out to be easy-to-find because $$256$$ is a perfect super-square (that is, $$4^4$$). I could only manage to find a formula for super square roots on Wikipedia using the Lambert W function. I originally did not know what this Lambert W function was, but soon after I found out it's just the inverse function of $$xe^x$$. This expression stated that $$ssrt(x) = ln(x)/W(ln(x))$$.

However, I cannot find a formula or even approximation for higher super-roots at all. I tried using linear approximation along with a Euclidean distance version of linear approximation on $$^3\sqrt{2}_s$$ and both were very innacurate. I could eventually calculate $$^3\sqrt{2}_s$$ to 42 digits as $$1.4766843373578699470892355853738898365517$$ using this high-precision calculator and binary guess-and-check methods, but I still can't seem to find a generic formula or approximation for these super-roots which have orders higher than $$2$$, so I would like to know if anyone knows of one.

• I've explored this a bit in 2017 and have made a *.pdf of the results. If you find this helpful I might extract something of it into a new answer. See go.helms-net.de/math/tetdocs/Wexzal_Superroot.pdf (the term "wexzal" is just taken from some old article of two mathematicians - Fantini and another one - who made a very nice study about the second superroot, you can google it, as old as is it is written in ascii) – Gottfried Helms Aug 6 '19 at 7:10
• There has also been a discussion of this in the tetration-forum. One member even proposed he has found a solution for the $n=3$-case. See math.eretrandre.org/tetrationforum/… as a start – Gottfried Helms Aug 8 '19 at 5:05
• A further link into the tetration forum with a short index to various discussions math.eretrandre.org/tetrationforum/… – Gottfried Helms Aug 8 '19 at 5:25
• A recursive approximation-procedure has been provided by the new user @Tomasz (who has not enough reputation to post a comment). He says, there is "a recursive rule to calculate rank $3$ super root from $a$. The proof can be shown in a few lines". Iterate "$x(n+1) = \exp(W(W(x(n) \cdot \ln (a))))$." (I assume, $W()$ is the Lambert-W() function). Unfortunately he gives no initial definition for $x(0)$, but I've tried this simply with $x(0)=a$ and it seems to work: Let $s=\lim_{n \to \infty} x(n)$ then $s^{s^s} = a$. (This recursion was given without proof) – Gottfried Helms Oct 22 '20 at 17:42
• I should add to my previous comment, that this seems to be simply the procedure that @OscarLanzi describes in his answer, and tested with his initial values I reproduce perfectly his iterative numbers. – Gottfried Helms Oct 23 '20 at 10:39

One thing we can do is render the third superroot as an iteration of second superroots.

Let $$x^{x^x}=a$$. Then

$$(x^x)^{(x^x)}=a^x$$

and we take two square superroots to get a fixed point iteration:

$$\color{blue}{x=\sqrt{\sqrt{a^x}_s}_s}$$

If we put $$a>1$$ and $$x=1$$ on the right side, we get

$$x=\sqrt{\sqrt{a}_s}_s$$

as our next iteration, and this will lie between $$1$$ and the true third superroot. This suggests that the scheme will be convergent.

Let's try this on the third superroot of $$16$$. Put $$x=1$$ on the right side of the fixed point relation above. Then we get

$$x=\sqrt{\sqrt{16}_s}_s=1.769\text{ (4 digits)}$$

Next time around we get to the same accuracy

$$x=\sqrt{\sqrt{134.9}_s}_s=1.958\text{ (4 digits)}$$

We seem to be converging to $$2$$ fairly rapidly.

For the third superroot of 2 the corresponding iterations are $$1.000, 1.380, 1.459, 1.473, 1.476, \color{blue}{1.477},...$$.

For smaller $$a$$ values the performance deteriorates, an outcome connected with square superroots being no longer single valued when the argument drops below 1. If $$a=0.9$$ and an initial value $$x=1$$ is tried, the iteration converges to about $$0.890$$ if we take square superroots $$\ge 1/e$$; but this does not work for $$a=0.8$$ due to one of the square superroots being undefined over the reals.

Just for sake of simplicity, it is possible to compute $$\sqrt[n]a_s$$ very easily when $$a>1$$ or $$n$$ is odd and $$a>0$$ using bisection, which doesn't involve any hard calculations (so this is something doable with a calculator that only has exponentiation, and perhaps some paper).

As long as you have $$x\le\sqrt[n]a_s\le y$$, we can iteratively consider $$[(x+y)/2]\widehat~\widehat~n$$, and replace either $$x$$ or $$y$$ with $$(x+y)/2$$.

To find initial values to work with, note that $$a\ge\sqrt[n]a_s$$. From there, repeatedly divide by 2 until you find a value where $$x\le\sqrt[n]a_s$$, and set $$y$$ to be $$2x$$. Then apply the above.

Here is a simple program that implements this. It takes 53 steps to compute the result

$$\sqrt2_s\approx1.4766843373578697$$

It is also notable that the above program may fail on larger inputs, since computing things like 256^^3 may not be reasonably doable. For a human, however, if the radicand is small enough, it is easy to find initial values manually.

Since some programming languages give $$a^b=\infty$$ if $$a$$ and $$b$$ are too large, and $$\infty$$ is considered greater than anything, there is no need to worry about larger inputs. This allows for the computation of things such as $$\sqrt{256}_s$$ in only 60 steps. See here.

Improvement:

One can employ a faster root-finding algorithm that works well against functions with large second derivatives to get faster practical use. The idea is similar to the above, where we define a new function

$$f(x)=x\widehat~\widehat~n-a$$

which we want to find the root of. We use the initial boundaries

$$(x_l,x_u)=(\min(1,a),\max(1,a))$$

or $$x_u=2$$ if $$f(2)>0$$, and then iteratively define the bisection and regula falsi midpoints:

$$x_\mathrm B:=\frac{x_u+x_l}2$$

$$x_\mathrm{RF}:=\frac{x_uf(x_l)-x_lf(x_u)}{f(x_l)-f(x_u)}=\frac{\frac{x_u}{f(x_u)}-\frac{x_l}{f(x_l)}}{\frac1{f(x_u)}-\frac1{f(x_l)}}\stackrel{f(x_u)\to\infty}\longrightarrow x_l$$

where $$x_\mathrm{RF}:=x_l$$ is used if $$f(x_u)$$ is infinite in double precision. We then define our new endpoint as

$$x_r:=\begin{cases}\frac12(x_\mathrm B+x_\mathrm{RF}),&\text{steps}\bmod3\ne0\text{ or }x_u-x_l>0.5\text{ or }|f(x_\mathrm{RF})|>0.5\\x_\mathrm{RF},&\text{else}\end{cases}$$

$$x_l:=x_r\text{ if }f(x_r)<0$$

$$x_u:=x_r\text{ if }f(x_r)>0$$

which converges roughly twice as fast as bisection at worst and slightly faster than linear convergence asymptotically. See here for an implementation of the above, which computes $$\sqrt2_s$$ and $$\sqrt{256}_s$$ in 9 and 21 iterations respectively.

• $x=a^{1/a}$ is already bounded, and while $\exp(W(\log(a)))$ is a nicer upper bound, it involves functions which are not so basic, which likely are just as hard to compute as what I've given. – Simply Beautiful Art Aug 8 '19 at 18:52
• That isn't really an issue since $\ln(a)$ works as an upper bound for $a>e$, for example. But note the algorithms in my answer already $\mathcal O(\ln(a))$. – Simply Beautiful Art Aug 8 '19 at 19:10
• I made a significant error in my comments. For $a>e$ all superroots lay above $e^{1/e}$ although then $a^{1/a}<e^{1/e}$. So my thoughts about the lower bounds for the search intervals were incorrect. I'll thus delete my earlier comments, rethinking. (The series solutions in my own answer are correct anyway) – Gottfried Helms Aug 8 '19 at 22:30

There is also an (individual) powerseries solution (Puisieux-series) for each $$n$$ separately. Unfortunately that series have a little radius of convergence (if nonzero at all), but might be summable using Euler-summation. I'll give an example for $$n=3$$. (More examples are in my small treatize on my webspace)

Let's define our basic function $$v = f_3(u) = u \cdot \exp( u \cdot \exp(u))$$ This has a power series that Pari/GP would easily display to many, say $$64$$, coefficients. It is Lagrange-invertible, say $$u = g_3(v) = g_{3,1} v + g_{3,2}/2! v^2 + g_{3,3}/3! v^3 + ... \qquad \qquad \\\ \qquad \qquad\qquad \qquad = v - 2/2! v^2 + 3/3! v^3 + 20/4! v^4 - 295/5! v^5 + 1554/6! v^6 + O(v^7)$$ Let's call the $$n$$'th superroot-function $$r_n(y)$$ such that $$x = r_n(y) \to y = \;^n x$$ then with $$n=3$$ we have asymptotically (or with small range of convergence) $$r_3(y) = \exp(g_3(\log(y)))$$

By visual impression using the first 64 coefficients it seems, as if $$|y| \lt \exp(1/2)$$ might give convergence. (Update: checked with $$\small y=\exp(\exp(-1))\approx 1.444$$ and it seems to be computable without Euler-summation at all)
Using Euler-summation the radius of convergence seems to be extensible - but I do not have the exact estimate for the growth of the coefficients. Let's see the summation of $$g_3(v)$$ for $$y=\exp(1/2)$$ such that $$v = \log(y) =1/2$$

partial sums of series g_3(1/2)
index  direct par-   Eulersummati- x=exp(g_3(v))
tial sums     on order 0.5
0              0             0  1.000000000
1   0.5000000000  0.2500000000  1.284025417
2   0.2500000000  0.3125000000  1.366837941
3   0.3125000000  0.3203125000  1.377558184
4   0.3645833333  0.3196614583  1.376661628
5   0.2877604167  0.3198649089  1.376941739
6   0.3214843750  0.3205749512  1.377919773
7   0.3465603299  0.3209987217  1.378503819
8   0.2958108569  0.3211143463  1.378663217
9   0.3245950850  0.3211192669  1.378670001
10   0.3403287830  0.3211180191  1.378668280
11   0.2981245348  0.3211251293  1.378678083
12   0.3269177171  0.3211317745  1.378687245
13   0.3373635102  0.3211344667  1.378690956
14   0.2982079698  0.3211348582  1.378691496
...     ...             ...           ...
55   0.1639008657  0.3211350702  1.378691789
56   0.2430865079  0.3211350702  1.378691789
57   0.5951181182  0.3211350702  1.378691789
58  0.09819587958  0.3211350702  1.378691789
59   0.2552736284  0.3211350702  1.378691789
60   0.6589999027  0.3211350702  1.378691789
61  0.01025715717  0.3211350702  1.378691789
62   0.2806542398  0.3211350702  1.378691789
63   0.7342007399  0.3211350702  1.378691789
64  -0.1064953930  0.3211350702  1.378691789
65   0.3257696112  0.3211350702  1.378691789
66   0.8210482830  0.3211350702  1.378691789


giving the result

\\ v = 1/2
\\ y = exp(v) = 1.648721271
\\ x = 1.378691789...  from the last entry in the above (exp() of Euler-summation)
x^x^x     \\ check result of exp(Eulersummation)
%368 = 1.648721271
x^x^x- y  \\ check error
%369 = -1.699715302 E-22


Using Euler-summation of higher order and still using $$64$$ coefficients I could approximate $$x$$ for $$y=16$$ thus $$v=\log(16) \approx 2.772588722$$ to a couple of digits

 partial sums of series g_3(v)
index    direct par-   Eulersummati- x=exp(g_3(v))
tial sums     on order 5.8
0                0             0  1.000000000
1      2.772588722  0.4138192123  1.512583646
2     -4.914659500  0.5946280226  1.812356664
3      5.742129363  0.6381957710  1.893062278
4      54.98695039  0.6360439550  1.888993136
5     -347.7931741  0.6382424852  1.893150713
6      632.6698779  0.6531183317  1.921523444
7      4675.314968  0.6695911376  1.953438470
8     -40693.25817  0.6786063818  1.971128816
9      101996.7689  0.6804959536  1.974856927
10      534495.9511  0.6805394727  1.974942873
11     -5898689.882  0.6823339559  1.978490056
12      18438758.44  0.6856034795  1.984969362
13      67398852.83  0.6883322388  1.990393263
14     -950277060.1  0.6894759540  1.992671008
...      ...            ...           ...
55  -2.664574273E40  0.6931453682  1.999996375
56   9.487882296E39  0.6931457660  1.999997171
57   9.002509460E41  0.6931461912  1.999998021
58  -6.072183170E42  0.6931464204  1.999998480
59   6.149361593E42  0.6931464569  1.999998553
60   1.803356223E44  0.6931464662  1.999998571
61  -1.371747030E45  0.6931465681  1.999998775
62   2.215487571E45  0.6931467349  1.999999109
63   3.558078850E46  0.6931468656  1.999999370
64  -3.073667769E47  0.6931469109  1.999999461
65   6.704426035E47  0.6931469096  1.999999458
66   6.882982587E48  0.6931469259  1.999999491


We see (and conclude for the general trend), that direct evaluation of the series $$g_3(v)$$ leads to unbounded partial sums, but which might be handled successfully by Eulersummation of appropriate order. For negative $$v$$ this might look different because Euler-summation can only be successful if the summands-to-be-partial-summed are (at least roughly) alternating in signs.

The same procedure can be done for higher $$n$$, getting simply different power series, and for $$n \to \infty$$ we get a well known powerseries with simple set of coefficients, see my small treatize which I have linked to.

update 2 Just for fun I tried $$y=256$$, $$v \approx 5.545177444$$. Using Eulersummation of order $$11.37$$ I got the following, with now $$128$$ terms of the power series needed:

partial sums of series g_3(v)
index      direct par-   Eulersummati- x=exp(g_3(v))
tial sums     on order 11.37
0                 0             0  1.000000000
1       5.545177444  0.4482762688  1.565611165
2      -25.20381545  0.6593619374  1.933558210
3       60.05049546  0.7137176680  2.041567036
4       847.9676320  0.7103553682  2.034714203
5      -12040.99635  0.7136655574  2.041460652
6       50708.63898  0.7376632728  2.091043603
7       568167.2105  0.7661483392  2.151463567
8      -11046187.51  0.7827563692  2.187493503
9       62011106.34  0.7863107520  2.195282527
10       504890269.0  0.7863264545  2.195316998
11   -1.267027432E10  0.7906504964  2.204830194
12    8.701591401E10  0.7990721357  2.223476886
13    4.880970073E11  0.8065436624  2.240151869
14   -1.618550515E13  0.8098085843  2.247477743
...        ...              ...          ...
115   6.351017290E120  0.8292221643  2.291535609
116   1.111485492E123  0.8292222056  2.291535703
117  -1.527759513E124  0.8292223349  2.291536000
118   3.224076332E124  0.8292224524  2.291536269
119   1.850204173E126  0.8292224818  2.291536336
120  -2.865500371E127  0.8292224486  2.291536260
121   9.519211594E127  0.8292224379  2.291536236
122   3.024663081E129  0.8292224940  2.291536364
123  -5.316849919E130  0.8292225794  2.291536560
124   2.371536870E131  0.8292226279  2.291536671
125   4.834907255E132  0.8292226204  2.291536654
126  -9.764577483E133  0.8292225984  2.291536603
127   5.423935796E134  0.8292226105  2.291536631
128   7.506467361E135  0.8292226592  2.291536743
129  -1.775490700E137  0.8292227055  2.291536849
130   1.177132710E138  0.8292227171  2.291536875


Meaning the result

\\ y = 256
\\ v = log(y) ~ 5.545
\\ x = 2.291536875...  from the last entry in the above (exp() of Euler-summation)
x^x^x     \\ check result of exp(Eulersummation)
%465 = 255.9990213


Using $$256$$ series coefficients and increased internal precision I got the result correct to even the eleventh digit.
So this as an initial solution and then Newton-interpolation (perhaps somehow following the idea of Oscar Lanzi) one might arrive at arbitrary accuracy - again: I did not check the lower bound for $$y$$...

Just a long comment, no answer

@SimplyBeautifulMath's simple answer triggered me to have another look at the searching/solver ansatz.
It is not so elegant as SBM's procedere, but fast enough and also easy. Here is the Pari/GP-code

itet(b,h=1)=my(x=1);for(k=1,h,x=b^x);x  \\ integer tetration to height h

{ sroot(m,maxorder=1)=my(logm,vsroots,ub);
logm=log(m) ;              \\ make the log a constant
vsroots=vectorv(maxorder); \\ column-vector of superroots of order 1 to (maxorder)
vsroots=m;              \\ define sr of order 1 to be m
vsroots= exp(LambertW(logm));  \\ use LambertW for direct estimate
\\            for sr of order 2
for(o=3,maxorder,          \\ for all requested higher orders...
ub=vsroots[o-1];   \\ use previous sr as upper bound for current sr
vsroots[o]=solve(sr=1,ub, log(sr) * itet(sr,o-1)-logm) \\ solve
\\ for identity of log (sr^^o) = log(m)
);
return(vsroots);}
\\ do the computation to get a table of 20 x 10 entries.
tmp= Mat(vector(10,e,sroot(10^e,20)~))


The idea is here to use a "good" estimate for the upper-bound for $$sr$$ of order $$o$$ to avoid numerical overflow in the search-routine - just take the $$sr$$ of smaller order. Using $$\log(sr \uparrow \uparrow o) = \log(m)$$ reduces again the problem of numerical overflow.

This made the following table for 1 to 20'th superroot of powers of 10.

 order       m -->
----+---------------------------------------------------------------------------------- ------------------
1         10        100       1000      10000     100000    1000000   10000000  100000000  1000000000  10000000000
2  2.5061841  3.5972850  4.5555357  5.4385827  6.2709196  7.0657967  7.8313895  8.5731845   9.2950869    10.000000
3  1.9235840  2.2127958  2.3849098  2.5072341  2.6019716  2.6791912  2.7443082  2.8005682   2.8500690    2.8942439
4  1.7343125  1.8662152  1.9319788  1.9743344  2.0050183  2.0288114  2.0480991  2.0642312   2.0780410    2.0900767
5  1.6440791  1.7188182  1.7522674  1.7726883  1.7869853  1.7978034  1.8064095  1.8134996   1.8194934    1.8246618
6  1.5924155  1.6399925  1.6597509  1.6714095  1.6794032  1.6853641  1.6900542  1.6938845   1.6970995    1.6998551
7  1.5594683  1.5920497  1.6048495  1.6122253  1.6172116  1.6208940  1.6237705  1.6261064   1.6280579    1.6297241
8  1.5369070  1.5603758  1.5692069  1.5742077  1.5775542  1.5800085  1.5819159  1.5834587   1.5847435    1.5858375
9  1.5206559  1.5382042  1.5445837  1.5481481  1.5505152  1.5522422  1.5535793  1.5546577   1.5555535    1.5563148
10  1.5084979  1.5220043  1.5267778  1.5294168  1.5311588  1.5324247  1.5334019  1.5341882   1.5348403    1.5353935
11  1.4991297  1.5097690  1.5134418  1.5154549  1.5167774  1.5177354  1.5184732  1.5190657   1.5195564    1.5199723
12  1.4917381  1.5002805  1.5031713  1.5047446  1.5057741  1.5065179  1.5070896  1.5075482   1.5079274    1.5082485
13  1.4857915  1.4927612  1.4950799  1.4963344  1.4971526  1.4977424  1.4981951  1.4985577   1.4988573    1.4991108
14  1.4809287  1.4866938  1.4885837  1.4896010  1.4902626  1.4907387  1.4911036  1.4913956   1.4916367    1.4918405
15  1.4768968  1.4817223  1.4832840  1.4841209  1.4846640  1.4850541  1.4853528  1.4855916   1.4857886    1.4859551
16  1.4735134  1.4775948  1.4789007  1.4795980  1.4800494  1.4803733  1.4806210  1.4808190   1.4809821    1.4811200
17  1.4706445  1.4741283  1.4752319  1.4758191  1.4761987  1.4764707  1.4766785  1.4768445   1.4769812    1.4770967
18  1.4681892  1.4711875  1.4721288  1.4726282  1.4729505  1.4731812  1.4733573  1.4734979   1.4736137    1.4737114
19  1.4660708  1.4686702  1.4694798  1.4699082  1.4701842  1.4703816  1.4705323  1.4706524   1.4707513    1.4708348
20  1.4642296  1.4664982  1.4671997  1.4675700  1.4678083  1.4679786  1.4681085  1.4682120   1.4682972    1.4683691


As already said, I find of course the procedere as shown by SBM much more elegant and also more interesting. This comment is just to have a table of values visible.