Another solution which does not involve the LambertW function involves generating the following series inverse, which I will explain in steps.
$$f=\sqrt{2(\exp(x)-x-1)};\;\;\;\;\text{xfixed}=f^{-1}(x)$$
$$\text{xfixed}= x - \frac{x^2}{6}+ \frac{x^3}{36} - \frac{x^4}{270} + \frac{x^5}{4320} + \frac{x^6}{17010} - \frac{139x^7}{5443200} ...$$
In pari-gp, this series would be "serreverse(sqrt((exp(x)-x-1)*2))". I don't have a closed form for the series. By plugging both square roots into this series, you can get either fixed point L for the base=b.
$$L_b = \frac{\text{xfixed}\left(\pm\sqrt{-2(\ln(\ln(b))+1)}\right)-\ln(\ln(b))}{\ln(b)}$$
One interesting thing about the square root radical is that if $b=\exp(1/e)$, then the radical is zero, which makes sense since $b=\exp(1/e)$ has only one fixed point. Using a 16 term series for xfixed, I got about 19 decimal digits of precision for the two fixed points of $b=\sqrt{2}$ So, why does this work? We must first show that iterating $$z \mapsto b^z$$ is exactly congruent with iterating
$$y \mapsto \exp(y)+k;\;\;k=\ln(\ln(b));\;\;\;y=z\ln(b)+\ln(\ln(b))$$
So then our fixed point problem for $z \mapsto b^z$ becomes equivalent to finding the fixed point for $y \mapsto \exp(y)+k$. If we have the fixed point for y, then we can get the fixed point for $z \mapsto b^z$ by using this equation based on the congruency.
$$z = \frac{y-\ln(ln(b))}{\ln(b)}$$
Next, we do some algebra to get the fixed point for
$$y = \exp(y)+k;\;\;\;k=\ln\ln(b))$$
$$y = k + 1 + y + \frac{y^2}{2} + \frac{y^3}{6} + \frac{y^4}{24} + \frac{y^5}{5!} ...$$
$$ -(k+1) = \frac{y^2}{2} + \frac{y^3}{6} + \frac{y^4}{4!} ...$$
$$ -2(k+1) = y^2 + \frac{2y^3}{6} + \frac{2y^4}{24} + \frac{2y^5}{5!} ...$$
Now we take the square root of both sides.
$$ \sqrt{-2(k+1)} = f(y) = y\cdot \sqrt{1 + \frac{2y}{6} + \frac{2y^2}{24} + \frac{2y^3}{5!} ...}$$
$$ \sqrt{-2(k+1)} = f(y) =
y + \frac{y^2}{6} + \frac{y^3}{36} + \frac{y^4}{270} + \frac{y^5}{2592} + \frac{17y^6}{544320} + \frac{11y^7}{5443200}+...$$
The square root term on the right has a nice defined formal Taylor series in y; which leads exactly to the $f$ in the 2nd line of my post. So then we can get k from y. We need the reverse, so we take the inverse of Taylor series for f.
$$ y = f^{-1}\left(\pm\sqrt{-2(k+1)}\right)=f^{-1}\left(\pm\sqrt{-2(\ln(ln(b))+1)}\right)$$
Now we have the fixed point $y=\exp(y)+k$, and the final step is to generate the fixed point $z=b^z$ from y using $z = \frac{y-\ln(ln(b))}{\ln(b)}$
