In the frame of some physics work, I have to solve for $x$ the equation $$H(x)=(x!)^k$$ where $H(x)$ is the hyperfactorial function and $k$ a postive real number which can be very large.
I wrote it as $$f(x)=k \qquad \text{where} \qquad f(x)=\frac{\log (H(x))}{\log (x!)}$$ From a numerical point of view, the problem is not difficult since $f(x)$ is "almost" a linear function of $x$. Some information at two points $$\lim_{x\to 1} \, f(x)=\frac{\log (2 \pi )-3}{2 (\gamma -1)}\approx 1.37437$$ $$\lim_{x\to 1} \, f'(x)=\frac{42+12 (\gamma -3) \gamma -3 \pi ^2+\left(\pi ^2-6\right) \log (2 \pi )}{24 (\gamma -1)^2}\approx 0.634375$$ $$\lim_{x\to 1} \, f''(x)\approx -0.021330$$ while $$\lim_{x\to \infty} \, f(x)=\infty \qquad \lim_{x\to \infty} \, f'(x)=\frac 12\qquad \lim_{x\to \infty} \,f^{(n)}(x)=0\quad \forall n>1$$ I hope that this is sufficient to justify the quasi-linearity of $f(x)$.
Concerning the unsigned curvature, $\kappa(0)\approx 0.012843$ and $\kappa(x)$ decreases very fast $(\kappa(10)\approx 0.001240, \kappa(100)\approx 0.000062)$.
Expanding $f(x)$ for infinite values of x $$f(x)=\frac{(2 \log (x)-1)x}{4( \log (x)-1)}+\frac{\log (x) (2 \log (x)-3-2 \log (2 \pi ))+\log (2 \pi )}{8 (\log (x)-1)^2}+\cdots$$ For any $x$, the second term is very small compared to the first (the maximum value of their ratio is $0.00324$ at $x \sim 43$. So, ignoring it, the equation becomes $$\frac{(2 \log (x)-1)x}{4( \log (x)-1)}=k$$ and $$\frac{f(x)}{\frac{(2 \log (x)-1)x}{4( \log (x)-1)} }=1+\frac{2 \log ^2(x)-\log \left(\frac{4 }{3}\pi ^2\right) \log (x)+\log (2 \pi )}{2 x \left(2 \log ^2(x)-3 \log (x)+1\right)}+\cdots$$
If $x$ is really large, a very crude estimate could be $x_0=2k$. This is not too bad for Newton method (I give below the iterates for $k=1234$ $$\left( \begin{array}{cc} n & x_n \\ 0 & 2468.000000 \\ 1 & 2297.505131 \\ 2 & 2297.548546 \end{array} \right)$$ while is exact solution of the complete equation should be $2297.186319$.
To have a better approximation, letting $x=e^y$, we should have $$e^{-y}=\frac 1{2k}\frac{ y-\frac12}{ y-1}$$ the solution of which being given in terms of the generalized Lambert function (have a look at equation $(4)$); this is nice to know but not very practical.
Just for the art for art's sake, is there any way to generate some better estimates ?
Any idea or suggestion would really be welcome.
My post being already too long, I prefer to add an answer to it rather than to edit it.
I started with a different approach considering instead that I am looking for the zero of function $$g(x)=\log (H(x))-k \log (x!)$$ Using $x_0=2k$, the first iterate of Newton method is given by $$x_1=2k-\frac{2 \log (H(2 k))-2 k \log ((2k)!)}{2 \log ((2k)!)+4 k-2 k \psi(2 k+1)+1-\log (2 \pi )}$$ Now, using expansion for large values of $k$, I end with $$\color{blue}{x_1^*= 2 k-\frac{2 k+\log (2 k)-\log (2 \pi )}{2 \log (2 k)}+\frac{2 \log (2 k)+1}{4 \log ^2(2 k)}}$$
For the test example $(k=1234)$, this gives $x_1=2309.706772$, $x_1^*=2309.706724$ while the exact solution is $2297.186319$.
From a numerical point of view, it is better to consider $g(x)$ rather than $f(x)$ since, for any of the proposed guesses $x_*$ $f(x_*) \, f''(x_*) <0$ means that, by Darboux theorem, we should have one overshoot of the solution while $g(x_*) \, g''(x_*) <0$ guarantees no overshoot at all.
Update
Doing the same work with the approximating function mentioned in the post
$$h(x)=(2 \log (x)-1)x-{4k( \log (x)-1)}$$ which is much simpler to handle, using high order methods, I obtained things such $$x_2=2k \frac{\sum_{i=1}^9 a_i t^i}{\sum_{i=1}^9 b_i t^i} \qquad \text{where} \qquad t={\log(2k)}$$
The $a_i$'s correspond to the sequence $$\{667,-111648,-156240,1348032,803040,-4435200,564480,4515840,-3225600,645120\}$$
and the $b_i$'s to $$\{-4509,-78720,69552,1099392,-272160,-3548160,1693440,3225600,-2903040,645120\} $$
From those, truncated series could easily be obtained by long division; for example $$x_2=2k\left(1-\frac{1}{2 t}-\frac{1}{4 t^2}-\frac{3}{8 t^3}-\frac{1}{2 t^4}-\frac{77}{96t^5}+O\left(\frac{1}{t^6}\right)\right)$$
Using the terms of the present table, for the test case we should get $x_2=2297.54854638212$ while the solution of $h(x)=0$ is $2297.54854638189$ and the exact solution of $g(x)=0$ is $2297.186319$.