As a small part in a statistical thermodynamics project, I need to compute the inverse of the hyperfactorial function.
So,as I wrote it, I need to find the zero of function
$$f(x)=\log (H(x))-k$$ for which
$$f'(x)=\log (\Gamma (x+1))+x+\frac{1}{2} (1-\log (2 \pi ))\qquad \text{and} \qquad f''(x)=\psi (x+1)+1$$
Since $k$ is large, for the estimate of the starting guess, I used the asymptotics
$$\log (H(x))=\frac{1}{4} x^2 (2 \log (x)-1)+\frac{1}{12} (6 x+1) \log (x)+\log (A)+\sum_{n=1}^\infty a_n x^{-2n}$$ where the $a_n$ form the sequence
$$\left\{\frac{1}{720},-\frac{1}{5040},\frac{1}{10080},-\frac{1}{9504},\frac{691}{360
3600},-\frac{1}{1872},\frac{3617}{1713600},-\frac{43867}{3907008},\frac{174611}{
2257200}\right\}$$
The estimate was made using the first term only
$$\frac{1}{4} x^2 (2 \log (x)-1)=k \implies x_0=\sqrt{\frac{4 k}{W\left(\frac{4 k}{e}\right)}}$$ The good point is that $f(x_0) >0$ and $f''(x_0)>0$ which means that, by Darboux theorem, Newton method would converge without any overshoot of the solution.
For sure, using floating point arithmetics, I cannot compute $H(x)$ and I just used the expansion in which the series has been truncated to the very firts terms but the derivative was exact. However, no approximation for the derivatives.
Using the above, I computed the first iterate of Newton method $(x_1)$ as well as the the first iterate of Halley method $(x_2)$.
Using $k=2^p$, here are some results
$$\left(
\begin{array}{ccccc}
p & x_0 & x_1 & x_2 & \text{exact} \\
1 & 2.7733509 & 2.3214362 & 2.2551702 & 2.2442276 \\
2 & 3.3553862 & 2.8968477 & 2.8436979 & 2.8372181 \\
3 & 4.1586005 & 3.6933378 & 3.6514727 & 3.6477083 \\
4 & 5.2543815 & 4.7827661 & 4.7502650 & 4.7481083 \\
5 & 6.7413690 & 6.2640778 & 6.2391502 & 6.2379290 \\
6 & 8.7556108 & 8.2734629 & 8.2545399 & 8.2538554 \\
7 & 11.484401 & 10.998235 & 10.983995 & 10.983615 \\
8 & 15.185387 & 14.695981 & 14.685344 & 14.685135 \\
9 & 20.213017 & 19.721051 & 19.713156 & 19.713041 \\
10 & 27.055187 & 26.561232 & 26.555402 & 26.555340 \\
11 & 36.384023 & 35.888542 & 35.884255 & 35.884222 \\
12 & 49.126276 & 48.629637 & 48.626495 & 48.626477 \\
13 & 66.560960 & 66.063447 & 66.061152 & 66.061143 \\
14 & 90.454838 & 89.956673 & 89.955000 & 89.954995 \\
15 & 123.25055 & 122.75190 & 122.75068 & 122.75068 \\
16 & 168.32793 & 167.82892 & 167.82804 & 167.82804 \\
17 & 230.36727 & 229.86799 & 229.86735 & 229.86735 \\
18 & 315.85443 & 315.35496 & 315.35449 & 315.35449 \\
19 & 433.78360 & 433.28399 & 433.28365 & 433.28365 \\
20 & 596.63558 & 596.13586 & 596.13561 & 596.13561 \\
21 & 821.73989 & 821.24009 & 821.23991 & 821.23991 \\
22 & 1133.1726 & 1132.6727 & 1132.6726 & 1132.6726 \\
23 & 1564.4008 & 1563.9009 & 1563.9009 & 1563.9009
\end{array}
\right)$$
Just remember that $H(1500) \sim 2.894 \times 10^{3331194}$.
My question is : could be proposed a simpler approximation of the inverse of the hyperfactorial in the same spirit as for the inverse of the factorial function (see here) that is to say without any iteration ?
Edit
In the same spirit of what he already did for the inverse factorial, @Gary proposed here a agnificent solution to the problem.
Written as
$$x \sim \sqrt{\frac{e t}{W(t)}+\frac{1}{12}}-\frac{1}{2} \qquad \text{with} \qquad t=\frac{8(k-\log (A))+1}{2 e}$$
Just to give an idea, I produce below the "bad" results (again for $k=2^p$)
$$\left(
\begin{array}{ccc}
p & \text{approximation} & \text{exact} \\
1 & \color{red}{2.244}1282 & 2.2442276 \\
2 & \color{red}{2.837}1718 & 2.8372181 \\
3 & \color{red}{3.647}6879 & 3.6477083 \\
4 & \color{red}{4.748}0997 & 4.7481083 \\
5 & \color{red}{6.23792}53 & 6.2379288 \\
6 & \color{red}{8.25385}39 & 8.2538553 \\
7 & \color{red}{10.983615} & 10.983615
\end{array}
\right)$$
In fact, @Gary was too modest since the difference between the two series is $$\frac 1{480x^2}\left(1-\frac 1 x+O\left(\frac{1}{x^2}\right) \right)$$
Update
If we consider the new expansion added by @Gary in comments, the difference between the two series is
$$\frac {103}{725760 x^6}\left(1-\frac 3 x+O\left(\frac{1}{x^2}\right) \right)$$
Best Answer
You may check that $$ \log H(x) = \frac{1}{4}\left( x^2 + x + \frac{1}{6} \right)\log \left( x^2 + x + \frac{1}{6} \right) - \frac{1}{4}\left( x^2 + x + \frac{1}{6} \right) - \frac{1}{8} + \log A + \mathcal{O}\!\left( \frac{1}{x} \right) $$ as $x\to +\infty$. This is because the difference between this approximation and the one you gave is $\mathcal{O}(1/x)$. Thus $$ \frac{4}{\mathrm{e}}\log \left( \frac{H(x)e^{1/8}}{A} \right) = \frac{1}{\mathrm{e}}\left( x^2 + x + \frac{1}{6} \right)\log \left( \frac{1}{\mathrm{e}}\left( x^2 + x + \frac{1}{6} \right)\right) + \mathcal{O}\!\left( \frac{1}{x} \right), $$ and hence $$ x^2 + x + \frac{1}{6} \approx \frac{4\log \left( \frac{H(x)\mathrm{e}^{1/8}}{A} \right)}{W\!\left( \frac{4}{\mathrm{e}}\log \left( \frac{H(x)\mathrm{e}^{1/8} }{A} \right) \right)}. $$ Solving for $x$ yields $$ x \approx - \frac{1}{2} + \sqrt {\frac{4\log \left( \frac{H(x)\mathrm{e}^{1/8}}{A} \right)}{W\!\left( \frac{4}{\mathrm{e}}\log \left( \frac{H(x)\mathrm{e}^{1/8} }{A} \right) \right)} + \frac{1}{12}} \;. $$