7.5 Generating a Prime Number (Testing for Primality)

7.5.1 Problem

You need to generate a random prime number or test to see if a number is prime.

7.5.2 Solution

Use the routines provided by your arbitrary-precision math library, or generate a random odd number and use the Rabin-Miller primality test to see whether the number generated is actually prime.

7.5.3 Discussion

Good arbitrary-precision math libraries have functions that can automatically generate primes and determine to a near certainty whether a number is prime. In addition, these libraries should have functionality that produces "safe" primes (that is, a prime whose value minus 1 divided by 2 is also prime). You should also be able to ask for a prime that gives a particular remainder when you divide that prime by a particular number. The last two pieces of functionality are useful for generating parameters for Diffie-Hellman key exchange.

The OpenSSL functionality for generating and testing primes is discussed in Recipe 7.4.

The most common way primes are generated is by choosing a random odd number of the desired bit length from a secure pseudo-random source (we discuss pseudo-randomness in depth in Recipe 11.1). Generally, the output of the random number generator will have the first and last bits set. Setting the last bit ensures that the number is odd; no even numbers are primes. Setting the first bit ensures that the generated number really is of the desired bit length.

When generating RSA keys, people usually set the first two bits of all their potential primes. That way, if you multiply two primes of the same bit length together, they'll produce a result that's exactly twice the bit length. When people talk about the "bit length of an RSA key," they're generally talking about the size of such a product.

For determining whether a number is prime, most people use the Rabin-Miller test, which can determine primality with high probability. Every time you run the Rabin-Miller test and the test reports the number "may be prime," the actual probability of the number being prime increases dramatically. By the time you've run five iterations and have received "may be prime" every time, the odds of the random value's not being prime aren't worth worrying about.

If you are generating a prime number for use in Diffie-Hellman key exchange (i.e., a "safe" prime), you should test the extra conditions before you even check to see if the number itself is prime because doing so will speed up tests.

We provide the following code that implements Rabin-Miller on top of the OpenSSL BIGNUM library, which almost seems worthless, because if you're using OpenSSL, it already contains this test as an API function (again, see Recipe 7.4). However, the OpenSSL BIGNUM API is straightforward. It should be easy to take this code and translate it to work with whatever package you're using for arbitrary precision math.

Do note, though, that any library you use is likely already to have a function that performs this work for you.

In this code, we explicitly attempt division for the first 100 primes, although we recommend trying more primes than that. (OpenSSL itself tries 2,048, a widely recommended number.) We omit the additional primes for space reasons, but you can find a list of those primes on this book's web site. In addition, we use spc_rand( ) to get a random binary value. See Recipe 11.2 for a discussion of this function.

#include <stdlib.h>
#include <openssl/bn.h>
   
#define NUMBER_ITERS    5
#define NUMBER_PRIMES   100
   
static unsigned long primes[NUMBER_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, 103, 107, 109, 113, 127, 131, 
  137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 
  227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 
  313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 
  419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 
  509, 521, 523, 541
};
   
static int is_obviously_not_prime(BIGNUM *p);
static int passes_rabin_miller_once(BIGNUM *p);
static unsigned int calc_b_and_m(BIGNUM *p, BIGNUM *m);
   
int spc_is_probably_prime(BIGNUM *p) {
  int i;
  if (is_obviously_not_prime(p)) return 0;
  for (i = 0;  i < NUMBER_ITERS;  i++)
    if (!passes_rabin_miller_once(p))
      return 0;
  return 1;
}
BIGNUM *spc_generate_prime(int nbits) {
  BIGNUM        *p = BN_new(  );
  unsigned char binary_rep[nbits / 8];
   
  /* This code assumes we'll only ever want to generate primes with the number of
   * bits a multiple of eight!
   */
  if (nbits % 8 || !p) abort(  );
   
  for (;;) {
    spc_rand(binary_rep, nbits / 8);
  
    /* Set the two most significant and the least significant bits to 1. */
    binary_rep[0] |= 0xc0;
    binary_rep[nbits / 8 - 1] |= 1;
   
    /* Convert this number to its BIGNUM representation */
    if (!BN_bin2bn(binary_rep, nbits / 8, p)) abort(  );
   
    /* If you're going to test for suitability as a Diffie-Hellman prime, do so
     * before calling spc_is_probably_prime(p).
     */
    if (spc_is_probably_prime(p)) return p;
  }
}
   
/* Try simple division with all our small primes.  This is, for each prime, if it
 * evenly divides p, return 0.  Note that this obviously doesn't work if we're
 * checking a prime number that's in the list!
 */
static int is_obviously_not_prime(BIGNUM *p) {
  int i;
   
  for (i = 0;  i < NUMBER_PRIMES;  i++)
    if (!BN_mod_word(p, primes[i])) return 1;
  return 0;
}
   
static int passes_rabin_miller_once(BIGNUM *p) {
  BIGNUM       a, m, z, tmp;
  BN_CTX       *ctx;
  unsigned int b, i;
   
  /* Initialize a, m, z and tmp properly. */
  BN_init(&a);
  BN_init(&m);
  BN_init(&z);
  BN_init(&tmp);
   
  ctx = BN_CTX_new(  );
  b = calc_b_and_m(p, &m);
   
  /* a is a random number less than p: */
  if (!BN_rand_range(&a, p)) abort(  );
  
  /* z = a^m mod p. */
  if (!BN_mod_exp(&z, &a, &m, p, ctx)) abort(  );
   
  /* if z = 1 at the start, pass. */
  if (BN_is_one(&z)) return 1;
   
  for (i = 0;  i < b;  i++) {
    if (BN_is_one(&z)) return 0;
    
    /* if z = p-1, pass! */
    BN_copy(&tmp, &z);
    if (!BN_add_word(&tmp, 1)) abort(  );
    if (!BN_cmp(&tmp, p)) return 1;
   
    /* z = z^2 mod p */
    BN_mod_sqr(&tmp, &z, p, ctx);
    BN_copy(&z, &tmp);
  }
   
  /* if z = p-1, pass! */
  BN_copy(&tmp, &z);
  if (!BN_add_word(&tmp, 1)) abort(  );
  if (!BN_cmp(&tmp, p)) return 1;
 
  /* Fail! */ 
  return 0;
}
   
/* b = How many times does 2 divide p - 1?   This gets returned.
 * m is (p-1)/(2^b).
 */
static unsigned int calc_b_and_m(BIGNUM *p, BIGNUM *x) {
  unsigned int b;
   
  if (!BN_copy(x, p)) abort(  );
  if (!BN_sub_word(x, 1))  abort(  );
   
  for (b = 0;  !BN_is_odd(x);  b++)
    BN_div_word(x, 2);
  return b;
}

7.5.4 See Also

Recipe 7.4, Recipe 11.1, Recipe 11.2