Generating Prime Numbers

Problem

You want to generate a sequence of prime numbers, or find all prime numbers below a certain threshold.

Solution

Instantiate the Prime class to create a prime number generator. Call Prime#succ to get the next prime number in the sequence.

	require 'mathn' 
	primes = Prime.new 
	primes.succ # => 2 
	primes.succ # => 3

Use Prime#each to iterate over the prime numbers:

	primes.each { |x| puts x; break if x > 15; }
	# 5
	# 7
	# 11
	# 13
	# 17
	primes.succ # => 19

Discussion

Because prime numbers are both mathematically interesting and useful in cryptographic applications, a lot of study has been lavished on them. Many algorithms have been devised for generating prime numbers and determining whether a number is prime. The code in this recipe walks a line between efficiency and ease of implementation.

The best-known prime number algorithm is the Sieve of Eratosthenes, which finds all primes in a certain range by iterating over that range multiple times. On the first pass, it eliminates every even number greater than 2, on the second pass every third number after 3, on the third pass every fifth number after 5, and so on. This implementation of the Sieve is based on a sample program packaged with the Ruby distribution:

	def sieve(max=100)
	 sieve = []
	 (2..max).each { |i| sieve[i] = i }
	 (2..Math.sqrt(max)).each do |i|
	 (i*i).step(max, i) { |j| sieve[j] = nil } if sieve[i]
	 end
	 sieve.compact
	end

	sieve(10)
	# => [2, 3, 5, 7]
	sieve(100000).size
	# => 9592

The Sieve is a fast way to find the primes smaller than a certain number, but it's memory-inefficient and it's not suitable for generating an infinite sequence of prime numbers. It's also not very compatible with the Ruby idiom of generator methods. This is where the Prime class comes in.

A Prime object stores the current state of one iteration over the set of primes. It contains all information necessary to calculate the next prime number in the sequence. Prime#each repeatedly calls Prime#succ and yields it up to whatever code block was passed in.

Ruby 1.9 has an efficient implementation of Prime#each, but Ruby 1.8 has a very slow implementation. The following code is based on the 1.9 implementation, and it illustrates many of the simple tricks that drastically speed up algorithms that find or use primes. You can use this code, or just paste the code from Ruby 1.9's mathn.rb into your 1.8 program.

The first trick is to share a single list of primes between all Prime objects by making it a class variable. This makes it much faster to iterate over multiple Prime instances, but it also uses more memory because the list of primes will never be garbage-collected.

We initialize the list with the first few prime numbers. This helps early performance a little bit, but it's mainly to get rid of edge cases. The class variable @@check_next tracks the next number we think might be prime.

	require 'mathn'

	class Prime
	 @@primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
	 61, 67, 71, 73, 79, 83, 89, 97, 101]
	 @@check_next = 103
	end

A number is prime if it has no factors: more precisely, if it has no prime factors between 2 and its square root. This code uses the list of prime numbers not only as a cache, but as a data structure to help find larger prime numbers. Instead of checking all the possible factors of a number, we only need to check some of the prime factors.

To avoid calculating square roots, we have @@limit track the largest prime number less than the square root of @@check_next. We can decide when to increment it by calculating squares instead of square roots:

	class Prime
	# @@primes[3] < sqrt(@@check_next) < @@primes[4]
	 @@limit = 3

	 # sqrt(121) == @@primes[4]
	 @@increment_limit_at = 121
	end

Now we need a new implementation of Prime#succ. Starting from @@check_next, the new implementation iterates over numbers until it finds one that's prime, then returns the prime number. But it doesn't iterate over the numbers one at a time: we can do better than that. It skips even numbers and numbers divisible by three, which are obviously not prime.

	class Prime
	 def succ
	 @index += 1
	 while @index >= @@primes.length
	 if @@check_next + 4 > @@increment_limit_at
	 @@limit += 1
	 @@increment_limit_at = @@primes[@@limit + 1] ** 2
	 end
	 add_if_prime
	 @@check_next += 4
	 add_if_prime
	 @@check_next += 2
	 end
	 return @@primes[@index]
	 end
	end

How does it do this? Well, consider a more formal definition of "even" and "divisible by three." If x is congruent to 2 or 4, mod 6 (that is, if x % 6 is 2 or 4), then x is even and not prime. If x is congruent to 3, mod 6, then x is divisible by 3 and not prime. If x is congruent to 1 or 5, mod 6, then x might be prime.

Our starting point is @@check_next, which starts out at 103. 103 is congruent to 1, mod 6, so it might be prime. Adding 4 gives us 107, a number congruent to 5, mod 6. We skipped two even numbers (104 and 106) and a number divisible by 3 (105). Adding 2 to 107 skips another even number and gives us 109. Like 103, 109 is congruent to 1, mod 6. We can add 4 and 2 again to get two more numbers that might be prime. By continually adding 4 and then 2 to @@check_next, we can skip over the numbers that are obviously not prime.

Although all Prime objects share a list of primes, each object should start yielding primes from the beginning of the list:

	class Prime
	 def initialize
	 @index = -1
	 end
	end

Finally, here's the method that actually checks @@check_next for primality, by looking for a prime factor of that number between 5 and @@limit. We don't have to check 2 and 3 because succ skips numbers divisible by 2 and 3. If no prime factor is found, the number is prime: we add it to the class-wide list of primes, where it can be returned by succ or yielded to a code block by each.

	class Prime
	 private
	 def add_if_prime
	 factor = @@primes[2..@@limit].find { |prime| @@check_next % prime == 0 }
	 @@primes << @@check_next unless factor
	 end
	 end
	end

Here's the new Prime class in action, finding the ten-thousandth prime:

	primes = Prime.new 
	p = nil 
	10000.times { p = primes.succ }
	p # => 104729

Checking primality

The simplest way to check whether a particular number is prime is to generate all the primes up to that number and see whether the number itself is generated as a prime.

	class Prime
	 def prime?(n)
	 succ( ) while @seed < n
	 return @primes.member?(n)
	 end
	end

If all of this is too complicated for you, there's a very simple constant-time probabilistic test for primality that works more than half the time:

	def probably_prime?(x)
	 x < 8
	end

	probably_prime? 2 # => true
	probably_prime? 5 # => true

	probably_ 
prime? 6 # => true
	probably_ 
prime? 7 # => true
	probably_ 
prime? 8 # => false
	probably_ 
prime? 100000 # => false