11
$\begingroup$

What is the number of digits (IntegerLength) of the factorial of 12 345 678 987 654 321? The number of zeros at the end of this factorial was calculated and it is huge: exactly 3 086 419 746 913 569 zeros. The calculation of the number of digits for this factorial if done by simple multiplication, will take about 199 years in my current machine which has a I7 with 6 cores and 64GB of memory. Yes, I know, this number will not fit in my memory, even though I have a SSD with 1TB extra memory which will also not be able to store this gigantic number. So, any ideas?

PS: the program to calculate the number of zeros at the end of a factorial is:

Fnz[n_] := Module[{f, s, p, z},
  p = 5;
  f = IntegerDigits[n, 5];
  s = Total[f];
  z = (n - s)/(p - 1);
  Print[n, "! has ", z, " zeros at the end."];
  ]

Here the result:

 Fnz[12345678987654321]

12345678987654321! has 3 086 419 746 913 569 zeros at the end.

$\endgroup$
  • 2
    $\begingroup$ Thanks for the question. One of my pet peeves, however, is people describing their problems as "hard", when what they mean is "I can't work out how to do this" - which is implied, because you posted the question here! :) $\endgroup$ – Patrick Stevens Dec 30 '15 at 18:12
  • $\begingroup$ Patrick, factoring a large number n (with say 300 digits) is also a hard math problem even though everybody knows how to do it. Just try every prime until you reach the Sqrt[n], right? Same thing here. You see, I could calculate the number of zeros at the end of the number WITHOUT having it in my hand to count the zeros. That´s is what I´m talking about. $\endgroup$ – Giorgio Dec 30 '15 at 18:31
  • 1
    $\begingroup$ Factoring is in general difficult, but factoring a factorial is not :) $\endgroup$ – Patrick Stevens Dec 30 '15 at 18:53
  • 1
    $\begingroup$ @Giorgio One possibility: buy 3.10^4 - 1 additional 1To SSD drives. More seriously, did you check Chip Hurst's answer (which gives another result)? $\endgroup$ – anderstood Dec 30 '15 at 19:19
  • 1
    $\begingroup$ What do you mean by "program"? You can always pack the code in my answer into it's own function. Something like Ndf[x_] := Ceiling[LogGamma[N[x + 1]]/Log[10]], and now Ndf[12345678987654321] outputs 193299018111544064. $\endgroup$ – Chip Hurst Dec 30 '15 at 19:47
21
$\begingroup$

We can exploit the built in LogGamma:

x = 12345678987654321;
Ceiling[LogGamma[N[x + 1]]/Log[10]]
193299018111544064

Edit, Addressing precision:

We have naively for $n > 1$, that $n! < n^n$. Taking logs of both sides gives the (not very tight) bound $\log\Gamma(x + 1) < x \log(x)$ for $x > 1$.

This means if we want the number of digits of $n!$ for some large $n$, asking for $\log_{10}(n \log(n))$ number of digits in $\log\Gamma$ will be more than enough.

Here it is in a function.

Ndf[x_] := 
  With[{p = 20 + Max[1, Ceiling[Log10[x Log[x]]]]},
    Ceiling[LogGamma[N[x + 1, p]]/Log[10]]
  ]

Ndf[12345678987654321]
193299018111544064
Ndf[10^30]
29565705518096748172348871081099
$\endgroup$
  • 1
    $\begingroup$ Correct, but at the edge of what can be done with machine precision. $\endgroup$ – Daniel Lichtblau Dec 30 '15 at 19:23
  • $\begingroup$ Yep, I suppose that is worth mentioning. $\endgroup$ – Chip Hurst Dec 30 '15 at 19:24
  • 3
    $\begingroup$ 12345678987654321+1 can just barely be represented exactly as an IEEE double-precision float, but only because it's even! If it were odd, that last 1 bit would be lost or rounded up. There is not even the potential of a loss of precision if you just do Ceiling[LogGamma[12345678987654321 + 1]/Log[10]]. $\endgroup$ – Mark Adler Dec 30 '15 at 22:28
  • $\begingroup$ @MarkAdler That's a good point about using exact inputs, the only issue is that there is a gray area where Mathematica will try to auto evaluate LogGamma, but it would require more memory than most computers have. For example In[2]:= MemoryConstrained[LogGamma[10^10], 10^8] Out[2]= $Aborted. It looks like it stops auto evaluating around 10^14. $\endgroup$ – Chip Hurst Dec 31 '15 at 0:17
  • 1
    $\begingroup$ I don't know what the gray area is, but this doesn't fall into it. With the Ceiling, Mathematica uses less than 4K of memory to do the exact calculation. $\endgroup$ – Mark Adler Dec 31 '15 at 0:22
13
$\begingroup$

I wondered if Chip's answer was exactly correct, given Daniel's comment about machine precision. So I did it a little differently with higher precision in a way that gives good confidence in the result. (It turns out that the machine precision answer is correct.)

LogGamma can be expanded around infinity thusly:

Series[LogGamma[z], {z, Infinity, 4}]

to give:

$$z \left(-\log \left(\frac{1}{z}\right)-1\right)+\frac{1}{2} \left(\log \left(\frac{1}{z}\right)+\log (2 \pi )\right)+\frac{1}{12 z}-\frac{1}{360 z^3}+O\left(\left(\frac{1}{z}\right)^5\right)$$

For large $z$ that should converge very rapidly. In fact, the $1\over z$ terms can be dropped, since $O\left(1\over z\right)$ will be way less than one digit. For really big z, like this one, this is a sufficient approximation for LogGamma to count digits (unless you happen to be right on the hairy edge of an integer):

$$\left(z-{1\over 2}\right)\log z-z+{1\over 2}\log\left(2\pi\right)$$

(Simplified a little.) You can add the $1\over 12z$ if you like, but it won't make any difference here.

Evaluating that for this case, I indeed get:

(z - 1/2) Log[z] - z + 1/2 Log[2 \[Pi]] /. 
  z -> 12345678987654321 + 1 // Ceiling[#/Log[10]] &
193299018111544064

Looking at the digits after the decimal:

(z - 1/2) Log[z] - z + 1/2 Log[2 \[Pi]] /. 
  z -> 12345678987654321 + 1 // N[#/Log[10], 20] &
1.9329901811154406373*10^17

So it's not on the hairy edge since the ...63.73 was rounded up to ...64, and in fact adding the $1\over z$ terms don't change the result.

Interestingly, that formula is not the log of Stirling's formula, and the log of Stirling's formula does not give the correct answer. (It's low by eight digits.)

$\endgroup$
  • $\begingroup$ Since the number of digits in $n!$ is bounded above by $2n\log_{10}(n)$, we could always do Ceiling[LogGamma[N[x + 1, Ceiling[2x Log10[x]]]]/Log[10]] and be assured that our answer is correct. (I may edit my post with something similar later today.) $\endgroup$ – Chip Hurst Dec 30 '15 at 21:47
  • 3
    $\begingroup$ If you use the 2nd version of Stirling's approximation at en.wikipedia.org/wiki/Stirling%27s_approximation which has corrections of order $\ln(n)$, n! ~ Sqrt[2*Pi*n] (n/E)^n, it comes out right, though: Ceiling@N[(Log[Sqrt[2 \[Pi] n]] + n (Log[n] - 1))/Log[10] /. n -> 12345678987654321, 20] yields 193299018111544064. $\endgroup$ – evanb Dec 30 '15 at 22:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.