14

In the frame on a research work, I needed to find quickly the smallest $p$ such that $$\large p!!\geq A$$ $A$ being a very large number.

In a first step, I tried to solve the equality using the superb approximation @robjohn produced eight years ago; it write

$$ n!=a^n \,A \qquad \implies \qquad n\sim e\,a \exp\Bigg(W\left(\log \left(\frac{A^2}{2 \pi a}\right)\right) \Bigg)-\frac 12 \tag 1$$ where $W(.)$ is the principal branch of Lambert function.

For the first case

$$(2n)!!=2^n\, \Gamma (n+1)$$ it is simple since $a=\frac 12$, leading to $$n_1 \sim e\,\exp\Bigg(W\left(\log \left(\frac{A^2}{\pi}\right)\right) \Bigg)-1 \quad \text{with}\qquad A=(2n)!!$$

This gives very decent results $$\left( \begin{array}{cc} p & \text{approximation}\\ 2 & 1.95022 \\ 4 & 3.97945 \\ 6 & 5.98782 \\ 8 & 7.99159 \\ 10 & 9.99369 \\ 20 & 19.9974 \\ 30 & 29.9984 \\ 40 & 39.9989 \\ 50 & 49.9992 \\ 60 & 59.9993 \\ 70 & 69.9994 \\ 80 & 79.9995 \\ 90 & 89.9996 \\ 100 & 99.9996 \\ \end{array} \right)$$

For the second case $$(2n+1)!!=\frac{2^{n+1}}{\sqrt{\pi }} \Gamma \left(n+\frac{3}{2}\right)$$

it is a bit more delicate since $a=\frac 12 \left(\frac{\sqrt{\pi }}{2}\right)^{\frac{1}{n}} $. Being short of idea, I approximated it by $a\sim\frac 12 \left(\frac{\sqrt{\pi}}{2}\right)^{\frac{1}{n_1}} $ and used $(1)$ with $A=(2n+1)!!$

Disclaimer the following has been modified after @njuffa's answer

$$\large \color{blue}{p\sim e\,\exp\Bigg[W\left(\frac{1}{e}\log \left(\frac{A^2}{\alpha}\right)\right) \Bigg]-1}$$

In the first version, I used $\color{red}{\alpha=e}$ while, from @njuffa's interesting approach and solution $\color{red}{\alpha=\sqrt{2\pi}}$ (which differ by $8$%)

Some results for the equation $$\color{red}{\huge p!!=10^k}$$

$$\left( \begin{array}{ccc} k & \text{estimate} & \text{solution} & \text{estimate} \\ & \text{with }\alpha=e & &\text{with }\alpha=\sqrt{2\pi} \\ 1 & 4.33968 & 4.34257 & 4.38794 \\ 2 & 6.79424 & 6.95182 & 6.83367 \\ 3 & 8.90847 & 9.04794 & 8.94379 \\ 4 & 10.8393 & 10.9684 & 10.8721 \\ 5 & 12.6494 & 12.7450 & 12.6804 \\ 6 & 14.3714 & 14.3729 & 14.4010 \\ 7 & 16.0253 & 15.9778 & 16.0538 \\ 8 & 17.6241 & 17.6253 & 17.6518 \\ 9 & 19.1773 & 19.2586 & 19.2043 \\ 10 & 20.6916 & 20.7763 & 20.7179 \\ 20 & 34.4678 & 34.4898 & 34.4905 \\ 30 & 46.8288 & 46.9065 & 46.8497 \\ 40 & 58.3999 & 58.4041 & 58.4198 \\ 50 & 69.4384 & 69.4640 & 69.4574 \\ 60 & 80.0819 & 80.0500 & 80.1003 \\ 70 & 90.4152 & 90.4213 & 90.4331 \\ 80 & 100.495 & 100.515 & 100.513 \\ 90 & 110.362 & 110.359 & 110.380 \\ 100 & 120.047 & 120.017 & 120.064 \\ 200 & 210.335 & 210.328 & 210.350 \\ 300 & 293.630 & 293.628 & 293.644 \\ 400 & 372.898 & 372.950 & 372.912 \\ 500 & 449.397 & 449.419 & 449.410 \\ 600 & 523.811 & 523.795 & 523.824 \\ 700 & 596.570 & 596.593 & 596.583 \\ 800 & 667.964 & 667.942 & 667.977 \\ 900 & 738.203 & 738.186 & 738.215 \\ 1000 & 807.443 & 807.459 & 807.455 \\ \end{array} \right)$$

So, at the time being, we have approximations for the inverse of the factorial, the hyperfactorial and the double factorial functions

My questions

  • Could we improve this one ?
  • Could we generalize the problem to multiple factorials ?

Any idea or suggestion would be very welcome.

  • The first case can be further improved for larger n by adding the asymptotic correction term of +1/(25*p). p_approx is plenty good enough for this purpose. The corrected p value is then at least an order of magnitude more accurate for p>40. – Martin Brown Oct 15 '25 at 08:55
  • @MartinBrown. Are you suggesting $n_2=n_1+\frac{1}{25n_1}$ ? – Claude Leibovici Oct 16 '25 at 10:45
  • Yes although I don't know how to typeset it here. I think a similar although more ad hoc correction term also applies to @Njuffa's answer for `n>10000 (approx 0.0599/n1^1.097) with the appropriate sign. – Martin Brown Oct 16 '25 at 12:54
  • @MartinBrown. Could you explain where $\frac 1{25 p}$ is precisely coming from ? I am very interested. Thanks & cheers :-) – Claude Leibovici Oct 16 '25 at 12:59
  • I think it is somehow related to using the higher order Stirling approximation for n! although why it is 1/(25n) and not 1/(24n) is a mystery to me. I should caution that I am a physicist rather than a mathematician but I see patterns in numbers. – Martin Brown Oct 16 '25 at 13:20
  • @njuffa thank you for the clarification of how it arises. I'm sure your formula is right (that hybrid of competing high order terms would explain the slight overshoot). A similar although much more empirical and grubby tweak also works for your rather fine answer below. Unclear if it holds good for larger n though. – Martin Brown Oct 16 '25 at 21:31

2 Answers2

9

Using $(5.11.7)$, we obtain $$ \Gamma \!\left( n + \tfrac{3}{2} \right) = \Gamma \!\left( n + 1 + \tfrac{1}{2} \right) \sim \left( \frac{n + 1}{\rm e} \right)^{n + 1} \sqrt {2\pi } $$ as $n\to+\infty$. Thus, the equation $A=(2n+1)!!$ may be approximated via $$ A \approx \left( \frac{2(n + 1)}{\rm e} \right)^{n + 1} \sqrt 2 \Longleftrightarrow \frac{2}{\rm e}\log \left( \frac{A}{\sqrt 2 } \right) \approx \frac{2(n + 1)}{\rm e} \log \left( \frac{2(n + 1)}{\rm e} \right). $$ Using the principal branch of the Lambert $W$-function, we can solve this for $n$ as follows: $$ n \approx \frac{\rm e}{2}\exp \left( W\!\left( \frac{2}{\rm e}\log \left( \frac{A}{\sqrt 2 } \right) \right)\right) - 1. $$

FD_bfa
  • 4,927
  • 4
  • 17
  • 55
Gary
  • 37,398
3

When I did this a while ago, I started with Stirling's original approximation to $n!$ (what we now call Stirling's approximation is apparently due to De Moivre, as I learned at History of Science and Mathematics Stackexchange): $$n! \approx \left(\frac{n+\frac{1}{2}}{e}\right)^{n+\frac{1}{2}}\sqrt{2 \pi}$$ This is readily invertible via the Lambert $W$ function . This yielded: $$\exp\left(W\left(\frac{1}{e}\ln\left(\frac{n!}{\sqrt{2 \pi}}\right)\right)+1\right)-\frac{1}{2} \approx n$$

For $k$ even, I then approximated $$k!! \left(\frac{k!!}{\sqrt{\frac{\pi}{2} k }}\right)\approx k!$$

With this, the approximation for the inverse of the double factorial $k!!$ for even $k$ became

$$\exp\left(W\left(\frac{1}{e}\ln\left(\frac{(k!!)^{2}}{\pi}\right)\right)+1\right)-1 \approx k$$

Likewise, for $k$ odd, I approximated $$k!! \left(\frac{k!!}{\sqrt{\frac{2}{\pi} k }}\right)\approx k!$$

With this, the approximation for the inverse of the double factorial $k!!$ for odd $k$ became

$$\exp\left(W\left(\frac{1}{e}\ln\left(\frac{(k!!)^{2}}{2}\right)\right)+1\right)-1 \approx k$$

The approximation formulas for odd and even $k$ differ only in the divisor of the innermost term. To generate an approximate inverse of the double factorial $k!!$ for any $k$, I therefore simply "split the difference" between the divisors of $\pi$ and $2$, which causes error to oscillate evenly between overestimation and underestimation:

$$\exp \left(W \left(\frac{1}{e} \ln \left(\frac{(k!!)^{2}}{\sqrt{2 \pi}}\right) \right) + 1\right) - 1 \approx k$$

The following table shows the numerical behavior of this blended approximation: $$ \begin{array}{|r|c|l|} \hline k & \mathrm{approx.}\ \mathrm{inverse}\ \mathrm{of}\ k!! & \mathrm{relative}\ \mathrm{error}\\ \hline 9 & 8.89448 &-1.17249 \cdot 10^{-2} \\ \hline 10 & 10.0877& +8.77103 \cdot 10^{-3} \\ \hline 99 & 98.9506 & -4.98935 \cdot 10^{-4} \\ \hline 100 & 100.049 & +4.85642 \cdot 10^{-4} \\ \hline 999 & 998.967 & -3.27436 \cdot 10^{-5} \\ \hline 1000 & 1000.03 & +3.26577 \cdot 10^{-5} \\ \hline 9999 & 9998.98 & -2.45192 \cdot 10^{-6} \\ \hline 10000 & 10000.0 & +2.45129 \cdot 10^{-6} \\ \hline 99999 & 99999.0 & -1.96123 \cdot 10^{-7} \\ \hline 100000 & 100000.& +1.96118 \cdot 10^{-7} \\ \hline \end{array} $$

njuffa
  • 2,333