17

I need to find the smallest value of $n$ such that $$\frac{a^n}{n!}\leq 10^{-k}$$ in which $a$ and $k$ are given (these can be large numbers).

I set the problem as : solve for $n$ the equation $$n!=a^n\, 10^k$$ I used for $n!$ Stirling approximation in which I ignored the $\sqrt n$ term. This gives an upper bound of the solution I rewrote as $$n_0=a\,e\, \frac A {W(A)}$$ $W$ being Lambert function and $$A=\frac{k \log (10)- \log (\sqrt{2 \pi} )}{a\,e }$$ which is not too bad (for example, using $k=1000$, $a=100$, the approximation gives $n_0\approx 1402.65$, the algebraic solution being $n\approx 1401.27$).

At this point, I could go backwards and find the solution. This is what I should call brute force.

For sure, starting from this estimate $n_0$, I could perform one iteration of Newton method and get $$n_1=n_0-\frac{\log (\Gamma (n_0+1))-(n_0 \log (a)+k \log (10))}{\psi (n_0+1)-\log (a)}$$ which is still an upper bound (Darboux theorem). Similarly, I could do the same using Stirling approximation and get $$n_1=n_0-\frac{\left(n_0+\frac{1}{2}\right) \log (n_0)-n_0 (1+\log (a))-k \log (10)+ \log (\sqrt{ 2 \pi })}{\frac{1}{2 n_0}+\log (n_0)-\log (a)}$$

However, using the rigorous formulation, this makes the second calculation quite expensive for little improvement and I wonder if something simpler could be considered.

Edit

It is sure that for a given value of $k$, the problem would be much simpler since I could perform a least square fit for a model $$n=\frac{\alpha}{W\big(\frac{\beta}a\big)}$$ and the results are quite good; for example, if $k=1000$, $\alpha=2298.64$, $\beta=845.965$ (to be compared to $\alpha_0=2301.67$ and $\beta_0=846.736$ from the initial model). For $a=100$, this would give $n=1401.28$ which is the answer. The problem is that many $k$'s have to be considered and curve fit does not look to be a solution.

Claude Leibovici
  • 260,315
  • If a is not divisible by 10 then you can get the number of 5's in n! by the value of k. – user94300 Jun 21 '15 at 08:49
  • Could you please elaborate ? Thanks. – Claude Leibovici Jun 21 '15 at 08:50
  • If k is 10, then we know that the number 5 comes 10 times in n!. Hence n is between 45 and 49 (including both). It can be a start. – user94300 Jun 21 '15 at 08:54
  • You have $n \approx 1401.27$, but you also have $n!$. Does this mean you want to to use the Gamma function rather than the factorial? – wythagoras Jun 21 '15 at 09:19
  • @wythagoras. This is just to make the problem continuous. For sure the solution is $1402$. I just use Gamma for the Newton step. – Claude Leibovici Jun 21 '15 at 09:21
  • I realize that this is an old question which may no longer have any relevance for you. I am curious about the background and context. I recognize the connection to the problem of evaluating the exponential function. Normally I would only use a Taylor series on a compact interval and then extend to the positive real line via repeated squaring. In this context your problem comes up and is rapidly solved by trial and error. So, if you have the time I would like to learn about your application. Kind regards. – Carl Christian Feb 07 '16 at 17:48
  • @CarlChristian. The problem is still alive ! In fact, it is a simplified version of a problem in molecular dynamics, statistical thermodynamics and molecular simulation. It is not so difficult except that the equation has to be solved zillions of time in a single simulation and any saving is important by the end. One of the issues is that $k$ varies very fast and can be very large (it is a stiff function of system temperature). Any idea and suggestion are welcome. Thanks for your interest. Cheers. – Claude Leibovici Feb 08 '16 at 06:19
  • @ClaudeLeibovici: Perhaps this idea can be used. Construct a table of $1/n!$ from $n=1$ to $n=170$ as $1/171!$ underflows in double precision. Given $a$, $k$ do a binary search in this table to solve the inequality. I would not program this search using loops but instead write all 8 comparisons out one after another ($2^7 < 170 \leq 2^8$). This can give longer, but also faster code as you remove the instructions necessary for loops and there is no branch misprediction. I think, but I am not certain, that vectorization is possible, so that you can do multiple values at the same time. – Carl Christian Feb 08 '16 at 10:16
  • @CarlChristian. The equation must be solved as $\log(n!)=\log(a^n 10^k)$ since $k$ and $n$ can be several thousands. Again, there is not much problem using Newton provided a good estimate. What I built ia an approximation of $\frac{W(A)} A$ which is a smooth function. The problem would be to find a better starting point to save a couple of iterations. – Claude Leibovici Feb 08 '16 at 10:53
  • @ClaudeLeibovici: Thank you very for explaining the problem to me. I will think about it, but I can not promise any improvement. I still recommend that you unroll the (Newton) loop completely, writing the code as a purely sequential list of arithmetic instructions. This is known as "loop free code", "straightline programs" or "the compiled code approach". I have used this approach to good effect when solving bond constraint equations in the context of MD using Newton's method and a sparse solver. Kind regards. – Carl Christian Feb 08 '16 at 11:12
  • @CarlChristian. It is funny to hear about "loop free code" ! It is one of the first things I learnt in 1960 : just repeat instructions, avoid functions and subroutines as much as you can and so on ! For some industrial codes I wrote, where efficiency and speed was required, I did that a lot (copy/paste is so easy). Even today, there are things I write in assembler ! Is there any place where I could have a look to your work ? Cheers. – Claude Leibovici Feb 08 '16 at 11:55
  • @ClaudeLeibovici: Our first paper on this was accepted by PPAM 2015 in December 2015 You can download a copy from my website. This is the direct link http://www8.cs.umu.se/~spock/my_papers/accelerating_newton.pdf I apologize in advance if I have missed central references. I have no contacts in industry and the people who do sparse solvers in academia have "long" ago moved away from loop free code for tiny systems to BLAS3 based code for gigantic sparse linear systems. Kind regards. – Carl Christian Feb 08 '16 at 18:48

1 Answers1

21

Stirling's Formula says $$ \log\left(\frac{a^n}{n!}\right) \sim -n\log\left(\frac n{ea}\right)-\frac12\log(2\pi n)\tag{1} $$ Multiplying by $-\frac1{ea}$ and subtracting $\frac1{2ea}\log(2\pi ea)$, we get $$ \begin{align} -\frac1{ea}\log\left(\frac{a^n}{n!}\sqrt{2\pi ea}\right) &\sim\frac{n}{ea}\log\left(\frac n{ea}\right)+\frac1{2ea}\log\left(\frac n{ea}\right)\\ &=\frac{n+\frac12}{ea}\log\left(\frac n{ea}\right)\\ &\approx\frac{n+\frac12}{ea}\log\left(\frac{n+\frac12}{ea}\right)-\frac1{2ea}\tag{2} \end{align} $$ Applying the Lambert W Function to $(2)$ gives $$ n\sim ea\exp\left(\operatorname{W}\left(-\frac1{ea}\log\left(\frac{a^n}{n!}\sqrt{2\pi a}\right)\right)\right)-\frac12\tag{3} $$ Plugging in $\log\left(\frac{a^n}{n!}\right)=-k\log(10)$ yields $$ \bbox[5px,border:2px solid #C0A000]{n\sim ea\exp\left(\operatorname{W}\left(\frac k{ea}\log(10)-\frac1{2ea}\log(2\pi a)\right)\right)-\frac12}\tag{4} $$ For $k=1000$ and $a=100$, $(4)$ gives $1401.274188$

robjohn
  • 345,667
  • This is an elegant and beautiful solution for sure ! – Claude Leibovici Dec 31 '16 at 08:38
  • 3
    I just finished a huge testing of your solution : I shall not ask you to guess what. It works like a charm ! Thanks again. – Claude Leibovici Dec 31 '16 at 14:04
  • 1
    You're welcome! I'm glad it works well. Other than Stirling, about the only approximation was $\log\left(1+\frac1{2n}\right)\approx\frac1{2n}$. There's not much to add inaccuracy. – robjohn Dec 31 '16 at 14:17
  • 1
    Happy New Year !. I suppose that you noticed that, setting $a=1$ your formule gives the inverse gamma function and the formula is very close to what Cantrell proposed. – Claude Leibovici Jan 01 '17 at 05:54
  • 2
    $\unicode{x1F389}\unicode{x1F38A}\unicode{x1F55B}\unicode{x1F4A5}\unicode{x1F37E}$ Happy New Year!! I just found David Cantrell's post – robjohn Jan 01 '17 at 07:19
  • 1
    More I see and use this answer, more I think it could be worth to publish it for the general solution of $a^n,n!=k$ which could write $$n=\frac{1}{2} \left(\frac{\log \left(\frac{a k^2}{2 \pi }\right)}{W\left(\frac{a \log \left(\frac{a k^2}{2 \pi }\right)}{2 e}\right)}-1\right)$$ – Claude Leibovici Nov 09 '17 at 08:06
  • May be, for the inverse factorial $(a=1)$, you could be interested by https://ir.lib.uwo.ca/cgi/viewcontent.cgi?article=7340&context=etd which is very recent. There is also this paper https://www.tandfonline.com/doi/full/10.1080/00029890.2018.1420983?scroll=top&needAccess=true which I cannot access. The journal seems to be dedicated to the topic. – Claude Leibovici Dec 01 '18 at 04:07
  • @ClaudeLeibovici: I believe $(2.12)$ on page $23$ of Properties and Computation of the Inverse of the Gamma function is identical to $(1)$ from this answer using $\exp(W(x))=\frac{x}{W(x)}$, which is the $a=1$ case of this answer. Thanks for the references! – robjohn Dec 01 '18 at 06:02
  • @ClaudeLeibovici: that is also $(67)$ on page $418$ of Gamma and Factorial in the Monthly. I was able to access this through UCLA. – robjohn Dec 01 '18 at 06:15
  • Yes, I noticed when I saw $(2.12)$. This is why I did send you my comment. What do you think about $(2.13)$ ? Unfortunately, I cannot access the journal. I have used your magnificent (this is the qualifier I use) approximation in several answers here on MSE. Cheers. – Claude Leibovici Dec 01 '18 at 06:39
  • $(2.13)$ seems a lot more accurate, especially for larger arguments. However, the simplicity of $(2.12)$ is attractive. – robjohn Dec 01 '18 at 06:49
  • Sorry to disturb you again. I have problems using equation $(2.13)$ of the thesis. It is written as $$y \sim \frac{1}{2}+\frac{1}{24 u (w+1)}-\frac{14 (w+1)^2+10 (w+1)+5}{5760 u^3 (w+1)^3}$$ with $u=\frac{\log \left(x\sqrt{2 \pi } \right)}{w}$ and $w=W\left(\frac{\log \left(\frac{x}{\sqrt{2 \pi }}\right)}{e}\right)$. If I make $x=10$, I get $0.503394$. Do you mind to compare with the paper ? Thanks a lot. – Claude Leibovici Dec 01 '18 at 16:12
  • 1
    Sorry for having written my comments in the wrong post. I did suspect this unpleasant situation. Having been myself twice the victim of plagiarism, I wrote to the editor to claim the anterioity of my work. In the first case, the editor added a note in the journal; in the second case, the paper was removed from the journal and all references on the Internet were removed. Be sure that I shall render unto Caesar the things which are Caesar's and only refer to your work for this topic. Cheers :-) – Claude Leibovici Dec 12 '19 at 03:42
  • @ClaudeLeibovici: I just did a Google search on "inverse gamma function" and found this answer by Gary posted on 6 Aug 2013, which predates my answer by over 3 years. So I cannot claim to be the first to have come up with this formula. – robjohn Dec 12 '19 at 08:15
  • Is Gary mentionded in the paper ? – Claude Leibovici Dec 12 '19 at 09:20
  • @ClaudeLeibovici: No. I don't see any references to StackExchange or MathOverflow in either paper. – robjohn Dec 12 '19 at 09:41
  • Past Monday, I gave a lecture in my former research group and I mentioned your solution for the equation. They were so impressed that they decided to use it systematically since they have a super-fast subroutine for the calculation of $W(x)$; as a result, no more need of any IF check in the loop and this saves a lot of CPU time for their applications. As an example, I did show them https://math.stackexchange.com/questions/3450135/approximate-frac-sin-pi-xx-with-an-error-of-less-than-10-6/3450206#3450206 In their names, and mine for sure, thank you very much. – Claude Leibovici Dec 20 '19 at 14:45