Saturday, October 20, 2012

Factorial of Large Numbers


Hi everyone! I am blogging from Coffee Klatch in San Dimas, CA! I was overdue for a latte. Their pumpkin brownies are to die for.

Factorial of Large Numbers

Calculators have an upper limit when working with numbers. Typically, the scientific calculators allow numbers up to 9.999999999 × 10^99. Anything over 10^100 is considered an overflow. Some of the higher graphing and CAS calculators allow for higher limits.

This means if you need to calculate anything over 69!, most calculators will give an "Overflow" error message.


Factorial:

The factorial of an integer n is the product:

n! = n × (n - 1) × (n - 2) × ... × 3 × 2 × 1



We can remedy this by using the laws of logarithms to calculate factorials of higher integers. Common logarithms (base 10) are used.

n! = n × (n - 1) × (n - 2) × ... × 2 × 1
log (n!) = log [ n × (n - 1) × (n - 2) × ... × 2 × 1 ]
log (n!) = log (n) + log (n - 1) + log (n - 2) + ... + log (2) + log (1)

log (n!) = ∑( log(k) from k = 2 to n)

Let T be this sum.

T = ∑(log (k) from k = 2 to n)


Note the following laws are used:

log (a × b) = log a + log b
log (1) = 0



If we take the anti-common logarithm of both sides, we get:

n! = 10^(∑( log(k) from k = 2 to n)
n! = 10^(T)
n! = 10^(frac(T) + int(T))
n! = 10^(frac(T)) × 10^(int(T))

Where:

(1) int(T) is the integer part of T. To find it, take the integer and ignore everything after the decimal point.

(2) frac(T) is the fractional part of T. To find this, take everything after the decimal point and ignore the integer part.

(3) int(T) + frac(T) = T

Example:
int(14.125) = 14
frac(14.125) = 0.125

A lot of programming and graphing scientific calculators have the int and frac functions. On most Hewlett Packard programming and scientific calculators, they are named IP and FP, respectively.

Example 1: 8!

Calculate 8! using the logarithm method.

T = ∑( log(k) for k = 2 to 8) = log 2 + log 3 + log 4 + log 5 + log 6 + log 7 + log 8 ≈ 4.605520523

In this case, int(T) = 4 and frac(T) = 0.605520523.

Then 8! ≈ 10^(0.65520523) × 10^4 ≈ 4.031999996 × 10^4

Because the accuracy of calculators and computer programs, I recommend that you round 10^(frac(T)) to a set amount of decimal places. I find that 5 decimal places works well. Then:

8! ≈ 4.03200 × 10^4

In fact, 8! = 40,320

Example 2: 111!

1. Calculate the sum ∑( log(k) from k=2 to 111). We find the answer to be 180.2462406.

2. Split the integer and fractional portions. Take the anti-common logarithm of the fractional portion. Hence, 10^0.2462406 ≈ 1.762952551. I round this answer to five decimal places to get 1.76295.

3. Answer: 111! ≈ 1.76295 × 10^180

In fact, using an HP 50g, 111! =
1,762,952,255,109,024,466,387,216,104,710,707,578,876,140,953,602,656,551,604,157,406,334,734,695,508,724,831,643,655,557,459,846,231,577,319,604,766,283,797,891,314,584,749,719,987,162,332,009,625,414,533,120,000,000,000,000,000,000,000,000

And yes, typing out this long number is not as easy it sounds.

Programming

Now we come to the programming portion of calculating n! using the logarithmic method.

This program was done on a HP 39gii, but can easily adapted on most graphing/programming calculators.

Program LFACT

EXPORT LFACT(N)
BEGIN
LOCAL T,M;
Σ(LOG(K)),K,2,N)→T;
ROUND(10^FRAC(T),5)→M;
MSGBOX(M + "10^" + INT(T));
RETURN {M, INT(T)};
END;


Here is a program in BASIC language. I used the App techBASIC from ByteWorks.
Everything after explanation points signify comments.

! Approximate Factorial (5 decimal mantissa plus exponent)
! Eddie Shore
! 10/20/2012

! To improve accuracy, set all variables as DOUBLE
01 DIM N as DOUBLE
02 DIM T as DOUBLE
03 DIM F as DOUBLE
04 DIM M as DOUBLE
05 DIM K as DOUBLE


10 PRINT "Factorial Program"
11 INPUT "N = "; N
12 T = 0

! Sum Loop
14 FOR K = 2 TO N
16 T = T + LOG10(K)
18 NEXT K

! Separation of T into its fractional and integer parts, we know T>0
20 M = INT(T)
22 F = T - M

! Round Fraction to five places - we have no ROUND function to work with
! 10^x when x<1 has a value less than 10
! Use string functions to round the answer to 5 decimal places (7 characters: the integer, decimal point, and five decimal places)
30 F = 10^F
32 F = VAL(LEFT(STR(F),7))

! Final results
40 PRINT "N! = "; F ; " x 10^"; M


Test Data:
8! = 4.0320 × 10^4
66! = 5.4434 × 10^92
122! = 9.87504 × 10^202
425! = 1.59775 × 10^934

Thank you as always, I appreciate the comments and compliments. Hope you enjoy the weekend.

Until next time,

Eddie


This blog is property of Edward Shore, 2012.



Casio fx-9750GIII and fx-CG 50: Playing Games with the Probability Simulation Mode

Casio fx-9750GIII and fx-CG 50: Playing Games with the Probability Simulation Mode The Probability Simulation add-in has six type...