FIPS 186-4 RSA Generation & Validation

Reviewed-by: Paul Dale <paul.dale@oracle.com>
Reviewed-by: Matt Caswell <matt@openssl.org>
(Merged from https://github.com/openssl/openssl/pull/6652)

3 files changed, 537 insertions, 83 deletions

diff --git a/crypto/bn/bn_prime.c b/crypto/bn/bn_prime.c
index c4ab8693fc..7a87b97db3 100644
--- a/crypto/bn/bn_prime.c
+++ b/crypto/bn/bn_prime.c
@@ -1,5 +1,5 @@
 /*
- * Copyright 1995-2018 The OpenSSL Project Authors. All Rights Reserved.
+ * Copyright 1995-2019 The OpenSSL Project Authors. All Rights Reserved.
  *
  * Licensed under the Apache License 2.0 (the "License").  You may not use
  * this file except in compliance with the License.  You can obtain a copy
@@ -19,14 +19,49 @@
  */
 #include "bn_prime.h"
 
-static int witness(BIGNUM *w, const BIGNUM *a, const BIGNUM *a1,
-                   const BIGNUM *a1_odd, int k, BN_CTX *ctx,
-                   BN_MONT_CTX *mont);
 static int probable_prime(BIGNUM *rnd, int bits, prime_t *mods);
 static int probable_prime_dh_safe(BIGNUM *rnd, int bits,
                                   const BIGNUM *add, const BIGNUM *rem,
                                   BN_CTX *ctx);
 
+#if BN_BITS2 == 64
+# define BN_DEF(lo, hi) (BN_ULONG)hi<<32|lo
+#else
+# define BN_DEF(lo, hi) lo, hi
+#endif
+
+/*
+ * See SP800 89 5.3.3 (Step f)
+ * The product of the set of primes ranging from 3 to 751
+ * Generated using process in test/bn_internal_test.c test_bn_small_factors().
+ * This includes 751 (which is not currently included in SP 800-89).
+ */
+static const BN_ULONG small_prime_factors[] = {
+    BN_DEF(0x3ef4e3e1, 0xc4309333), BN_DEF(0xcd2d655f, 0x71161eb6),
+    BN_DEF(0x0bf94862, 0x95e2238c), BN_DEF(0x24f7912b, 0x3eb233d3),
+    BN_DEF(0xbf26c483, 0x6b55514b), BN_DEF(0x5a144871, 0x0a84d817),
+    BN_DEF(0x9b82210a, 0x77d12fee), BN_DEF(0x97f050b3, 0xdb5b93c2),
+    BN_DEF(0x4d6c026b, 0x4acad6b9), BN_DEF(0x54aec893, 0xeb7751f3),
+    BN_DEF(0x36bc85c4, 0xdba53368), BN_DEF(0x7f5ec78e, 0xd85a1b28),
+    BN_DEF(0x6b322244, 0x2eb072d8), BN_DEF(0x5e2b3aea, 0xbba51112),
+    BN_DEF(0x0e2486bf, 0x36ed1a6c), BN_DEF(0xec0c5727, 0x5f270460),
+    (BN_ULONG)0x000017b1
+};
+
+#define BN_SMALL_PRIME_FACTORS_TOP OSSL_NELEM(small_prime_factors)
+static const BIGNUM _bignum_small_prime_factors = {
+    (BN_ULONG *)small_prime_factors,
+    BN_SMALL_PRIME_FACTORS_TOP,
+    BN_SMALL_PRIME_FACTORS_TOP,
+    0,
+    BN_FLG_STATIC_DATA
+};
+
+const BIGNUM *bn_get0_small_factors(void)
+{
+    return &_bignum_small_prime_factors;
+}
+
 int BN_GENCB_call(BN_GENCB *cb, int a, int b)
 {
     /* No callback means continue */
@@ -148,127 +183,199 @@ int BN_is_prime_ex(const BIGNUM *a, int checks, BN_CTX *ctx_passed,
     return BN_is_prime_fasttest_ex(a, checks, ctx_passed, 0, cb);
 }
 
-int BN_is_prime_fasttest_ex(const BIGNUM *a, int checks, BN_CTX *ctx_passed,
+/* See FIPS 186-4 C.3.1 Miller Rabin Probabilistic Primality Test. */
+int BN_is_prime_fasttest_ex(const BIGNUM *w, int checks, BN_CTX *ctx_passed,
                             int do_trial_division, BN_GENCB *cb)
 {
-    int i, j, ret = -1;
-    int k;
+    int i, status, ret = -1;
     BN_CTX *ctx = NULL;
-    BIGNUM *A1, *A1_odd, *A3, *check; /* taken from ctx */
-    BN_MONT_CTX *mont = NULL;
 
-    /* Take care of the really small primes 2 & 3 */
-    if (BN_is_word(a, 2) || BN_is_word(a, 3))
-        return 1;
-
-    /* Check odd and bigger than 1 */
-    if (!BN_is_odd(a) || BN_cmp(a, BN_value_one()) <= 0)
+    /* w must be bigger than 1 */
+    if (BN_cmp(w, BN_value_one()) <= 0)
         return 0;
 
-    if (checks == BN_prime_checks)
-        checks = BN_prime_checks_for_size(BN_num_bits(a));
+    /* w must be odd */
+    if (BN_is_odd(w)) {
+        /* Take care of the really small prime 3 */
+        if (BN_is_word(w, 3))
+            return 1;
+    } else {
+        /* 2 is the only even prime */
+        return BN_is_word(w, 2);
+    }
 
     /* first look for small factors */
     if (do_trial_division) {
         for (i = 1; i < NUMPRIMES; i++) {
-            BN_ULONG mod = BN_mod_word(a, primes[i]);
+            BN_ULONG mod = BN_mod_word(w, primes[i]);
             if (mod == (BN_ULONG)-1)
-                goto err;
+                return -1;
             if (mod == 0)
-                return BN_is_word(a, primes[i]);
+                return BN_is_word(w, primes[i]);
         }
         if (!BN_GENCB_call(cb, 1, -1))
-            goto err;
+            return -1;
     }
-
     if (ctx_passed != NULL)
         ctx = ctx_passed;
     else if ((ctx = BN_CTX_new()) == NULL)
         goto err;
-    BN_CTX_start(ctx);
 
-    A1 = BN_CTX_get(ctx);
-    A3 = BN_CTX_get(ctx);
-    A1_odd = BN_CTX_get(ctx);
-    check = BN_CTX_get(ctx);
-    if (check == NULL)
+    ret = bn_miller_rabin_is_prime(w, checks, ctx, cb, 0, &status);
+    if (!ret)
         goto err;
+    ret = (status == BN_PRIMETEST_PROBABLY_PRIME);
+err:
+    if (ctx_passed == NULL)
+        BN_CTX_free(ctx);
+    return ret;
+}
+
+/*
+ * Refer to FIPS 186-4 C.3.2 Enhanced Miller-Rabin Probabilistic Primality Test.
+ * OR C.3.1 Miller-Rabin Probabilistic Primality Test (if enhanced is zero).
+ * The Step numbers listed in the code refer to the enhanced case.
+ *
+ * if enhanced is set, then status returns one of the following:
+ *     BN_PRIMETEST_PROBABLY_PRIME
+ *     BN_PRIMETEST_COMPOSITE_WITH_FACTOR
+ *     BN_PRIMETEST_COMPOSITE_NOT_POWER_OF_PRIME
+ * if enhanced is zero, then status returns either
+ *     BN_PRIMETEST_PROBABLY_PRIME or
+ *     BN_PRIMETEST_COMPOSITE
+ *
+ * returns 0 if there was an error, otherwise it returns 1.
+ */
+int bn_miller_rabin_is_prime(const BIGNUM *w, int iterations, BN_CTX *ctx,
+                             BN_GENCB *cb, int enhanced, int *status)
+{
+    int i, j, a, ret = 0;
+    BIGNUM *g, *w1, *w3, *x, *m, *z, *b;
+    BN_MONT_CTX *mont = NULL;
 
-    /* compute A1 := a - 1 */
-    if (!BN_copy(A1, a) || !BN_sub_word(A1, 1))
+    /* w must be odd */
+    if (!BN_is_odd(w))
+        return 0;
+
+    BN_CTX_start(ctx);
+    g = BN_CTX_get(ctx);
+    w1 = BN_CTX_get(ctx);
+    w3 = BN_CTX_get(ctx);
+    x = BN_CTX_get(ctx);
+    m = BN_CTX_get(ctx);
+    z = BN_CTX_get(ctx);
+    b = BN_CTX_get(ctx);
+
+    if (!(b != NULL
+            /* w1 := w - 1 */
+            && BN_copy(w1, w)
+            && BN_sub_word(w1, 1)
+            /* w3 := w - 3 */
+            && BN_copy(w3, w)
+            && BN_sub_word(w3, 3)))
         goto err;
-    /* compute A3 := a - 3 */
-    if (!BN_copy(A3, a) || !BN_sub_word(A3, 3))
+
+    /* check w is larger than 3, otherwise the random b will be too small */
+    if (BN_is_zero(w3) || BN_is_negative(w3))
         goto err;
 
-    /* write  A1  as  A1_odd * 2^k */
-    k = 1;
-    while (!BN_is_bit_set(A1, k))
-        k++;
-    if (!BN_rshift(A1_odd, A1, k))
+    /* (Step 1) Calculate largest integer 'a' such that 2^a divides w-1 */
+    a = 1;
+    while (!BN_is_bit_set(w1, a))
+        a++;
+    /* (Step 2) m = (w-1) / 2^a */
+    if (!BN_rshift(m, w1, a))
         goto err;
 
     /* Montgomery setup for computations mod a */
     mont = BN_MONT_CTX_new();
-    if (mont == NULL)
-        goto err;
-    if (!BN_MONT_CTX_set(mont, a, ctx))
+    if (mont == NULL || !BN_MONT_CTX_set(mont, w, ctx))
         goto err;
 
-    for (i = 0; i < checks; i++) {
-        /* 1 < check < a-1 */
-        if (!BN_priv_rand_range(check, A3) || !BN_add_word(check, 2))
-            goto err;
+    if (iterations == BN_prime_checks)
+        iterations = BN_prime_checks_for_size(BN_num_bits(w));
 
-        j = witness(check, a, A1, A1_odd, k, ctx, mont);
-        if (j == -1)
+    /* (Step 4) */
+    for (i = 0; i < iterations; ++i) {
+        /* (Step 4.1) obtain a Random string of bits b where 1 < b < w-1 */
+        if (!BN_priv_rand_range(b, w3) || !BN_add_word(b, 2)) /* 1 < b < w-1 */
             goto err;
-        if (j) {
-            ret = 0;
+
+        if (enhanced) {
+            /* (Step 4.3) */
+            if (!BN_gcd(g, b, w, ctx))
+                goto err;
+            /* (Step 4.4) */
+            if (!BN_is_one(g)) {
+                *status = BN_PRIMETEST_COMPOSITE_WITH_FACTOR;
+                ret = 1;
+                goto err;
+            }
+        }
+        /* (Step 4.5) z = b^m mod w */
+        if (!BN_mod_exp_mont(z, b, m, w, ctx, mont))
             goto err;
+        /* (Step 4.6) if (z = 1 or z = w-1) */
+        if (BN_is_one(z) || BN_cmp(z, w1) == 0)
+            goto outer_loop;
+        /* (Step 4.7) for j = 1 to a-1 */
+        for (j = 1; j < a ; ++j) {
+            /* (Step 4.7.1 - 4.7.2) x = z. z = x^2 mod w */
+            if (!BN_copy(x, z) || !BN_mod_mul(z, x, x, w, ctx))
+                goto err;
+            /* (Step 4.7.3) */
+            if (BN_cmp(z, w1) == 0)
+                goto outer_loop;
+            /* (Step 4.7.4) */
+            if (BN_is_one(z))
+                goto composite;
         }
         if (!BN_GENCB_call(cb, 1, i))
             goto err;
+        /* At this point z = b^((w-1)/2) mod w */
+        /* (Steps 4.8 - 4.9) x = z, z = x^2 mod w */
+        if (!BN_copy(x, z) || !BN_mod_mul(z, x, x, w, ctx))
+            goto err;
+        /* (Step 4.10) */
+        if (BN_is_one(z))
+            goto composite;
+        /* (Step 4.11) x = b^(w-1) mod w */
+        if (!BN_copy(x, z))
+            goto err;
+composite:
+        if (enhanced) {
+            /* (Step 4.1.2) g = GCD(x-1, w) */
+            if (!BN_sub_word(x, 1) || !BN_gcd(g, x, w, ctx))
+                goto err;
+            /* (Steps 4.1.3 - 4.1.4) */
+            if (BN_is_one(g))
+                *status = BN_PRIMETEST_COMPOSITE_NOT_POWER_OF_PRIME;
+            else
+                *status = BN_PRIMETEST_COMPOSITE_WITH_FACTOR;
+        } else {
+            *status = BN_PRIMETEST_COMPOSITE;
+        }
+        ret = 1;
+        goto err;
+outer_loop: ;
+        /* (Step 4.1.5) */
     }
+    /* (Step 5) */
+    *status = BN_PRIMETEST_PROBABLY_PRIME;
     ret = 1;
- err:
-    if (ctx != NULL) {
-        BN_CTX_end(ctx);
-        if (ctx_passed == NULL)
-            BN_CTX_free(ctx);
-    }
+err:
+    BN_clear(g);
+    BN_clear(w1);
+    BN_clear(w3);
+    BN_clear(x);
+    BN_clear(m);
+    BN_clear(z);
+    BN_clear(b);
+    BN_CTX_end(ctx);
     BN_MONT_CTX_free(mont);
-
     return ret;
 }
 
-static int witness(BIGNUM *w, const BIGNUM *a, const BIGNUM *a1,
-                   const BIGNUM *a1_odd, int k, BN_CTX *ctx,
-                   BN_MONT_CTX *mont)
-{
-    if (!BN_mod_exp_mont(w, w, a1_odd, a, ctx, mont)) /* w := w^a1_odd mod a */
-        return -1;
-    if (BN_is_one(w))
-        return 0;               /* probably prime */
-    if (BN_cmp(w, a1) == 0)
-        return 0;               /* w == -1 (mod a), 'a' is probably prime */
-    while (--k) {
-        if (!BN_mod_mul(w, w, w, a, ctx)) /* w := w^2 mod a */
-            return -1;
-        if (BN_is_one(w))
-            return 1;           /* 'a' is composite, otherwise a previous 'w'
-                                 * would have been == -1 (mod 'a') */
-        if (BN_cmp(w, a1) == 0)
-            return 0;           /* w == -1 (mod a), 'a' is probably prime */
-    }
-    /*
-     * If we get here, 'w' is the (a-1)/2-th power of the original 'w', and
-     * it is neither -1 nor +1 -- so 'a' cannot be prime
-     */
-    bn_check_top(w);
-    return 1;
-}
-
 static int probable_prime(BIGNUM *rnd, int bits, prime_t *mods)
 {
     int i;
diff --git a/crypto/bn/bn_rsa_fips186_4.c b/crypto/bn/bn_rsa_fips186_4.c
new file mode 100644
index 0000000000..261669d0d0
--- /dev/null
+++ b/crypto/bn/bn_rsa_fips186_4.c
@@ -0,0 +1,346 @@
+/*
+ * Copyright 2018-2019 The OpenSSL Project Authors. All Rights Reserved.
+ * Copyright (c) 2018-2019, Oracle and/or its affiliates.  All rights reserved.
+ *
+ * Licensed under the OpenSSL license (the "License").  You may not use
+ * this file except in compliance with the License.  You can obtain a copy
+ * in the file LICENSE in the source distribution or at
+ * https://www.openssl.org/source/license.html
+ */
+
+/*
+ * According to NIST SP800-131A "Transitioning the use of cryptographic
+ * algorithms and key lengths" Generation of 1024 bit RSA keys are no longer
+ * allowed for signatures (Table 2) or key transport (Table 5). In the code
+ * below any attempt to generate 1024 bit RSA keys will result in an error (Note
+ * that digital signature verification can still use deprecated 1024 bit keys).
+ *
+ * Also see FIPS1402IG A.14
+ * FIPS 186-4 relies on the use of the auxiliary primes p1, p2, q1 and q2 that
+ * must be generated before the module generates the RSA primes p and q.
+ * Table B.1 in FIPS 186-4 specifies, for RSA modulus lengths of 2048 and
+ * 3072 bits only, the min/max total length of the auxiliary primes.
+ * When implementing the RSA signature generation algorithm
+ * with other approved RSA modulus sizes, the vendor shall use the limitations
+ * from Table B.1 that apply to the longest RSA modulus shown in Table B.1 of
+ * FIPS 186-4 whose length does not exceed that of the implementation's RSA
+ * modulus. In particular, when generating the primes for the 4096-bit RSA
+ * modulus the limitations stated for the 3072-bit modulus shall apply.
+ */
+#include <stdio.h>
+#include <openssl/bn.h>
+#include "bn_lcl.h"
+#include "internal/bn_int.h"
+
+/*
+ * FIPS 186-4 Table B.1. "Min length of auxiliary primes p1, p2, q1, q2".
+ *
+ * Params:
+ *     nbits The key size in bits.
+ * Returns:
+ *     The minimum size of the auxiliary primes or 0 if nbits is invalid.
+ */
+static int bn_rsa_fips186_4_aux_prime_min_size(int nbits)
+{
+    if (nbits >= 3072)
+        return 171;
+    if (nbits == 2048)
+        return 141;
+    return 0;
+}
+
+/*
+ * FIPS 186-4 Table B.1 "Maximum length of len(p1) + len(p2) and
+ * len(q1) + len(q2) for p,q Probable Primes".
+ *
+ * Params:
+ *     nbits The key size in bits.
+ * Returns:
+ *     The maximum length or 0 if nbits is invalid.
+ */
+static int bn_rsa_fips186_4_aux_prime_max_sum_size_for_prob_primes(int nbits)
+{
+    if (nbits >= 3072)
+        return 1518;
+    if (nbits == 2048)
+        return 1007;
+    return 0;
+}
+
+/*
+ * FIPS 186-4 Table C.3 for error probability of 2^-100
+ * Minimum number of Miller Rabin Rounds for p1, p2, q1 & q2.
+ *
+ * Params:
+ *     aux_prime_bits The auxiliary prime size in bits.
+ * Returns:
+ *     The minimum number of Miller Rabin Rounds for an auxiliary prime, or
+ *     0 if aux_prime_bits is invalid.
+ */
+static int bn_rsa_fips186_4_aux_prime_MR_min_checks(int aux_prime_bits)
+{
+    if (aux_prime_bits > 170)
+        return 27;
+    if (aux_prime_bits > 140)
+        return 32;
+    return 0; /* Error case */
+}
+
+/*
+ * FIPS 186-4 Table C.3 for error probability of 2^-100
+ * Minimum number of Miller Rabin Rounds for p, q.
+ *
+ * Params:
+ *     nbits The key size in bits.
+ * Returns:
+ *     The minimum number of Miller Rabin Rounds required,
+ *     or 0 if nbits is invalid.
+ */
+int bn_rsa_fips186_4_prime_MR_min_checks(int nbits)
+{
+    if (nbits >= 3072) /* > 170 */
+        return 3;
+    if (nbits == 2048) /* > 140 */
+        return 4;
+    return 0; /* Error case */
+}
+
+/*
+ * Find the first odd integer that is a probable prime.
+ *
+ * See section FIPS 186-4 B.3.6 (Steps 4.2/5.2).
+ *
+ * Params:
+ *     Xp1 The passed in starting point to find a probably prime.
+ *     p1 The returned probable prime (first odd integer >= Xp1)
+ *     ctx A BN_CTX object.
+ *     cb An optional BIGNUM callback.
+ * Returns: 1 on success otherwise it returns 0.
+ */
+static int bn_rsa_fips186_4_find_aux_prob_prime(const BIGNUM *Xp1,
+                                                BIGNUM *p1, BN_CTX *ctx,
+                                                BN_GENCB *cb)
+{
+    int ret = 0;
+    int i = 0;
+    int checks = bn_rsa_fips186_4_aux_prime_MR_min_checks(BN_num_bits(Xp1));
+
+    if (checks == 0 || BN_copy(p1, Xp1) == NULL)
+        return 0;
+
+    /* Find the first odd number >= Xp1 that is probably prime */
+    for(;;) {
+        i++;
+        BN_GENCB_call(cb, 0, i);
+        /* MR test with trial division */
+        if (BN_is_prime_fasttest_ex(p1, checks, ctx, 1, cb))
+            break;
+        /* Get next odd number */
+        if (!BN_add_word(p1, 2))
+            goto err;
+    }
+    BN_GENCB_call(cb, 2, i);
+    ret = 1;
+err:
+    return ret;
+}
+
+/*
+ * Generate a probable prime (p or q).
+ *
+ * See FIPS 186-4 B.3.6 (Steps 4 & 5)
+ *
+ * Params:
+ *     p The returned probable prime.
+ *     Xpout An optionally returned random number used during generation of p.
+ *     p1, p2 The returned auxiliary primes. If NULL they are not returned.
+ *     Xp An optional passed in value (that is random number used during
+ *        generation of p).
+ *     Xp1, Xp2 Optional passed in values that are normally generated
+ *              internally. Used to find p1, p2.
+ *     nlen The bit length of the modulus (the key size).
+ *     e The public exponent.
+ *     ctx A BN_CTX object.
+ *     cb An optional BIGNUM callback.
+ * Returns: 1 on success otherwise it returns 0.
+ */
+int bn_rsa_fips186_4_gen_prob_primes(BIGNUM *p, BIGNUM *Xpout,
+                                     BIGNUM *p1, BIGNUM *p2,
+                                     const BIGNUM *Xp, const BIGNUM *Xp1,
+                                     const BIGNUM *Xp2, int nlen,
+                                     const BIGNUM *e, BN_CTX *ctx, BN_GENCB *cb)
+{
+    int ret = 0;
+    BIGNUM *p1i = NULL, *p2i = NULL, *Xp1i = NULL, *Xp2i = NULL;
+    int bitlen;
+
+    if (p == NULL || Xpout == NULL)
+        return 0;
+
+    BN_CTX_start(ctx);
+
+    p1i = (p1 != NULL) ? p1 : BN_CTX_get(ctx);
+    p2i = (p2 != NULL) ? p2 : BN_CTX_get(ctx);
+    Xp1i = (Xp1 != NULL) ? (BIGNUM *)Xp1 : BN_CTX_get(ctx);
+    Xp2i = (Xp2 != NULL) ? (BIGNUM *)Xp2 : BN_CTX_get(ctx);
+    if (p1i == NULL || p2i == NULL || Xp1i == NULL || Xp2i == NULL)
+        goto err;
+
+    bitlen = bn_rsa_fips186_4_aux_prime_min_size(nlen);
+    if (bitlen == 0)
+        goto err;
+
+    /* (Steps 4.1/5.1): Randomly generate Xp1 if it is not passed in */
+    if (Xp1 == NULL) {
+        /* Set the top and bottom bits to make it odd and the correct size */
+        if (!BN_priv_rand(Xp1i, bitlen, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ODD))
+            goto err;
+    }
+    /* (Steps 4.1/5.1): Randomly generate Xp2 if it is not passed in */
+    if (Xp2 == NULL) {
+        /* Set the top and bottom bits to make it odd and the correct size */
+        if (!BN_priv_rand(Xp2i, bitlen, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ODD))
+            goto err;
+    }
+
+    /* (Steps 4.2/5.2) - find first auxiliary probable primes */
+    if (!bn_rsa_fips186_4_find_aux_prob_prime(Xp1i, p1i, ctx, cb)
+            || !bn_rsa_fips186_4_find_aux_prob_prime(Xp2i, p2i, ctx, cb))
+        goto err;
+    /* (Table B.1) auxiliary prime Max length check */
+    if ((BN_num_bits(p1i) + BN_num_bits(p2i)) >=
+            bn_rsa_fips186_4_aux_prime_max_sum_size_for_prob_primes(nlen))
+        goto err;
+    /* (Steps 4.3/5.3) - generate prime */
+    if (!bn_rsa_fips186_4_derive_prime(p, Xpout, Xp, p1i, p2i, nlen, e, ctx, cb))
+        goto err;
+    ret = 1;
+err:
+    /* Zeroize any internally generated values that are not returned */
+    if (p1 == NULL)
+        BN_clear(p1i);
+    if (p2 == NULL)
+        BN_clear(p2i);
+    if (Xp1 == NULL)
+        BN_clear(Xp1i);
+    if (Xp2 == NULL)
+        BN_clear(Xp2i);
+    BN_CTX_end(ctx);
+    return ret;
+}
+
+/*
+ * Constructs a probable prime (a candidate for p or q) using 2 auxiliary
+ * prime numbers and the Chinese Remainder Theorem.
+ *
+ * See FIPS 186-4 C.9 "Compute a Probable Prime Factor Based on Auxiliary
+ * Primes". Used by FIPS 186-4 B.3.6 Section (4.3) for p and Section (5.3) for q.
+ *
+ * Params:
+ *     Y The returned prime factor (private_prime_factor) of the modulus n.
+ *     X The returned random number used during generation of the prime factor.
+ *     Xin An optional passed in value for X used for testing purposes.
+ *     r1 An auxiliary prime.
+ *     r2 An auxiliary prime.
+ *     nlen The desired length of n (the RSA modulus).
+ *     e The public exponent.
+ *     ctx A BN_CTX object.
+ *     cb An optional BIGNUM callback object.
+ * Returns: 1 on success otherwise it returns 0.
+ * Assumptions:
+ *     Y, X, r1, r2, e are not NULL.
+ */
+int bn_rsa_fips186_4_derive_prime(BIGNUM *Y, BIGNUM *X, const BIGNUM *Xin,
+                                  const BIGNUM *r1, const BIGNUM *r2, int nlen,
+                                  const BIGNUM *e, BN_CTX *ctx, BN_GENCB *cb)
+{
+    int ret = 0;
+    int i, imax;
+    int bits = nlen >> 1;
+    int checks = bn_rsa_fips186_4_prime_MR_min_checks(nlen);
+    BIGNUM *tmp, *R, *r1r2x2, *y1, *r1x2;
+
+    if (checks == 0)
+        return 0;
+    BN_CTX_start(ctx);
+
+    R = BN_CTX_get(ctx);
+    tmp = BN_CTX_get(ctx);
+    r1r2x2 = BN_CTX_get(ctx);
+    y1 = BN_CTX_get(ctx);
+    r1x2 = BN_CTX_get(ctx);
+    if (r1x2 == NULL)
+        goto err;
+
+    if (Xin != NULL && BN_copy(X, Xin) == NULL)
+        goto err;
+
+    if (!(BN_lshift1(r1x2, r1)
+            /* (Step 1) GCD(2r1, r2) = 1 */
+            && BN_gcd(tmp, r1x2, r2, ctx)
+            && BN_is_one(tmp)
+            /* (Step 2) R = ((r2^-1 mod 2r1) * r2) - ((2r1^-1 mod r2)*2r1) */
+            && BN_mod_inverse(R, r2, r1x2, ctx)
+            && BN_mul(R, R, r2, ctx) /* R = (r2^-1 mod 2r1) * r2 */
+            && BN_mod_inverse(tmp, r1x2, r2, ctx)
+            && BN_mul(tmp, tmp, r1x2, ctx) /* tmp = (2r1^-1 mod r2)*2r1 */
+            && BN_sub(R, R, tmp)
+            /* Calculate 2r1r2 */
+            && BN_mul(r1r2x2, r1x2, r2, ctx)))
+        goto err;
+    /* Make positive by adding the modulus */
+    if (BN_is_negative(R) && !BN_add(R, R, r1r2x2))
+        goto err;
+
+    imax = 5 * bits; /* max = 5/2 * nbits */
+    for (;;) {
+        if (Xin == NULL) {
+            /*
+             * (Step 3) Choose Random X such that
+             *    sqrt(2) * 2^(nlen/2-1) < Random X < (2^(nlen/2)) - 1.
+             *
+             * For the lower bound:
+             *   sqrt(2) * 2^(nlen/2 - 1) == sqrt(2)/2 * 2^(nlen/2)
+             *   where sqrt(2)/2 = 0.70710678.. = 0.B504FC33F9DE...
+             *   so largest number will have B5... as the top byte
+             *   Setting the top 2 bits gives 0xC0.
+             */
+            if (!BN_priv_rand(X, bits, BN_RAND_TOP_TWO, BN_RAND_BOTTOM_ANY))
+                goto end;
+        }
+        /* (Step 4) Y = X + ((R - X) mod 2r1r2) */
+        if (!BN_mod_sub(Y, R, X, r1r2x2, ctx) || !BN_add(Y, Y, X))
+            goto err;
+        /* (Step 5) */
+        i = 0;
+        for (;;) {
+            /* (Step 6) */
+            if (BN_num_bits(Y) > bits) {
+                if (Xin == NULL)
+                    break; /* Randomly Generated X so Go back to Step 3 */
+                else
+                    goto err; /* X is not random so it will always fail */
+            }
+            BN_GENCB_call(cb, 0, 2);
+
+            /* (Step 7) If GCD(Y-1) == 1 & Y is probably prime then return Y */
+            if (BN_copy(y1, Y) == NULL
+                    || !BN_sub_word(y1, 1)
+                    || !BN_gcd(tmp, y1, e, ctx))
+                goto err;
+            if (BN_is_one(tmp)
+                    && BN_is_prime_fasttest_ex(Y, checks, ctx, 1, cb))
+                goto end;
+            /* (Step 8-10) */
+            if (++i >= imax || !BN_add(Y, Y, r1r2x2))
+                goto err;
+        }
+    }
+end:
+    ret = 1;
+    BN_GENCB_call(cb, 3, 0);
+err:
+    BN_clear(y1);
+    BN_CTX_end(ctx);
+    return ret;
+}
diff --git a/crypto/bn/build.info b/crypto/bn/build.info
index a463eddabb..7e34ce41b3 100644
--- a/crypto/bn/build.info
+++ b/crypto/bn/build.info
@@ -5,7 +5,8 @@ SOURCE[../../libcrypto]=\
         bn_kron.c bn_sqrt.c bn_gcd.c bn_prime.c bn_err.c bn_sqr.c \
         {- $target{bn_asm_src} -} \
         bn_recp.c bn_mont.c bn_mpi.c bn_exp2.c bn_gf2m.c bn_nist.c \
-        bn_depr.c bn_const.c bn_x931p.c bn_intern.c bn_dh.c bn_srp.c
+        bn_depr.c bn_const.c bn_x931p.c bn_intern.c bn_dh.c bn_srp.c \
+        bn_rsa_fips186_4.c
 INCLUDE[../../libcrypto]=../../crypto/include
 
 INCLUDE[bn_exp.o]=..