brainwagon

In fact, at the moment, it’s the largest known prime, with over 9.8 million digits. As part of my pi day celebration yesterday, I was trying to review how I might speed up my C code which calculates pi to large numbers of digits. Most of the fast ways rely on fast multiplication, utilizing the FFT algorithm. I wasn’t sure how that really worked.

So… I decided to write a program to learn and test my knowledge.

Rather than compute pi though, I thought I might try a somewhat different but similar task that relied on big number computation. I remembered that Fabrice Bellard had an obfuscated C code entry that used an integer Fast Fourier Transform to print out the largest (then) known prime. It’s really quite nifty. I decided to try to implement a similar thing, but rather than starting with his obfuscated code, I decided that I’d try to use the FFTW library to do the same.

It’s a tiny bit tricky: after all, we are using a floating point FFT to multiply large integers, and getting the precision issues nailed isn’t entirely obvious. My original idea was to represent the large numbers in base 100000 (giving five decimal digits per place) but for reasons which aren’t entirely clear to me, when computing the current largest prime, I seemed to run out of precision. Limiting myself to base 10000 seems to have cured the problem.

My program seems to work (at least minimally) now. I can compute Bellard’s number (2**6972593-1) in about 9.7 seconds, and the largest known prime (2**32582657-1) in about 53 seconds. Ironically, such numbers are really easy to output in base two: they are just lots of ones. But to convert them to base ten for formatting in the way us humans like is much more complicated.

The final output looks like:

 12,457,502,601,536,945,540,085,550,157,479,950,312,279,598,515,115,184
284,367,047,566,259,111,523,599,739,738,055,975,960,661,684,593,910,041
988,688,211,130,870,620,428,490,430,485,634,271,939,241,796,764,631,759
...
... 181631 lines deleted for brevity
...
826,726,604,958,937,322,582,512,072,612,621,443,114,535,641,869,584,273
577,446,330,457,465,821,333,212,445,737,104,635,692,000,092,659,011,752
880,154,053,967,871

I have double checked this against the official value on mersenne.org, and it matches precisely. I know that the program hides at least one gross inefficiency still: I suspect that I can actually compute the large prime in just about twelve seconds. But my head hurts for the moment, and I think I’ll pass for now.