What approximation techniques exist for the square super-root function?
Some numerical stabs in the dark yielded the following for an iterative approach:

We're looking for the solution y = f(x) where y^y = x.

$$y \ln y = \ln x$$

$$y = g(x,y) = e^{\frac{\ln x}{y}}$$

The value for $y$ is a fixed point of the above equation, and empirically it seems to converge for some values of $x$, but for larger values of $x$ it oscillates or diverges.

Then I tried an approach similar to Newton's iterative square-root:

$$y = \frac{y_{previous} + y_{*}}{2} = \frac{y + e^{\frac{\ln x}{y}}}{2}$$

where y* is supposed to represent a nonconverging but optimistic answer that maintains accuracy if you happen to guess an accurate initial value (in square root y2 = x, it's y* = x/y).

This appears to converge, but very slowly at the low end of $x$ (near $x_{min} = (\frac{1}{e})^{\frac{1}{e}}$ )

It also looks like a good initial guess is $y_0 = \ln(x)+1$.

So I figured maybe there's a better-converging solution:

$y = (1-a) \times y + a \times g(x,y)$ for some value of $a$ that is a function of $x$.

Then I found something interesting.

If I get a converging answer $y$ from the above approach for $y^y = x$, and then compute $y_2 = g(x,y+ \epsilon) = e^{\frac{\ln(x)}{y+\epsilon}}$, it appears as though $y_2-y$ = approximately $\epsilon \times (-\ln(y))$.... e.g. if we had a guess $y_1 = y + \epsilon$ for some unknown $\epsilon$, and computed $y_2 = g(x,y_1)$, then $(y_2 - y) \approx \epsilon \times (-\ln(y)) = (y_1 - y) \times (-\ln(y))$. (Just to clarify, I have no analysis to verify this, but the numbers just popped out of some numerical evaluation I performed.)

Solve for the linear terms in $y$, and you get $y = \dfrac{y_2 + \ln(y) \times y_1}{1 + \ln(y)}$... use $\ln(y_1)$ in place of $\ln(y)$ and you get this iterative approximation:

$$\begin{align*}
    y[n+1] &amp;= \frac{g(x,y[n]) + \ln(y[n])\times y[n]}{1+\ln(y[n])}\\
           &amp;= \frac{e^{\frac{\ln(x)}{y[n]}} + \ln(y[n])\times y[n]}{1+\ln(y[n])}
\end{align*}$$

This appears to work very well, with the initial guess $y = 1+\ln(x)$, and appears to converge within 4 or 5 iterations.

(Someone could probably show that this is equivalent to Newton-Raphson in some way, but I think it's beyond my capability.)