1
$\begingroup$

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)$$

$\endgroup$
2
$\begingroup$

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 + O\left( \frac{1}{x} \right) $$ as $x\to +\infty$. This is because the difference between this approximation and the one you gave is $O(1/x)$. Thus $$ \frac{4}{e}\log \left( \frac{H(x)e^{1/8}}{A} \right) = \frac{1}{e}\left( x^2 + x + \frac{1}{6} \right)\log \left( \frac{1}{e}\left( x^2 + x + \frac{1}{6} \right)\right) + O\left( \frac{1}{x} \right), $$ and hence $$ x^2 + x + \frac{1}{6} \approx \frac{4\log \left( \frac{H(x)e^{1/8}}{A} \right)}{W\left( \frac{4}{e}\log \left( \frac{H(x)e^{1/8} }{A} \right) \right)}. $$ Solving for $x$ yields $$ x \approx - \frac{1}{2} + \sqrt {\frac{4\log \left( \frac{H(x)e^{1/8}}{A} \right)}{W\left( \frac{4}{e}\log \left( \frac{H(x)e^{1/8} }{A} \right) \right)} + \frac{1}{12}} . $$

$\endgroup$
  • $\begingroup$ It is a real pleasure to have an answer from you. We spoke with @robjohn about your work math.stackexchange.com/questions/430167/…. I shall immediately try and be back to you. $\endgroup$ – Claude Leibovici Jan 13 '20 at 8:35
  • $\begingroup$ Super-extra-mega-magnificent is to weak as a qualifier for this solution ! You made my day. Thank you very much ! Cheers :-) $\endgroup$ – Claude Leibovici Jan 13 '20 at 8:53
  • $\begingroup$ You are correct about the error term. It seems to me that $x^2 +x +\frac{1}{6}$ is a natural variable as $$\log H(x) \sim \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 - \frac{1}{480\left( x^2 + x + \frac{1}{6} \right)}+ \frac{31}{90720\left( x^2 + x + \frac{1}{6} \right)^2 } + \cdots$$. $\endgroup$ – Gary Jan 13 '20 at 18:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.