Fun with Prime Numbers

Steve Litt writes "Fun With Prime Numbers contains a series of prime number finding algorithms starting with the most brute force imaginable, and working up to a paged algorithm capable of finding the first 1,716,050,469 primes in an hour and a half on a commodity machine. There are faster algorithms on the net, but these algorithms are within the reach of mere mortals and are fully explained."

Comments

Obligatory

[zedmelon]

This article will be a prime target for bad jokes.

Re:Obligatory

[Anonymous Coward]

We must not be divided on this issue.

Re:Obligatory

[dtfinch]

0x1337 is prime

Re:Obligatory

[Trejkaz]

Luckily when it comes to the quality of a Slashdot comment, bad jokes don't factor into it.

Re:Obligatory

[eclectro]

This article will be a prime target for bad jokes

It's true. Good jokes come alomg only so often, and we don't know when.

Re:Obligatory

[conan776]

God sieve us from these bad puns....

Re:Obligatory

[Surazal]

Not only don't we know when, we don't remember the last time, either.

Where is it?

[TamMan2000]

I got a 404...

Re:Where is it?

[zedmelon]

Google Cache

Re:Where is it?

[canwaf]

Sorry, 404 isn't a prime number. Try again.

Re:Where is it?

[Baelrun]

I got a 404.. That's not a prime number. No wonder the algorithm is so fast.

This begs the question:

[Anonymous Coward]

Why not just save primes on a disk instead of recalculating them all the time?

Re:This begs the question:

[Anonymous Coward]

Why waste disk space with primes instead of just recalculating them?

Re:This begs the question:

[psetzer]

Actually, in the final implementation, I thought that he does that. It creates a few thousand files, each one containing n primes, and then it goes through each previous file up to the square root of the largest number in the current file, and then it quits and goes to the next file. If you want to duplicate this guy's efforts, I think that the biggest help for this would be a REALLY fast HD array and lots of RAM.

Re:This begs the question:

[Dasein]

There are a lot of them.

As a matter of fact 2048 bits can represent numbers up to 2^2048 ~= 3.23x10^616. By the Prime Number Theorem there is about 2.27x10^613 primes in that space. If each prime required 2048 bits to store, then storing all the 2048 bit primes requires more petabytes then you and I are likely to ever see.

(Jeez I hope I didn't do the math wrong. The most convenient calc I have available to me as I write this in the Windows calc)

Re:This begs the question:

[Omniscient Ferret]

On my computer, generating primes from 2 to 2^32 takes a third of the time to unbzip2 the same list.

Re:This begs the question:

[maxwell demon]

Indeed, the number of particles in the universe is assumed to be approximately 10^85, so if you would store one bit per particle, you'd need about 5*10^530 universes to store all those primes. Hell, even if you stored at every single particle as many bits as there are particles in the universe, you'd still need 5*10^445 universes.

Of course that's assuming you stupidly write all those numbers individually down, using the whole 2048 bitf or each number. Now what if we compress the table?

The complete list can, of course, be stored in a bitfield with 2^2048 individual bits which tell if a given number is prime or not (e.g. 1 means the number is prime, 0 means it's not). Now, of those 3.23*10^616 bits, about 2.27*10^613 bits are 1. So a simple approximation would be to treat the bits as independent, with each one having a probability of (2.27*10^613)/(3.23*10^616)=7.03*10^-4 of being 1. So if we neglect all correlations of the bits, we get an information of about 0.0084 bits per bitfield bit, which means we can compress to about 2.7*10^614 bits. That is, even neglecting all correlations (even the trivial ones of all even numbers except 2 being non-prime) we get only about 12 bits per prime. That's substantially less than the 2048 bits per prime we'd need with straightforward storage. Of course, we would still need an insane number of universes to store our table :-)

Now, of course we know that we can "compress" the full table very well to a small prime-finding program which fits nicely on even a single floppy (but would take longer to "decompress" than we are willing to wait), so the list of primes doesn't actually contain too much information (indeed, the program can compress the whole infinite prime table - decompression of course needs unlimited ressources, then -, therefore the average information content of a single prime is obviously zero).

Re:This begs the question:

[StyroCupMan]

If each number only has a 0.000703 chance of being prime, we can simplify this whole calculation with this function (pseudo-coded):

boolean isPrime(int Number) {
return false;
}

That function is 92.97% accurate. That's an A-. Good enough for me. :)

What is special about prime numbers?

[Dancin_Santa]

I understand that number sequences like Fibonacci manifest themselves in Nature and others like pi provide a fairly decent random number generation method. However, aside from the interesting property that they can't be divided by anything other than themselves and 1, I do not 'get' what is interesting or useful about prime numbers.

But then again, I haven't studied mathematics to any great extent beyond the multi-var calc and linear algebra back in high school. That was such a long time ago. Any math Geniuses out there want to clue is in on why primes are interesting?

Re:What is special about prime numbers?

[MillionthMonkey]

Brief overview. There is also the Prime Spiral.

Re:What is special about prime numbers?

[Jazzer_Techie]

I do not 'get' what is interesting or useful about prime numbers.

Prime numbers are useful in encyption, but those are very large primes that are difficult to factor. To a mathematician, prime numbers are fascinating. The distribution of primes, which is specifially random but generally predictable (check out the MathWorld article for more details) is of particular interest.

For example, the Riemann Zeta Function of n, which is the infinite sum of all terms k^(-n)[k varies from one to infinity], is expressable as the infinite product (1-p1^(-n))*(1-p2^(-n))*(1-p3^(-n))... [where pn is the nth prime number]. This is a mind-boggling connection. For example Zeta(2)= (pi^2)/6 = (1 + 1/2 + 1/4 +1/8 ...)= ((1-2^(-2))*(1-3^(-2))*(1-5^(-2))...) There is therefore a connection between pi and the distribution of the primes. Math is crazy stuff.

chosing typesafety over speed?

[Anonymous Coward]

So the macros knock the time down from 4.49 to 4.15 seconds -- less than a 10% reduction. In my opinion, that's not enough benefit to give up the robustness and typechecking of functions.

I never thought that I'd see a C programmer type something like that.

The author isn't a lightweight ...

[xmas2003]

While the prime number page is a bit odd (obligatory half-hearted pun), the author is Steve Litt who wrote Samba Unleashed ... plus there's already 40+ posts on the article, yet his web site is still pretty snappy despite the /. crowd ...

BTW, the first bazillion prime numbers HAVE been calculated, so for those /.'ers with spare CPU cycles, consider something perhaps more worthwhile such as Folding@HOME

Encryption

[sbszine]

Read Penn Jillette's great explanation for details.

Re:Encryption

[sbszine]

Yeah, but that's an example of them being useful. It doesn't really make them esthetically "interesting".

In that case I refer you to the fundamental theorum of arithmetic, which I think is pretty interesting.

Have some more fun with primes

[pestilence4hr]

A very fast segmented sieve (to find primes in an arbitrary range) can be found here, by Tomás Oliveira e Silva.

Of tangential interest is his Goldbach conjecture verification project. You can download his client and contribute to the project. He's shooting for 10^18...

the factor command in Unix/Linux

[3770]

There is a command in many unices and Linux called factor.

If you want to be a true geek you can try it on your friends phone numbers and find out if it is a prime.

Then, the next time you talk, inform them that their phone number is a prime, or tell them of their phone numbers prime factorization, and enjoy watching them think that you are a super geek and a super genius.

For even better effect, pretend to count in your head before you tell them this.

Re:the factor command in Unix/Linux

[0racle]

You don't really have friends do you.

Re:the factor command in Unix/Linux

[bob beta]

If you want to be a true geek, you run factor against all the spare quartz crystals and oscillator blocks in your junkbox, to see which ones factor down to useful frequencies. 'A multiple of 9600 baud, hmmm....'

Re:This begs the question

[YGingras]

What's the 1,716,050,470th prime?

The 1,716,050,470th prime is 40,100,000,093.

Find out more on the Nth Prime Page.

Gaps between primes

[tbjw]

From the article:

who can be sure that there aren't two consecutive primes somewhere that are more than 65536 apart

Let's recall elementary number theory:

65537! + 2 = 2 * 3 * ... * 65537 + 2 = 2 * (3 * ... * 65537 + 1)
65537! + 3 = 2 * 3 * ... * 65537 + 3 = 3 * (2* 4 * ... * 65537 + 1)
...

65537! + 65537 = 2* 3 * ... * 65536 * 65537 + 65537 = 65537(65536!+1)

So the 65536 consecutive numbers 65537! +2 , ..., 65537! + 65537 are all nonprime (the first is divisible by 2, the second by 3, the thrid by 4 etc). Since we also know that there are infinitely many primes, there must be a prime greater than 65537! + 65537. The biggest prime less than 65537! + 2 and the least prime bigger than 65537! + 65537 are consecutive primes which are at least 65536 apart.

Of course, 65537! is a massive number, which is unlikely to crop up in these calculations. There may be other sufficiently large gaps between primes among smaller numbers.
Consider, for instance, the method above applied with 5 in the place of 65536. We see that {6! +2, ..., 6! + 5} ={722, 723, 724, 725, 726} are all nonprime, but also 24, 25, 26, 27, 28 are 5 consecutive nonprime numbers.
For where we can expect large 'prime-free' gaps, I'll defer to any number theorist.

Prime Numbers for cash!

[kuzb]

If you're interested in running a distributed.net-like application with the chance to make some money using your CPU cycles to help find prime numbers, have a look at here. They give some good information about what prime numbers are, and what they're good for.

Fun with primes?

[The Master Control P]

How about when you're on the vintage mainframe level in Tron 2.0, ask some random program to calculate the seventh even prime number. When he segfaults, you get access to the directory containing Tron Legacy. Just don't ask yourself that question...

Improvement

[acidblood]

Blockquoth the FA's author, regarding how to store a list of primes:

Keeping it in an array is simplest, but one must declare the array before finding primes. How big do you declare it.

Have a look here for a pretty tight upper bound on the number of primes up to a given number. Using an array, instead of a linked list, would probably lead to a small speed improvement on his code.

He could also use an std::vector from C++. As far as I can tell it's pretty easy to resize.

Re:A better use?

[benna]

encryption.

Not as much fun as

[Snoe]

The Prime Number Shitting Bear. Watching a console window spit out prime numbers might do it for the geeks, but everyone can loves the prime number shitting bear.

This is actually the MOST important thing to do!

[Theovon]

For those who scoff at this kind of stuff, we have to keep in mind that it's this sort of tinkering that results in some of the most amazingly useful and interesting stuff later.

I have one lament, which is that his code formatting (indenting) got munged, making it really hard to read the code. I'm really tired, so I just couldn't manage to read through all the C code. I hope he fixes the page so that I can read it more easily later.

Better Implementation

[jkleint]

As mentioned, DJB's primegen is much faster. Even his Sieve of Eratosthenes implementation will do the primes up to a billion in 1.38 seconds on a 2.4 GHz P4, which is roughly 100x faster than this guy's "all primes below 100 million in less than 13 seconds" (although he doesn't specify on what machine). There's no reason to be using huge memory "pages," you can sieve up to N with any size interval you want (like say, the size of your L1 d-cache) using only Sqrt(N) memory for the primes up to Sqrt(N).

I just RTFA...

[acidblood]

...and I'm going to list a few improvements here that I think the author missed.

In one of the early programs, the author had the idea of skipping all even numbers, which are of course divisible by 2. In the next program he found out how he could skip numbers divisible by 3. One can make a systematic method out of that -- in the literature it's known as adding a wheel to the sieve of Erathostenes. Thing is, he implemented the sieve of Erathostenes foregoing this improvement he had implemented earlier. This leads to a big improvement: when using a wheel with the first four primes 2,3,5,7, there are only 48 residue classes out of 210 = 2*3*5*7 possible that are not divisible by any of these four primes. In practice this means that, for each 210 numbers you'd consider, the sieve would only need be applied in 48 of them. That's a gain of 4.375. One can make a gain of nearly 5 using a wheel of the first five primes, but it begins to get complicated.

Another improvement would be to make the program cache-friendly. He sort of did this, down one level of the memory hierarchy, by paging the data to disk -- now he only needed to realize that, just as memory is much faster than disk, cache memory is faster than main memory. I've developed a sieve employing this optimization and it pays off indeed. This is considered a basic technique in the realm of high-performance computing. (The NFSNET, see my sig, uses a sieve including this and other optimizations, for instance).

I might get flamed for this, but he should consider assembly programming. There aren't many applications where assembly will pay off so handsomely as this one, and the `core' of the sieve is very small, so only a few lines of assembly are needed. Vector integer instructions, such as MMX and the integer instructions of SSE2, can be employed to very good advantage here. My sieve program also included this, and it gives a real nice boost -- and only about 25 lines of assembly were used, the rest was C.

Lastly there's the issue of representing the output. He stored the primes on disk by writing their values. A little compression can be realized in this representation: for instance he's printing out the primes in decimal, while hexadecimal would be more space-efficient. And instead of printing out the ASCII characters from 0 to 9, A to F, he could use a binary representation and read them through an hex dump -- that saves half the space, as a byte is used to store two hexadecimal digits instead of one.

A perfectly good representation that could save a lot of space (8 times, to be precise) would be a bitmap, which he used on some of his programs. That's where you represent primality or lack of it by setting and clearing bits of an integer. It's pretty easy to determine a prime from a bitmap value, so this representation is just as good as the previous one. One could also use the wheel idea I explained before to save even more space: a concrete example would be that for each block of 210 numbers, one needs only bitmaps for 48 of them, as the rest are divisible by the first four primes. In this example a savings of 4.375 would be realized. Of course, all these ideas would require special readers programs, instead of just opening up the files in a text editor -- a small price to pay for such a large space savings.

But in case only sequential reading of the primes is desired, there is an even more efficient representation: just write out the prime gaps, i.e. the difference between consecutive primes. One bit can be saved by realizing that, other than 3 - 2, all prime gaps are even, so one can just store the prime gap divided by two. One could really nitpick and point out that 0 is an impossible prime gap, so one could subtract two from all gaps, but this gain is too small to be worth the hassle. According to Thomas Nicely (incidentally the guy who discovered the bug on the Pentium FPU during his studies on computational number theory), one could store the prime gaps of all

Re:I just RTFA...

[coyote-san]

You overlooked another key optimization - no number can have prime factors larger than sqrt(N). You don't need to bother computing any of the intermediate primes. With a 210-into-48 encoding you can quickly determine primality of any 32-bit number with a precomputed block of only 2k! This is such a small buffer that you could easily bootstrap it the first time the function is called - you can bootstrap the process by hardcoding nothing more than the first 6 bytes.

This raises an interesting question - is it faster to load precomputed blocks from disk or to regenerate it in the L1 cache? One buffer would hold the first block, a second buffer would hold the 2nd-Nth block computed on demand, and a third buffer would hold the target block that's being updated by each successive block.

Re: I just RTFA...

[Omniscient Ferret]

That's where you represent primality or lack of it by setting and clearing bits of an integer.
Working out the mask for those bit operations can take enough time to cancel out the space benefits. I think you could optimize the different masks & their steppings.

... there is an even more efficient representation: just write out the prime gaps, i.e. the difference between consecutive primes.
I tried this, but I don't know where that program is now. I used 0 as an escape code, allowing me to either note the next delta in hex, or else the current number, to detect errors. I only did this for output, though.

I could see that being effective for the trial factor prime number table, too.

Oh, and the nfsnet thing that's your homepage is cool.

There's a primes program in (uh, on Debian) the bsdgames package; if memory serves, it does the wheel thing suggested.

A suggestion for the article:
If the trial factor prime number table gets big enough to swap (either out of cache, or out to disk), just segment the trial factors too.

This is something I've given a lot of thought to; I've been interested in prime numbers since I was a kid. I've thought about making CDs or DVDs full of primes.

Re:I just RTFA...

[acidblood]

The GMP library uses a probabilistic algorithm as documented here. This is not as bad as it sounds. For instance, Dan Bernstein suggests in this paper a probabilistic procedure which consists of some trial factoring, a strong probable prime test and a Lucas test. Quoting him:

There is no known example of a composite number p that slips past all of these tests; I'm confident that nobody will ever find one by accident.

Of course, if you absolutely have to be sure that you're dealing with a prime, you could use a probabilistic algorithm to sieve for probable primes, then apply the n-1 test to obtain a primality certificate for your probable prime (which then becomes a certified prime with 0% probability of being composite). Of course this would be a bit costly, but for small integers the n-1 test is hard to beat -- ECPP and APR-CL have just too much overhead at this range. But I'm pretty sure your application, whatever it is, can live with the small probability that a PrP is composite.

Re:I just RTFA...

[johntromp]

Covering 30 numbers per byte and precomputing all the bit-indices, one arrives at: http://www.cwi.nl/~tromp/pearls.html#sieve

Re:Mersenne

[acidblood]

I'd venture to say I know what I'm talking about. One can't publish a paper about large-scale distributed primality proofs without knowing these basics. Here's a sample of my work. Now don't get me wrong, I'm not trying to sound like a smart-ass, but I'm not going to let someone question my knowledge in this field (which I have studied for years) and remain quiet.

Now it seems people here are not familiar with the terminology associated with Mersennes, so I'll give a crash course.

A Mersenne number M_n is of the form 2^n - 1. If 2^n - 1 is prime, then M_n is a Mersenne prime. If 2^n - 1 is composite, then M_n is a Mersenne composite. The integer n in this expression is called the exponent. When I talk about `exponents of Mersenne composites', I'm refering to some prime n for which M_n = 2^n - 1 is composite.

Now if 2^n - 1 is prime, then n must be prime, but the converse is not always true, the smallest counterexample being 2^11 - 1 = 2047 = 23*89. Turns out that for prime exponents n = 24,036,583, only in 41 cases is 2^n - 1 actually prime. A list of these cases can be found here. For the other 1,509,222 prime values of n = 24,036,583, 2^n - 1 is composite. This translates to a compositeness ratio of 99.9973%.

I hope everything is clearer now.

How to prove that all odd numbers are prime

[dargaud]

"It was mentioned on CNN that the new prime number discovered recently is four times bigger than the previous record." --John Blasik

"You know what seems odd to me? Numbers that aren't divisible by two." --Michael Wolf.

"I don't get even, I get odder."

How to prove that all odd numbers are prime? Well, this problem has different solutions whether you are a:

• Mathematician: 1 is prime, 3 is prime, 5 is prime, 7 is prime, and by induction we have that all the odd integers are prime.
• Physicist: 1 is prime, 3 is prime, 5 is prime, 7 is prime, 9 is an experimental error...
• Engineer: 1 is prime, 3 is prime, 5 is prime, 7 is prime, 9 is prime...
• Chemist: 1 prime, 3 prime, 5 prime... hey, let's publish!
• Modern physicist using renormalization: 1 is prime, 3 is prime, 5 is prime, 7 is prime, 9 is ... 9/3 is prime, 11 is prime, 13 is prime, 15 is ... 15/3 is prime, 17 is prime, 19 is prime, 21 is ... 21/3 is prime...
• Quantum Physicist: All numbers are equally prime and non-prime until observed.
• Professor: 1 is prime, 3 is prime, 5 is prime, 7 is prime, and the rest are left as an exercise for the student.
• Confused Undergraduate: Let p be any prime number larger than 2. Then p is not divisible by 2, so p is odd. QED
• Measure nontheorist: There are exactly as many odd numbers as primes (Euclid, Cantor), and exactly one even prime (namely 2), so there must be exactly one odd nonprime (namely 1).
• Cosmologist: 1 is prime, yes it is true....
• Computer Scientist: 1 is prime, 10 is prime, 11 is prime, 101 is prime...
• Programmer: 1 is prime, 3 is prime, 5 is prime, 7 is prime, 9 will be fixed in the next release, ...
• C programmer: 01 is prime, 03 is prime, 05 is prime, 07 is prime, 09 is really 011 which everyone knows is prime, ...
• BASIC programmer: What's a prime?
• COBOL programmer: What's an odd number?
• Windows programmer: 1 is prime. Wait...
• Mac programmer: Now why would anyone want to know about that? That's not user friendly. You don't worry about it, we'll take care of it for you.
• Bill Gates: 1. No one will ever need any more than 1.
• ZX-81 Computer Programmer: 1 is prime, 3 is prime, Out of Memory.
• Pentium owner: 1 is prime, 3 is prime, 5 is prime, 7 is prime, 8.9999978 is prime...
• GNU programmer: % prime
  usage: prime [-nV] [--quiet] [--silent] [--version] [-e script] --catenate --concatenate | c --create | d --diff --compare | r --append | t --list | u --update | x -extract --get [ --atime-preserve ] [ -b, --block-size N ] [ -B, --read-full-blocks ] [ -C, --directory DIR ] [--checkpoint ] [ -f, --file [HOSTNAME:]F ] [ --force-local ] [ -F, --info-script F --new-volume-script F ] [-G, --incremental ] [ -g, --listed-incremental F ] [ -h, --dereference ] [ -i, --ignore-zeros ] [ --ignore-failed-read ] [ -k, --keep-old-files ] [ -K, --starting-file F ] [ -l, --one-file-system ] [ -L, --tape-length N ] [ -m, --modification-time ] [ -M, --multi-volume ] [ -N, --after-date DATE, --newer DATE ] [ -o, --old-archive, --portability ] [ -O, --to-stdout ] [ -p, --same-permissions, --preserve-permissions ] [ -P, --absolute-paths ] [ --preserve ] [ -R, --record-number ] [ [-f script-file] [--expression=script] [--file=script-file] [file...]
  prime: you must specify exactly one of the r, c, t, x, or d options
  For more information, type "prime --help''
• Unix programmer: 1 is prime, 3 is prime, 5 is prime, 7 is prime, ...
  Segmentation fault, Core dumped.
• Computer programmer: 1 is prime, 3 is prime, 5 is prime

Re:This begs the question

[bertas28]

This somehow reminds me of a series of memorandum I came across in Paradox, a magazine published by the Melbourne Uni Math Society (MUMS).

It consists of a bunch of proofs that there is no largest prime - the list of proofs is entitled "15 good reasons why Pure Mathematics is not taught to first year students." My favourites are:
Proof by example:
"Let x be the largest prime. Then x=91 but 91+6=97 which is prime. Therefore 91 cannot be the largest prime number. Therefore there is no largest prime."

Proof by intuition:
"Prime numbers are integers that can be divided by themselves only; prime numbers are odd with the exception of 2. By intuition as n->infinity, there will always be an odd number cannot be divided by another number besides itself."

Proof by experimental data:
"Suppose n is the highest prime. Then 2n-1 is also prime. But 2n-1>n so there is no highest prime. (Check: 2*2-1=3, 2*3-1=5, 2*5-1=11, 2*11-1=23, so true.)"

But the greatest of them all:
Proof by having no idea what a prime is:
"Say the largest prime is x, then 2x is also a prime since the statement is true for all natural numbers."

Go to http://www.ms.unimelb.edu.au/~paradox/archive/ for the magazine, these proofs appeared in issue 2 of 2002.