Using Gosper’s factorial approximation to determine the maximal $n$ in $n! = 2^p$

approximationcalculusfactoriallogarithmsrecreational-mathematics

I'm trying to use the Gosper's factorial approximation to get a rough value of the maximal integer $n$ in $$n!\approx n^ne^{-n}\sqrt{\left(2n+\frac13\right)\pi} = 2^p$$ while avoiding to calculate the loop $F=\prod_{i=1}^{n}$ while $F$ is smaller than $2^p$ (and handling overflows). $p$ itself is a positive integer being a power of $2$.

This is for a computer algorithm that works with on integer lengths, that can be $p = 16$, $32$, $64$, $128$, $256$… ($p$ is a constant, known before the algorithm is evaluated, and is itself $p=2^b$, $b$ number of bits). The actual (unsigned) integer maximal value is actually $2^p-1$, but, to simplify, let's use $2^p$.

Starting from
$$n^ne^{-n}\sqrt{\left(2n+\frac13\right)\pi} = 2^p$$
I thought going with the (natural) $\log$ could be helpful
$$n\log(n)-n+\dfrac 12\log{\left((2n+\frac13)\pi\right)} = p\log (2)$$
but from here…

Maybe there is a better idea to approximate $n$ that costs less than $O(p)$? (approximation of the actual time complexity)

(Even if all $n$ values up to the max $b$ could be pre-calculated, let's assume the general case since we're on Math SE)

Best Answer

I think that it is simpler to just consider the equation in $n$ $$n!=2^p-1$$

For the inverse of the factorial function, @Gary proposed in year $2013$ a superb approximation (have a look here). Applied to your case, it will give $$n \sim \frac{{\log \left( {\frac{2^p-1}{{\sqrt {2\pi } }}} \right)}}{{W\left( {\frac{1}{e}\log \left( {\frac{2^p-1}{{\sqrt {2\pi } }}} \right)} \right)}} - \frac{1}{2}\tag 1$$ where $W(.)$ is Lambert function. Up to you to decide if you want to use $\lceil n \rceil$ or $\lfloor n \rfloor$.

I give below the estimated value of $n$ as well as the exact solution for a few values of $p$

$$\left( \begin{array}{ccc} p & \text{estimate} & \text{solution} \\ 2 & 2.39249 & 2.40587 \\ 4 & 3.67329 & 3.68024 \\ 8 & 5.42819 & 5.43213 \\ 16 & 8.22333 & 8.22553 \\ 32 & 12.8558 & 12.8570 \\ 64 & 20.6665 & 20.6671 \\ 128 & 34.0398 & 34.0401 \\ 256 & 57.2588 & 57.2590 \\ 512 & 98.0766 & 98.0767 \\ 1024 & 170.624 & 170.624 \end{array} \right)$$

If you do not have access to Lambert function, just approximate it using $$W(t)\approx L_1-L_2+\frac{L_2}{L_1}+\frac{L_2(L_2-2)}{2L_1^2}+\frac{L_2(2L_2^2-9L_2+6)}{6L_1^3}+\cdots$$ where $L_1=\log(t)$ and $L_2=\log(L_1)$

Edit

If you want a real shortcut approximation, use for $p=2^k$ $$n = a + b \,k^c$$ A quick and dirty regression for $1 \leq k \leq 20$ (with $R^2=0.999994$) $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & 0.658561 & 0.004641 & \{0.649405,0.667716\} \\ b & 0.261512 & 0.001185 & \{0.259174,0.263851\} \\ c & 1.234300 & 0.001444 & \{1.231450,1.237150\} \\ \end{array}$$