# (Numerical Analysis) What is the largest double float represented for the gamma function and $n!$

Consider that \begin{align} \Gamma(n+1) = n! \end{align} for any integers. I then got the following two questions:

1. What is the largest value of $$n$$ for which $$Γ(n+1)$$ and $$n!$$ can be exactly represented by a double-precision floating-point number?
2. What is the largest value of $$n$$ for which $$Γ(n+1)$$ and $$n!$$ can be approximately represented by a double-precision floating-point number that does not overflow?

So I have been working on this question for some time now but I do not seem to be satisfied by my own answer. My thoughts so far on (1.) is that since we start getting non representable integers after $$2^{53}$$ my initial guess would be the largest $$n$$ such that $$n!=\Gamma(n+1)<2^{53}$$. I found this value to be $$n=18$$ through MatLab simulations, I do wonder though if there are still some precisely represented for $$n>18$$, and how would I go about testing that?

I have run $$\Gamma(n+1) - n!$$ for several $$n>18$$ and they start getting some serious errors but only after $$n=22$$. It makes intuitive sense to me that the gamma function should produce more errors as it is a an approximation of an integral, but how am I to know that $$23!$$ is not also represented?

I also considered finding zeros for $$0=n! - (1+f)\cdot2^e$$ w.r.t $$f$$, note that $$(1+f)$$ is the mantissa and $$2^e$$ is the radix, such that the second term is in fact the positive normalized representation of the double float. My thought is that it would only be able to find a zero in case $$f$$ was a representable as a double float, however, I figured that this could not be done like this as $$n!$$ would always be either true or a rounding to a double float, so in any case MatLab would always be able to find a representable $$f$$ that produced a zero.

I am not sure how to continue with this exercise, any good pointers?

• The Gamma function is not evaluated by numerical integration. In fact, the accuracy might be better for larger values.
– user16034
Commented Feb 22, 2022 at 22:18

## 2 Answers

Nice one. Using the standard ieee 754 double precision floating point format, we can exactly represent integers less than $$2^{53}$$, multiplied by a power of 2. When we try to find the largest n where n! can be represented that way, then for example the factor 12 of 12! equals $$3 \cdot 2^2$$. If my calculator gets it right then the largest n is 22, with $$22! ≈ 1.124 \cdot 10^{21}$$; 22! has the factor $$2^{19}$$, and $$22! / 2^{19} < 2^{51}$$.

Consider the table below, showing naturals and the binary representation of their factorial with trailing bits stripped off.

IEEE single/double/quadruple precision are able to represent $$13!$$/$$22!/37!$$ exactly.

$$\begin{array} &0! & 1 & 1 \\ 1! & 1 & 1 \\ 2! & 1 & 1 \\ 3! & 2 & 11 \\ 4! & 2 & 11 \\ 5! & 4 & 1111 \\ 6! & 6 & 101101 \\ 7! & 9 & 100111011 \\ 8! & 9 & 100111011 \\ 9! & 12 & 101100010011 \\ 10! & 14 & 11011101011111 \\ 11! & 18 & 100110000100010101 \\ 12! & 19 & 1110010001100111111 \\ 13! & 23 & 10111001100101000110011 \\ 14! & 26 & 10100010011000011101100101 \\ 15! & 30 & 100110000011101110111011101011 \\ 16! & 30 & 100110000011101110111011101011 \\ 17! & 34 & 1010000110111111011101110110011011 \\ 18! & 37 & 1011010111110111011001100101001110011 \\ 19! & 41 & 11011000000101011100100110000011010001001 \\ 20! & 44 & 10000111000011011001110111110010000010101101 \\ 21! & 48 & 101100010100000111011111010011011010111000110001 \\ 22! & 51 & 111100111011101010010011000010101100111110000011011 \\ 23! & 56 & 10101111001011100001100110101111110001010010011001101101 \\ 24! & 58 & 1000001101100010100100110100001111010011110111001101000111\\ \cdots\end{array}$$

Regarding the non-overflowing value, it suffices to accumulate the base-2 logarithms of the naturals until you reach $$128/1024/16384$$. This can be done by brute-force. We have

$$\begin{array} &\log_234! \sim 127.79512061296911\\ \log_2170! \sim 1019.3694529277266\\ \log_21754! \sim 16378.090470312889\\ \end{array}$$ and larger values will overflow.

My old TI SR-56 supported factorials up to $$69!$$ (overflow at $$10^{100}$$).