Factorial(1000) as number

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?

Answer 1

[Walter Roberson]

vpa(factorial(sym(1000)))

ans = 

4.02387260077093773543702433923e+2567

Answer 2

[William Rose]

@Khannanov Shamil,

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?

@Walter Roberson answered the question.

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.

Answer 3

[the cyclist]

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?

Answer 4

[William Rose]

@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.