Now on ScienceBlogs: The Laboratory at Harvard

Seed Media Group

Built on Facts

An exploration of physics and the quest to understand our world.

Profile

profile.jpg Matt Springer is a graduate student of physics at Texas A&M university. He is also an occasional writer and tinkerer, and he is probably too curious for his own good.

Search

Feeds/Networks

http://www.wikio.com

Donate

Help Matt not starve! Use this link to amazon.com when you order from Amazon, and a fraction of the purchase price will be sent to me at zero cost to you. Much obliged, and thanks for your patronage.


Recent Posts

Recent Comments

Archives

Physics/Math Blogs

My Other-Than-Science Reading (variable, very incomplete)

« Ghosts and Kamikazes | Main | Backyard Nuclear Reactors »

Sunday Function

Category: Sunday Function
Posted on: November 9, 2008 10:00 AM, by Matt Springer

Math-averse readers! Do not be scared off! You can enjoy this entry even if as far as you're concerned the equations are pretty pictures of Cypriot syllabary.

Not long ago we looked at adding up lots of consecutive integers. Multiplying consecutive integers is also interesting, and not only that it has a tremendous number of uses all throughout physics and pure mathematics. The function that multiplies the integers from 1 to n is called the factorial function, and rather unusually it's denoted with an exclamation point, like this:

1.png

Here I'm using dots to denote multiplication, as is pretty common in physics. Now the factorial function gets tremendously huge tremendously fast. The factorial of one hundred is almost 10158, which is just mind-bogglingly enormous. Now if you want to find the factorial of one million, you're absolutely hosed. Your calculator will choke, and I'm pretty sure even sophisticated systems like Mathematica will start to smoke and sputter before finally giving you an "out of memory" error. There's just no good way to multiply a million integers together, let alone finding the factorial of the much larger numbers in physics that crop up all the time. Maybe we can find a clean approximation which is good for large n? We can.

Since n! is so unwieldy by virtue of its size, let's chop it down a notch by finding its natural logarithm. We'll use this old fact from high school algebra:

2.png

Knowing that, we can write the log of n!

3.png

That's not an approximation, so far it's exact. But exact though it is, it's not very useful. There's a trick you can use in calculus to approximate a sum with an integral, so we'll go ahead and use it.

4.png

Now evaluate the integral:

5.png

Now the number 1 is pretty much by definition negligible in comparison to large n, so we can drop it from the above equation. From there, take the exponential of both sides to recover an expression for the factorial instead of the logarithm of the factorial.

6.png

It's just an approximation, but it's a good one. Now it turns out that some careful finagling with the limits of integration in the integral (or alternately an entirely different argument involving the gamma function) can improve our estimate. I'll skip the derivation and present the final result. It's called Stirling's approximation, and it's our Sunday Function:

7.png

You might think that since the root n term grows without limit, our original approximation can't be right. But that's not the case. The nn term is so much more quickly growing that for many practical purposes it still works fine. After all, really this is an approximation of the logarithm of n factorial anyway. Neither our rough approximation nor the better one actually converges on n!, instead the absolute error merely grows more slowly than n!. That said, the relative error goes to zero using the root n term. It doesn't if you leave out the root n term.

Those slightly technical matters aside, both expressions are very useful. How about finding that factorial of one million using it? Let's do it:

8.png

Evaluating the first term is easy; the second will make your calculator choke. Some simple tricks involving logarithms will make it very easy though. I'll save those as an exercise for the alert reader and present the result:

9.png

Whew! Now that is a large number, containing more than five and a half million digits. Readers who have done the calculation themselves will have noticed that as we mentioned before the square root term affects the result in this stacked exponential notation not in the slightest. To get that answer took me about 60 seconds worth of writing down a couple lines to keep track of the logarithms and typing the figures into the calculator. To find the answer without Stirling's approximation would have been utterly beyond my grasp. Thanks, Stirling!

Share this: Stumbleupon Reddit Email + More

TrackBacks

TrackBack URL for this entry: http://scienceblogs.com/mt/pings/85382

Comments

1

For what it's worth, Mathematica actually calculates the approximation to be 8.2639*10^(5565708) in 1.735 seconds, and a separate run got the full value in 1.703 seconds.

Though, naturally it only showed me part of the full output. When I asked it for the full value, it broke into tears.

I heart Mathematica!

Posted by: Foxy | November 9, 2008 12:36 PM

2

Actually, check that, it gave the full output after a few seconds, no sweat : D

Posted by: Foxy | November 9, 2008 12:38 PM

3

That will teach me to doubt Mathematica! But we can crank up the pain and try again. How about the factorial of a googol (10^100)?

Using Stirling's approximation I get that it's about 10^(10^101.99), so clearly since there's something like a hundred googol zeros required to just to write the number down, full output on a computer screen would be challenging.

Posted by: Matt Springer | November 9, 2008 2:15 PM

4

The free software SAGE, on my Macbook Pro, gave 1000000! in about 15 seconds. Not blindingly fast, but it didn't choke and (did I say this) SAGE is free.
Works great with Latex too.

Posted by: dean | November 9, 2008 2:16 PM

5

'Overflow[]'

You win this one, Mr. Springer! But be careful, Mathematica and I - we'll be waiting : )

Posted by: Foxy | November 9, 2008 3:02 PM

6

Stirling's formula actually generalizes to a nice series. One interesting property is that the series as a whole never converges if you take all the terms but any truncated part of the series forms a really good approximation of n! for large n.

Stirling's formula has some nice applications. For example, it is helpful for getting tight bounds on Chebyshev's theorem on the distribution of the primes without much pain.

Posted by: Joshua Zelinsky | November 9, 2008 5:20 PM

7
That will teach me to doubt Mathematica!

Actually even the brute-force algorithm only takes half an hour on a fast CPU. Say you represent a big integer by an array of 32-bit unsigned integers (2^32 = 4 billion). Then the largest "big integer" you would need would have a length of (log 10^6!) / (log 2^32) ~ 578,000 machine integers (takes 2.2 MB). In the naive factorial algorithm that loops through the integers in the range 1..10^6 and multiplies them to an accumulator,

acc = 1
for i in 1..10^6 {acc = acc * i}

i never exceeds 10^6, so it fits in a single machine integer. So the operation "acc * i" is a multiplication of a big integer and a machine integer - we can do this by looping through the components of the big integer, multiplying by i, and adding the overflows to the left. This is about 2*n arithmetic operations (n the length of the big integer). So in all we need no more than 2 * 578,000 * 1,000,000 floating ops, about 1 trillion. A modern CPU at 20 billion ops/sec (dedicated ALU) would do this in 10^12 / (2*10^10 /s) = one minute, if it weren't for pointer arithmetic and memory access. In my attempt it took half an hour.

The answer is 1.C3EF * (10h)^4,622,221. (10h = 16d - so the exponent (10h)^4,622,221 is similar to 10^5,565,708).

Posted by: easily amused | November 9, 2008 10:14 PM

8

Well, that's amazing. I just did Stirling's approximation on a calculator and got

8.26393134 x 10^5,565,708.


I wrote a short computer program to sum: log10(2) + log10(3) + ... + log10(1000000)

and got 5,565,708.91718666 which means

8.26393052 x 10^5,565,708 (very close!)

It took only a fraction of a second to run. I don't know what Mathematica was doing for Foxy to take 1.7 seconds.

Posted by: Tim Gaede | November 10, 2008 11:31 AM

9

Actually there is a difference in the computation speed of Factorial in Mathematica depending on when you attempt to do the numerical (i.e., floating point) evaluation. So, on my machine for

N[Factorial[(10^6)]]

it takes 0.76 seconds.

But

Factorial[N[10^6]]

takes 0.000093 seconds, giving the same result in both cases. This all has to do with how Mathematica attempts to address exact calculation versus numerical ones. And hence the order in which you set things up can make a difference.

Posted by: David Reiss | November 10, 2008 5:23 PM

Post a Comment

(Email address is entirely optional, but a consistent email - fake is fine - helps the system identify repeat commenters as not spam.)





ScienceBlogs

Search ScienceBlogs:

Go to:

Advertisement
Follow ScienceBlogs on Twitter
Visit the Collective Imagination blog
Advertisement
Enter to win

© 2006-2009 Seed Media Group LLC. ScienceBlogs is a registered trademark of Seed Media Group. All rights reserved.

Sites by Seed Media Group: Seed Media Group | ScienceBlogs | SEEDMAGAZINE.COM