# Factorial(1000) as number

13 views (last 30 days)
Khannanov Shamil on 14 Apr 2021
Answered: William Rose on 15 Apr 2021
I want to be able to see Factorial(1000) in the form of
x.xxxxe+yyy
but it only returns
inf
How can I display such large numbers?

Walter Roberson on 14 Apr 2021
vpa(factorial(sym(1000)))
ans =
4.02387260077093773543702433923e+2567

William Rose on 15 Apr 2021
You want to calculate
n! ~ sqrt(2*pi*n) * (n/e)^n (Stirling's approximation)
but it is too big when n>170, so what do you do?
I also worked out the following answer, which is to figure out the integer base-10 exponent, and the corresponding mantissa, using the standard rules for powers and logarithms:
M1=sqrt(2*pi*n); E1=n*log10(n/exp(1));
M2 = M1*10^(E1-floor(E1)); E2=floor(n*log10(n/exp(1)));
while M2>10, M2=M2/10; E2=E2+1; end
fprintf('n! ~= %.3fx10^%d\n',M2,E2);
Line 1 computes the mantissa and base-10 exponent, but the exponent is generally not an integer.
Line 2 makes the exponent an integer and adjust the mantissa accoridngly.
Line 3 divides the mantissa by 10, and raises the exponent by 1, repeatedly, until the mantissa is <10.
Line 4 prints the result.
When you compare the result above to the actual factorial, for 10<=n<=170, you see that the Stirling approximation is good, and gets better as n gets bigger. For n=1000, the above code gives
n! ~= 4.024x10^2567
which matches the value given by @Walter Roberson.

the cyclist on 14 Apr 2021
The largest finite floating-point number in double precision is
realmax
ans = 1.7977e+308
For reference
factorial(170)
ans = 7.2574e+306
and
factorial(171)
ans = Inf
Often, applications using these very large numbers end up dividing them by each other, so perhaps you can use logarithms for intermediate calculations? Perhaps you could take advantage of Stirling's approximation somehow?

William Rose on 14 Apr 2021
@Khannanov Shamil, I agree with the suggesiton to use Sritlin's approximation and just report the log of the factorial:
where the term O(ln n) can be ignored.
Or use this form of Stirling:
and ignore the term O(1/n). Don't try to actually calculate n! using the formula above, because it will be too big, but you can work out the mantissa and the base-10 exponent, using the usual rules of algebra, and print the result using fprintf() intelligently.