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$

1 Answer 1

Reset to default
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$
3
  • $\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, 2020 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, 2020 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, 2020 at 18:01

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .