Last year, I wrote a post about a program I wrote to compute a really large prime number. . At the time, 2^{32482657}-1 was “king” of the primes, having over 9.8 million digits. Today, one year later, it turns out that the Great Internet Mersenne Prime Search has discovered two even larger primes. I dusted off my code, and was pleased to find that it still worked, generating the nearly 13 million digits of the new king of primes in roughly two minutes…

Here are the digits of 2^{43112609}-1:

316,470,269,330,255,923,143,453,723, 949,337,516,054,106,188,475,264,644,140, 304,176,732,811,247,493,069,368,692,043, 185,121,611,837,856,726,816,539,985,465, 097,356,123,432,645,179,673,853,590,577,432597 lines deleted...915,084,069,991,159,279,791,908,398,130, 223,304,824,083,119,093,195,998,014,562, 456,347,941,202,195,900,928,079,670,729, 447,921,616,491,887,478,265,780,022,181, 166,697,152,511

Addendum: Okay, it dawns on me that I never explained why this is difficult. Here’s the thing: Mersenne primes are really easy to write in binary: they are just a long string on ones. But to write them out in base 10 like we all know and love is more difficult. To do it, I actually do all my arithmetic in base 1000 (a convenient choice, because we want to print out the number with commas in place). Computing the prime isn’t very hard really, a^{b} is 1 if b equals 0, is a * a^{b-1} if b is odd, and is a^{b/2} squared if b is even. Then, you only need to do log2(b) multiplies to get your (in this case very big) number. The real trick comes in doing the multiply: the obvious algorithm takes O(n^{2}) time to multiply two numbers. But a faster way is to use the FFT. The trick is realizing that multiplying numbers using the old scheme is basically a lot like convolution. Convolution can be changed to simple multiplication by using the Fast Fourier Transform. So, here’s the basic idea. Suppose you have a number represented in base 1000. Pad it with an equal number of zeros. Then take the Fast Fourier Transform. Square each resulting bin. Then inverse FFT it. You’ll be left with some bins with values > 1000. But you can normalize those by dealing with the carries. And you are done! The FFT only takes O(nlogn) time, the normalization is linear, so you are much, much faster in the limit doing multiplies this way.

There are a few issues with precision but overall, it works quite well.