10957b409SSimon J. Gerraty /*
20957b409SSimon J. Gerraty * Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
30957b409SSimon J. Gerraty *
40957b409SSimon J. Gerraty * Permission is hereby granted, free of charge, to any person obtaining
50957b409SSimon J. Gerraty * a copy of this software and associated documentation files (the
60957b409SSimon J. Gerraty * "Software"), to deal in the Software without restriction, including
70957b409SSimon J. Gerraty * without limitation the rights to use, copy, modify, merge, publish,
80957b409SSimon J. Gerraty * distribute, sublicense, and/or sell copies of the Software, and to
90957b409SSimon J. Gerraty * permit persons to whom the Software is furnished to do so, subject to
100957b409SSimon J. Gerraty * the following conditions:
110957b409SSimon J. Gerraty *
120957b409SSimon J. Gerraty * The above copyright notice and this permission notice shall be
130957b409SSimon J. Gerraty * included in all copies or substantial portions of the Software.
140957b409SSimon J. Gerraty *
150957b409SSimon J. Gerraty * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
160957b409SSimon J. Gerraty * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
170957b409SSimon J. Gerraty * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
180957b409SSimon J. Gerraty * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
190957b409SSimon J. Gerraty * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
200957b409SSimon J. Gerraty * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
210957b409SSimon J. Gerraty * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
220957b409SSimon J. Gerraty * SOFTWARE.
230957b409SSimon J. Gerraty */
240957b409SSimon J. Gerraty
250957b409SSimon J. Gerraty #include "inner.h"
260957b409SSimon J. Gerraty
270957b409SSimon J. Gerraty /*
280957b409SSimon J. Gerraty * Make a random integer of the provided size. The size is encoded.
290957b409SSimon J. Gerraty * The header word is untouched.
300957b409SSimon J. Gerraty */
310957b409SSimon J. Gerraty static void
mkrand(const br_prng_class ** rng,uint32_t * x,uint32_t esize)320957b409SSimon J. Gerraty mkrand(const br_prng_class **rng, uint32_t *x, uint32_t esize)
330957b409SSimon J. Gerraty {
340957b409SSimon J. Gerraty size_t u, len;
350957b409SSimon J. Gerraty unsigned m;
360957b409SSimon J. Gerraty
370957b409SSimon J. Gerraty len = (esize + 31) >> 5;
380957b409SSimon J. Gerraty (*rng)->generate(rng, x + 1, len * sizeof(uint32_t));
390957b409SSimon J. Gerraty for (u = 1; u < len; u ++) {
400957b409SSimon J. Gerraty x[u] &= 0x7FFFFFFF;
410957b409SSimon J. Gerraty }
420957b409SSimon J. Gerraty m = esize & 31;
430957b409SSimon J. Gerraty if (m == 0) {
440957b409SSimon J. Gerraty x[len] &= 0x7FFFFFFF;
450957b409SSimon J. Gerraty } else {
460957b409SSimon J. Gerraty x[len] &= 0x7FFFFFFF >> (31 - m);
470957b409SSimon J. Gerraty }
480957b409SSimon J. Gerraty }
490957b409SSimon J. Gerraty
500957b409SSimon J. Gerraty /*
510957b409SSimon J. Gerraty * This is the big-endian unsigned representation of the product of
520957b409SSimon J. Gerraty * all small primes from 13 to 1481.
530957b409SSimon J. Gerraty */
540957b409SSimon J. Gerraty static const unsigned char SMALL_PRIMES[] = {
550957b409SSimon J. Gerraty 0x2E, 0xAB, 0x92, 0xD1, 0x8B, 0x12, 0x47, 0x31, 0x54, 0x0A,
560957b409SSimon J. Gerraty 0x99, 0x5D, 0x25, 0x5E, 0xE2, 0x14, 0x96, 0x29, 0x1E, 0xB7,
570957b409SSimon J. Gerraty 0x78, 0x70, 0xCC, 0x1F, 0xA5, 0xAB, 0x8D, 0x72, 0x11, 0x37,
580957b409SSimon J. Gerraty 0xFB, 0xD8, 0x1E, 0x3F, 0x5B, 0x34, 0x30, 0x17, 0x8B, 0xE5,
590957b409SSimon J. Gerraty 0x26, 0x28, 0x23, 0xA1, 0x8A, 0xA4, 0x29, 0xEA, 0xFD, 0x9E,
600957b409SSimon J. Gerraty 0x39, 0x60, 0x8A, 0xF3, 0xB5, 0xA6, 0xEB, 0x3F, 0x02, 0xB6,
610957b409SSimon J. Gerraty 0x16, 0xC3, 0x96, 0x9D, 0x38, 0xB0, 0x7D, 0x82, 0x87, 0x0C,
620957b409SSimon J. Gerraty 0xF7, 0xBE, 0x24, 0xE5, 0x5F, 0x41, 0x04, 0x79, 0x76, 0x40,
630957b409SSimon J. Gerraty 0xE7, 0x00, 0x22, 0x7E, 0xB5, 0x85, 0x7F, 0x8D, 0x01, 0x50,
640957b409SSimon J. Gerraty 0xE9, 0xD3, 0x29, 0x42, 0x08, 0xB3, 0x51, 0x40, 0x7B, 0xD7,
650957b409SSimon J. Gerraty 0x8D, 0xCC, 0x10, 0x01, 0x64, 0x59, 0x28, 0xB6, 0x53, 0xF3,
660957b409SSimon J. Gerraty 0x50, 0x4E, 0xB1, 0xF2, 0x58, 0xCD, 0x6E, 0xF5, 0x56, 0x3E,
670957b409SSimon J. Gerraty 0x66, 0x2F, 0xD7, 0x07, 0x7F, 0x52, 0x4C, 0x13, 0x24, 0xDC,
680957b409SSimon J. Gerraty 0x8E, 0x8D, 0xCC, 0xED, 0x77, 0xC4, 0x21, 0xD2, 0xFD, 0x08,
690957b409SSimon J. Gerraty 0xEA, 0xD7, 0xC0, 0x5C, 0x13, 0x82, 0x81, 0x31, 0x2F, 0x2B,
700957b409SSimon J. Gerraty 0x08, 0xE4, 0x80, 0x04, 0x7A, 0x0C, 0x8A, 0x3C, 0xDC, 0x22,
710957b409SSimon J. Gerraty 0xE4, 0x5A, 0x7A, 0xB0, 0x12, 0x5E, 0x4A, 0x76, 0x94, 0x77,
720957b409SSimon J. Gerraty 0xC2, 0x0E, 0x92, 0xBA, 0x8A, 0xA0, 0x1F, 0x14, 0x51, 0x1E,
730957b409SSimon J. Gerraty 0x66, 0x6C, 0x38, 0x03, 0x6C, 0xC7, 0x4A, 0x4B, 0x70, 0x80,
740957b409SSimon J. Gerraty 0xAF, 0xCA, 0x84, 0x51, 0xD8, 0xD2, 0x26, 0x49, 0xF5, 0xA8,
750957b409SSimon J. Gerraty 0x5E, 0x35, 0x4B, 0xAC, 0xCE, 0x29, 0x92, 0x33, 0xB7, 0xA2,
760957b409SSimon J. Gerraty 0x69, 0x7D, 0x0C, 0xE0, 0x9C, 0xDB, 0x04, 0xD6, 0xB4, 0xBC,
770957b409SSimon J. Gerraty 0x39, 0xD7, 0x7F, 0x9E, 0x9D, 0x78, 0x38, 0x7F, 0x51, 0x54,
780957b409SSimon J. Gerraty 0x50, 0x8B, 0x9E, 0x9C, 0x03, 0x6C, 0xF5, 0x9D, 0x2C, 0x74,
790957b409SSimon J. Gerraty 0x57, 0xF0, 0x27, 0x2A, 0xC3, 0x47, 0xCA, 0xB9, 0xD7, 0x5C,
800957b409SSimon J. Gerraty 0xFF, 0xC2, 0xAC, 0x65, 0x4E, 0xBD
810957b409SSimon J. Gerraty };
820957b409SSimon J. Gerraty
830957b409SSimon J. Gerraty /*
840957b409SSimon J. Gerraty * We need temporary values for at least 7 integers of the same size
850957b409SSimon J. Gerraty * as a factor (including header word); more space helps with performance
860957b409SSimon J. Gerraty * (in modular exponentiations), but we much prefer to remain under
870957b409SSimon J. Gerraty * 2 kilobytes in total, to save stack space. The macro TEMPS below
880957b409SSimon J. Gerraty * exceeds 512 (which is a count in 32-bit words) when BR_MAX_RSA_SIZE
890957b409SSimon J. Gerraty * is greater than 4464 (default value is 4096, so the 2-kB limit is
900957b409SSimon J. Gerraty * maintained unless BR_MAX_RSA_SIZE was modified).
910957b409SSimon J. Gerraty */
920957b409SSimon J. Gerraty #define MAX(x, y) ((x) > (y) ? (x) : (y))
930957b409SSimon J. Gerraty #define ROUND2(x) ((((x) + 1) >> 1) << 1)
940957b409SSimon J. Gerraty
950957b409SSimon J. Gerraty #define TEMPS MAX(512, ROUND2(7 * ((((BR_MAX_RSA_SIZE + 1) >> 1) + 61) / 31)))
960957b409SSimon J. Gerraty
970957b409SSimon J. Gerraty /*
980957b409SSimon J. Gerraty * Perform trial division on a candidate prime. This computes
990957b409SSimon J. Gerraty * y = SMALL_PRIMES mod x, then tries to compute y/y mod x. The
1000957b409SSimon J. Gerraty * br_i31_moddiv() function will report an error if y is not invertible
1010957b409SSimon J. Gerraty * modulo x. Returned value is 1 on success (none of the small primes
1020957b409SSimon J. Gerraty * divides x), 0 on error (a non-trivial GCD is obtained).
1030957b409SSimon J. Gerraty *
1040957b409SSimon J. Gerraty * This function assumes that x is odd.
1050957b409SSimon J. Gerraty */
1060957b409SSimon J. Gerraty static uint32_t
trial_divisions(const uint32_t * x,uint32_t * t)1070957b409SSimon J. Gerraty trial_divisions(const uint32_t *x, uint32_t *t)
1080957b409SSimon J. Gerraty {
1090957b409SSimon J. Gerraty uint32_t *y;
1100957b409SSimon J. Gerraty uint32_t x0i;
1110957b409SSimon J. Gerraty
1120957b409SSimon J. Gerraty y = t;
1130957b409SSimon J. Gerraty t += 1 + ((x[0] + 31) >> 5);
1140957b409SSimon J. Gerraty x0i = br_i31_ninv31(x[1]);
1150957b409SSimon J. Gerraty br_i31_decode_reduce(y, SMALL_PRIMES, sizeof SMALL_PRIMES, x);
1160957b409SSimon J. Gerraty return br_i31_moddiv(y, y, x, x0i, t);
1170957b409SSimon J. Gerraty }
1180957b409SSimon J. Gerraty
1190957b409SSimon J. Gerraty /*
1200957b409SSimon J. Gerraty * Perform n rounds of Miller-Rabin on the candidate prime x. This
1210957b409SSimon J. Gerraty * function assumes that x = 3 mod 4.
1220957b409SSimon J. Gerraty *
1230957b409SSimon J. Gerraty * Returned value is 1 on success (all rounds completed successfully),
1240957b409SSimon J. Gerraty * 0 otherwise.
1250957b409SSimon J. Gerraty */
1260957b409SSimon J. Gerraty static uint32_t
miller_rabin(const br_prng_class ** rng,const uint32_t * x,int n,uint32_t * t,size_t tlen,br_i31_modpow_opt_type mp31)1270957b409SSimon J. Gerraty miller_rabin(const br_prng_class **rng, const uint32_t *x, int n,
1280957b409SSimon J. Gerraty uint32_t *t, size_t tlen, br_i31_modpow_opt_type mp31)
1290957b409SSimon J. Gerraty {
1300957b409SSimon J. Gerraty /*
1310957b409SSimon J. Gerraty * Since x = 3 mod 4, the Miller-Rabin test is simple:
1320957b409SSimon J. Gerraty * - get a random base a (such that 1 < a < x-1)
1330957b409SSimon J. Gerraty * - compute z = a^((x-1)/2) mod x
1340957b409SSimon J. Gerraty * - if z != 1 and z != x-1, the number x is composite
1350957b409SSimon J. Gerraty *
1360957b409SSimon J. Gerraty * We generate bases 'a' randomly with a size which is
1370957b409SSimon J. Gerraty * one bit less than x, which ensures that a < x-1. It
1380957b409SSimon J. Gerraty * is not useful to verify that a > 1 because the probability
1390957b409SSimon J. Gerraty * that we get a value a equal to 0 or 1 is much smaller
1400957b409SSimon J. Gerraty * than the probability of our Miller-Rabin tests not to
1410957b409SSimon J. Gerraty * detect a composite, which is already quite smaller than the
1420957b409SSimon J. Gerraty * probability of the hardware misbehaving and return a
1430957b409SSimon J. Gerraty * composite integer because of some glitch (e.g. bad RAM
1440957b409SSimon J. Gerraty * or ill-timed cosmic ray).
1450957b409SSimon J. Gerraty */
1460957b409SSimon J. Gerraty unsigned char *xm1d2;
1470957b409SSimon J. Gerraty size_t xlen, xm1d2_len, xm1d2_len_u32, u;
1480957b409SSimon J. Gerraty uint32_t asize;
1490957b409SSimon J. Gerraty unsigned cc;
1500957b409SSimon J. Gerraty uint32_t x0i;
1510957b409SSimon J. Gerraty
1520957b409SSimon J. Gerraty /*
1530957b409SSimon J. Gerraty * Compute (x-1)/2 (encoded).
1540957b409SSimon J. Gerraty */
1550957b409SSimon J. Gerraty xm1d2 = (unsigned char *)t;
1560957b409SSimon J. Gerraty xm1d2_len = ((x[0] - (x[0] >> 5)) + 7) >> 3;
1570957b409SSimon J. Gerraty br_i31_encode(xm1d2, xm1d2_len, x);
1580957b409SSimon J. Gerraty cc = 0;
1590957b409SSimon J. Gerraty for (u = 0; u < xm1d2_len; u ++) {
1600957b409SSimon J. Gerraty unsigned w;
1610957b409SSimon J. Gerraty
1620957b409SSimon J. Gerraty w = xm1d2[u];
1630957b409SSimon J. Gerraty xm1d2[u] = (unsigned char)((w >> 1) | cc);
1640957b409SSimon J. Gerraty cc = w << 7;
1650957b409SSimon J. Gerraty }
1660957b409SSimon J. Gerraty
1670957b409SSimon J. Gerraty /*
1680957b409SSimon J. Gerraty * We used some words of the provided buffer for (x-1)/2.
1690957b409SSimon J. Gerraty */
1700957b409SSimon J. Gerraty xm1d2_len_u32 = (xm1d2_len + 3) >> 2;
1710957b409SSimon J. Gerraty t += xm1d2_len_u32;
1720957b409SSimon J. Gerraty tlen -= xm1d2_len_u32;
1730957b409SSimon J. Gerraty
1740957b409SSimon J. Gerraty xlen = (x[0] + 31) >> 5;
1750957b409SSimon J. Gerraty asize = x[0] - 1 - EQ0(x[0] & 31);
1760957b409SSimon J. Gerraty x0i = br_i31_ninv31(x[1]);
1770957b409SSimon J. Gerraty while (n -- > 0) {
1780957b409SSimon J. Gerraty uint32_t *a, *t2;
1790957b409SSimon J. Gerraty uint32_t eq1, eqm1;
1800957b409SSimon J. Gerraty size_t t2len;
1810957b409SSimon J. Gerraty
1820957b409SSimon J. Gerraty /*
1830957b409SSimon J. Gerraty * Generate a random base. We don't need the base to be
1840957b409SSimon J. Gerraty * really uniform modulo x, so we just get a random
1850957b409SSimon J. Gerraty * number which is one bit shorter than x.
1860957b409SSimon J. Gerraty */
1870957b409SSimon J. Gerraty a = t;
1880957b409SSimon J. Gerraty a[0] = x[0];
1890957b409SSimon J. Gerraty a[xlen] = 0;
1900957b409SSimon J. Gerraty mkrand(rng, a, asize);
1910957b409SSimon J. Gerraty
1920957b409SSimon J. Gerraty /*
1930957b409SSimon J. Gerraty * Compute a^((x-1)/2) mod x. We assume here that the
1940957b409SSimon J. Gerraty * function will not fail (the temporary array is large
1950957b409SSimon J. Gerraty * enough).
1960957b409SSimon J. Gerraty */
1970957b409SSimon J. Gerraty t2 = t + 1 + xlen;
1980957b409SSimon J. Gerraty t2len = tlen - 1 - xlen;
1990957b409SSimon J. Gerraty if ((t2len & 1) != 0) {
2000957b409SSimon J. Gerraty /*
2010957b409SSimon J. Gerraty * Since the source array is 64-bit aligned and
2020957b409SSimon J. Gerraty * has an even number of elements (TEMPS), we
2030957b409SSimon J. Gerraty * can use the parity of the remaining length to
2040957b409SSimon J. Gerraty * detect and adjust alignment.
2050957b409SSimon J. Gerraty */
2060957b409SSimon J. Gerraty t2 ++;
2070957b409SSimon J. Gerraty t2len --;
2080957b409SSimon J. Gerraty }
2090957b409SSimon J. Gerraty mp31(a, xm1d2, xm1d2_len, x, x0i, t2, t2len);
2100957b409SSimon J. Gerraty
2110957b409SSimon J. Gerraty /*
2120957b409SSimon J. Gerraty * We must obtain either 1 or x-1. Note that x is odd,
2130957b409SSimon J. Gerraty * hence x-1 differs from x only in its low word (no
2140957b409SSimon J. Gerraty * carry).
2150957b409SSimon J. Gerraty */
2160957b409SSimon J. Gerraty eq1 = a[1] ^ 1;
2170957b409SSimon J. Gerraty eqm1 = a[1] ^ (x[1] - 1);
2180957b409SSimon J. Gerraty for (u = 2; u <= xlen; u ++) {
2190957b409SSimon J. Gerraty eq1 |= a[u];
2200957b409SSimon J. Gerraty eqm1 |= a[u] ^ x[u];
2210957b409SSimon J. Gerraty }
2220957b409SSimon J. Gerraty
2230957b409SSimon J. Gerraty if ((EQ0(eq1) | EQ0(eqm1)) == 0) {
2240957b409SSimon J. Gerraty return 0;
2250957b409SSimon J. Gerraty }
2260957b409SSimon J. Gerraty }
2270957b409SSimon J. Gerraty return 1;
2280957b409SSimon J. Gerraty }
2290957b409SSimon J. Gerraty
2300957b409SSimon J. Gerraty /*
2310957b409SSimon J. Gerraty * Create a random prime of the provided size. 'size' is the _encoded_
2320957b409SSimon J. Gerraty * bit length. The two top bits and the two bottom bits are set to 1.
2330957b409SSimon J. Gerraty */
2340957b409SSimon J. Gerraty static void
mkprime(const br_prng_class ** rng,uint32_t * x,uint32_t esize,uint32_t pubexp,uint32_t * t,size_t tlen,br_i31_modpow_opt_type mp31)2350957b409SSimon J. Gerraty mkprime(const br_prng_class **rng, uint32_t *x, uint32_t esize,
2360957b409SSimon J. Gerraty uint32_t pubexp, uint32_t *t, size_t tlen, br_i31_modpow_opt_type mp31)
2370957b409SSimon J. Gerraty {
2380957b409SSimon J. Gerraty size_t len;
2390957b409SSimon J. Gerraty
2400957b409SSimon J. Gerraty x[0] = esize;
2410957b409SSimon J. Gerraty len = (esize + 31) >> 5;
2420957b409SSimon J. Gerraty for (;;) {
2430957b409SSimon J. Gerraty size_t u;
2440957b409SSimon J. Gerraty uint32_t m3, m5, m7, m11;
2450957b409SSimon J. Gerraty int rounds, s7, s11;
2460957b409SSimon J. Gerraty
2470957b409SSimon J. Gerraty /*
2480957b409SSimon J. Gerraty * Generate random bits. We force the two top bits and the
2490957b409SSimon J. Gerraty * two bottom bits to 1.
2500957b409SSimon J. Gerraty */
2510957b409SSimon J. Gerraty mkrand(rng, x, esize);
2520957b409SSimon J. Gerraty if ((esize & 31) == 0) {
2530957b409SSimon J. Gerraty x[len] |= 0x60000000;
2540957b409SSimon J. Gerraty } else if ((esize & 31) == 1) {
2550957b409SSimon J. Gerraty x[len] |= 0x00000001;
2560957b409SSimon J. Gerraty x[len - 1] |= 0x40000000;
2570957b409SSimon J. Gerraty } else {
2580957b409SSimon J. Gerraty x[len] |= 0x00000003 << ((esize & 31) - 2);
2590957b409SSimon J. Gerraty }
2600957b409SSimon J. Gerraty x[1] |= 0x00000003;
2610957b409SSimon J. Gerraty
2620957b409SSimon J. Gerraty /*
2630957b409SSimon J. Gerraty * Trial division with low primes (3, 5, 7 and 11). We
2640957b409SSimon J. Gerraty * use the following properties:
2650957b409SSimon J. Gerraty *
2660957b409SSimon J. Gerraty * 2^2 = 1 mod 3
2670957b409SSimon J. Gerraty * 2^4 = 1 mod 5
2680957b409SSimon J. Gerraty * 2^3 = 1 mod 7
2690957b409SSimon J. Gerraty * 2^10 = 1 mod 11
2700957b409SSimon J. Gerraty */
2710957b409SSimon J. Gerraty m3 = 0;
2720957b409SSimon J. Gerraty m5 = 0;
2730957b409SSimon J. Gerraty m7 = 0;
2740957b409SSimon J. Gerraty m11 = 0;
2750957b409SSimon J. Gerraty s7 = 0;
2760957b409SSimon J. Gerraty s11 = 0;
2770957b409SSimon J. Gerraty for (u = 0; u < len; u ++) {
2780957b409SSimon J. Gerraty uint32_t w, w3, w5, w7, w11;
2790957b409SSimon J. Gerraty
2800957b409SSimon J. Gerraty w = x[1 + u];
2810957b409SSimon J. Gerraty w3 = (w & 0xFFFF) + (w >> 16); /* max: 98302 */
2820957b409SSimon J. Gerraty w5 = (w & 0xFFFF) + (w >> 16); /* max: 98302 */
2830957b409SSimon J. Gerraty w7 = (w & 0x7FFF) + (w >> 15); /* max: 98302 */
2840957b409SSimon J. Gerraty w11 = (w & 0xFFFFF) + (w >> 20); /* max: 1050622 */
2850957b409SSimon J. Gerraty
2860957b409SSimon J. Gerraty m3 += w3 << (u & 1);
2870957b409SSimon J. Gerraty m3 = (m3 & 0xFF) + (m3 >> 8); /* max: 1025 */
2880957b409SSimon J. Gerraty
2890957b409SSimon J. Gerraty m5 += w5 << ((4 - u) & 3);
2900957b409SSimon J. Gerraty m5 = (m5 & 0xFFF) + (m5 >> 12); /* max: 4479 */
2910957b409SSimon J. Gerraty
2920957b409SSimon J. Gerraty m7 += w7 << s7;
2930957b409SSimon J. Gerraty m7 = (m7 & 0x1FF) + (m7 >> 9); /* max: 1280 */
2940957b409SSimon J. Gerraty if (++ s7 == 3) {
2950957b409SSimon J. Gerraty s7 = 0;
2960957b409SSimon J. Gerraty }
2970957b409SSimon J. Gerraty
2980957b409SSimon J. Gerraty m11 += w11 << s11;
2990957b409SSimon J. Gerraty if (++ s11 == 10) {
3000957b409SSimon J. Gerraty s11 = 0;
3010957b409SSimon J. Gerraty }
3020957b409SSimon J. Gerraty m11 = (m11 & 0x3FF) + (m11 >> 10); /* max: 526847 */
3030957b409SSimon J. Gerraty }
3040957b409SSimon J. Gerraty
3050957b409SSimon J. Gerraty m3 = (m3 & 0x3F) + (m3 >> 6); /* max: 78 */
3060957b409SSimon J. Gerraty m3 = (m3 & 0x0F) + (m3 >> 4); /* max: 18 */
3070957b409SSimon J. Gerraty m3 = ((m3 * 43) >> 5) & 3;
3080957b409SSimon J. Gerraty
3090957b409SSimon J. Gerraty m5 = (m5 & 0xFF) + (m5 >> 8); /* max: 271 */
3100957b409SSimon J. Gerraty m5 = (m5 & 0x0F) + (m5 >> 4); /* max: 31 */
3110957b409SSimon J. Gerraty m5 -= 20 & -GT(m5, 19);
3120957b409SSimon J. Gerraty m5 -= 10 & -GT(m5, 9);
3130957b409SSimon J. Gerraty m5 -= 5 & -GT(m5, 4);
3140957b409SSimon J. Gerraty
3150957b409SSimon J. Gerraty m7 = (m7 & 0x3F) + (m7 >> 6); /* max: 82 */
3160957b409SSimon J. Gerraty m7 = (m7 & 0x07) + (m7 >> 3); /* max: 16 */
3170957b409SSimon J. Gerraty m7 = ((m7 * 147) >> 7) & 7;
3180957b409SSimon J. Gerraty
3190957b409SSimon J. Gerraty /*
3200957b409SSimon J. Gerraty * 2^5 = 32 = -1 mod 11.
3210957b409SSimon J. Gerraty */
3220957b409SSimon J. Gerraty m11 = (m11 & 0x3FF) + (m11 >> 10); /* max: 1536 */
3230957b409SSimon J. Gerraty m11 = (m11 & 0x3FF) + (m11 >> 10); /* max: 1023 */
3240957b409SSimon J. Gerraty m11 = (m11 & 0x1F) + 33 - (m11 >> 5); /* max: 64 */
3250957b409SSimon J. Gerraty m11 -= 44 & -GT(m11, 43);
3260957b409SSimon J. Gerraty m11 -= 22 & -GT(m11, 21);
3270957b409SSimon J. Gerraty m11 -= 11 & -GT(m11, 10);
3280957b409SSimon J. Gerraty
3290957b409SSimon J. Gerraty /*
3300957b409SSimon J. Gerraty * If any of these modulo is 0, then the candidate is
3310957b409SSimon J. Gerraty * not prime. Also, if pubexp is 3, 5, 7 or 11, and the
3320957b409SSimon J. Gerraty * corresponding modulus is 1, then the candidate must
3330957b409SSimon J. Gerraty * be rejected, because we need e to be invertible
3340957b409SSimon J. Gerraty * modulo p-1. We can use simple comparisons here
3350957b409SSimon J. Gerraty * because they won't leak information on a candidate
3360957b409SSimon J. Gerraty * that we keep, only on one that we reject (and is thus
3370957b409SSimon J. Gerraty * not secret).
3380957b409SSimon J. Gerraty */
3390957b409SSimon J. Gerraty if (m3 == 0 || m5 == 0 || m7 == 0 || m11 == 0) {
3400957b409SSimon J. Gerraty continue;
3410957b409SSimon J. Gerraty }
3420957b409SSimon J. Gerraty if ((pubexp == 3 && m3 == 1)
343*cc9e6590SSimon J. Gerraty || (pubexp == 5 && m5 == 1)
344*cc9e6590SSimon J. Gerraty || (pubexp == 7 && m7 == 1)
345*cc9e6590SSimon J. Gerraty || (pubexp == 11 && m11 == 1))
3460957b409SSimon J. Gerraty {
3470957b409SSimon J. Gerraty continue;
3480957b409SSimon J. Gerraty }
3490957b409SSimon J. Gerraty
3500957b409SSimon J. Gerraty /*
3510957b409SSimon J. Gerraty * More trial divisions.
3520957b409SSimon J. Gerraty */
3530957b409SSimon J. Gerraty if (!trial_divisions(x, t)) {
3540957b409SSimon J. Gerraty continue;
3550957b409SSimon J. Gerraty }
3560957b409SSimon J. Gerraty
3570957b409SSimon J. Gerraty /*
3580957b409SSimon J. Gerraty * Miller-Rabin algorithm. Since we selected a random
3590957b409SSimon J. Gerraty * integer, not a maliciously crafted integer, we can use
3600957b409SSimon J. Gerraty * relatively few rounds to lower the risk of a false
3610957b409SSimon J. Gerraty * positive (i.e. declaring prime a non-prime) under
3620957b409SSimon J. Gerraty * 2^(-80). It is not useful to lower the probability much
3630957b409SSimon J. Gerraty * below that, since that would be substantially below
3640957b409SSimon J. Gerraty * the probability of the hardware misbehaving. Sufficient
3650957b409SSimon J. Gerraty * numbers of rounds are extracted from the Handbook of
3660957b409SSimon J. Gerraty * Applied Cryptography, note 4.49 (page 149).
3670957b409SSimon J. Gerraty *
3680957b409SSimon J. Gerraty * Since we work on the encoded size (esize), we need to
3690957b409SSimon J. Gerraty * compare with encoded thresholds.
3700957b409SSimon J. Gerraty */
3710957b409SSimon J. Gerraty if (esize < 309) {
3720957b409SSimon J. Gerraty rounds = 12;
3730957b409SSimon J. Gerraty } else if (esize < 464) {
3740957b409SSimon J. Gerraty rounds = 9;
3750957b409SSimon J. Gerraty } else if (esize < 670) {
3760957b409SSimon J. Gerraty rounds = 6;
3770957b409SSimon J. Gerraty } else if (esize < 877) {
3780957b409SSimon J. Gerraty rounds = 4;
3790957b409SSimon J. Gerraty } else if (esize < 1341) {
3800957b409SSimon J. Gerraty rounds = 3;
3810957b409SSimon J. Gerraty } else {
3820957b409SSimon J. Gerraty rounds = 2;
3830957b409SSimon J. Gerraty }
3840957b409SSimon J. Gerraty
3850957b409SSimon J. Gerraty if (miller_rabin(rng, x, rounds, t, tlen, mp31)) {
3860957b409SSimon J. Gerraty return;
3870957b409SSimon J. Gerraty }
3880957b409SSimon J. Gerraty }
3890957b409SSimon J. Gerraty }
3900957b409SSimon J. Gerraty
3910957b409SSimon J. Gerraty /*
3920957b409SSimon J. Gerraty * Let p be a prime (p > 2^33, p = 3 mod 4). Let m = (p-1)/2, provided
3930957b409SSimon J. Gerraty * as parameter (with announced bit length equal to that of p). This
3940957b409SSimon J. Gerraty * function computes d = 1/e mod p-1 (for an odd integer e). Returned
3950957b409SSimon J. Gerraty * value is 1 on success, 0 on error (an error is reported if e is not
3960957b409SSimon J. Gerraty * invertible modulo p-1).
3970957b409SSimon J. Gerraty *
3980957b409SSimon J. Gerraty * The temporary buffer (t) must have room for at least 4 integers of
3990957b409SSimon J. Gerraty * the size of p.
4000957b409SSimon J. Gerraty */
4010957b409SSimon J. Gerraty static uint32_t
invert_pubexp(uint32_t * d,const uint32_t * m,uint32_t e,uint32_t * t)4020957b409SSimon J. Gerraty invert_pubexp(uint32_t *d, const uint32_t *m, uint32_t e, uint32_t *t)
4030957b409SSimon J. Gerraty {
4040957b409SSimon J. Gerraty uint32_t *f;
4050957b409SSimon J. Gerraty uint32_t r;
4060957b409SSimon J. Gerraty
4070957b409SSimon J. Gerraty f = t;
4080957b409SSimon J. Gerraty t += 1 + ((m[0] + 31) >> 5);
4090957b409SSimon J. Gerraty
4100957b409SSimon J. Gerraty /*
4110957b409SSimon J. Gerraty * Compute d = 1/e mod m. Since p = 3 mod 4, m is odd.
4120957b409SSimon J. Gerraty */
4130957b409SSimon J. Gerraty br_i31_zero(d, m[0]);
4140957b409SSimon J. Gerraty d[1] = 1;
4150957b409SSimon J. Gerraty br_i31_zero(f, m[0]);
4160957b409SSimon J. Gerraty f[1] = e & 0x7FFFFFFF;
4170957b409SSimon J. Gerraty f[2] = e >> 31;
4180957b409SSimon J. Gerraty r = br_i31_moddiv(d, f, m, br_i31_ninv31(m[1]), t);
4190957b409SSimon J. Gerraty
4200957b409SSimon J. Gerraty /*
4210957b409SSimon J. Gerraty * We really want d = 1/e mod p-1, with p = 2m. By the CRT,
4220957b409SSimon J. Gerraty * the result is either the d we got, or d + m.
4230957b409SSimon J. Gerraty *
4240957b409SSimon J. Gerraty * Let's write e*d = 1 + k*m, for some integer k. Integers e
4250957b409SSimon J. Gerraty * and m are odd. If d is odd, then e*d is odd, which implies
4260957b409SSimon J. Gerraty * that k must be even; in that case, e*d = 1 + (k/2)*2m, and
4270957b409SSimon J. Gerraty * thus d is already fine. Conversely, if d is even, then k
4280957b409SSimon J. Gerraty * is odd, and we must add m to d in order to get the correct
4290957b409SSimon J. Gerraty * result.
4300957b409SSimon J. Gerraty */
4310957b409SSimon J. Gerraty br_i31_add(d, m, (uint32_t)(1 - (d[1] & 1)));
4320957b409SSimon J. Gerraty
4330957b409SSimon J. Gerraty return r;
4340957b409SSimon J. Gerraty }
4350957b409SSimon J. Gerraty
4360957b409SSimon J. Gerraty /*
4370957b409SSimon J. Gerraty * Swap two buffers in RAM. They must be disjoint.
4380957b409SSimon J. Gerraty */
4390957b409SSimon J. Gerraty static void
bufswap(void * b1,void * b2,size_t len)4400957b409SSimon J. Gerraty bufswap(void *b1, void *b2, size_t len)
4410957b409SSimon J. Gerraty {
4420957b409SSimon J. Gerraty size_t u;
4430957b409SSimon J. Gerraty unsigned char *buf1, *buf2;
4440957b409SSimon J. Gerraty
4450957b409SSimon J. Gerraty buf1 = b1;
4460957b409SSimon J. Gerraty buf2 = b2;
4470957b409SSimon J. Gerraty for (u = 0; u < len; u ++) {
4480957b409SSimon J. Gerraty unsigned w;
4490957b409SSimon J. Gerraty
4500957b409SSimon J. Gerraty w = buf1[u];
4510957b409SSimon J. Gerraty buf1[u] = buf2[u];
4520957b409SSimon J. Gerraty buf2[u] = w;
4530957b409SSimon J. Gerraty }
4540957b409SSimon J. Gerraty }
4550957b409SSimon J. Gerraty
4560957b409SSimon J. Gerraty /* see inner.h */
4570957b409SSimon J. Gerraty uint32_t
br_rsa_i31_keygen_inner(const br_prng_class ** rng,br_rsa_private_key * sk,void * kbuf_priv,br_rsa_public_key * pk,void * kbuf_pub,unsigned size,uint32_t pubexp,br_i31_modpow_opt_type mp31)4580957b409SSimon J. Gerraty br_rsa_i31_keygen_inner(const br_prng_class **rng,
4590957b409SSimon J. Gerraty br_rsa_private_key *sk, void *kbuf_priv,
4600957b409SSimon J. Gerraty br_rsa_public_key *pk, void *kbuf_pub,
4610957b409SSimon J. Gerraty unsigned size, uint32_t pubexp, br_i31_modpow_opt_type mp31)
4620957b409SSimon J. Gerraty {
4630957b409SSimon J. Gerraty uint32_t esize_p, esize_q;
4640957b409SSimon J. Gerraty size_t plen, qlen, tlen;
4650957b409SSimon J. Gerraty uint32_t *p, *q, *t;
4660957b409SSimon J. Gerraty union {
4670957b409SSimon J. Gerraty uint32_t t32[TEMPS];
4680957b409SSimon J. Gerraty uint64_t t64[TEMPS >> 1]; /* for 64-bit alignment */
4690957b409SSimon J. Gerraty } tmp;
4700957b409SSimon J. Gerraty uint32_t r;
4710957b409SSimon J. Gerraty
4720957b409SSimon J. Gerraty if (size < BR_MIN_RSA_SIZE || size > BR_MAX_RSA_SIZE) {
4730957b409SSimon J. Gerraty return 0;
4740957b409SSimon J. Gerraty }
4750957b409SSimon J. Gerraty if (pubexp == 0) {
4760957b409SSimon J. Gerraty pubexp = 3;
4770957b409SSimon J. Gerraty } else if (pubexp == 1 || (pubexp & 1) == 0) {
4780957b409SSimon J. Gerraty return 0;
4790957b409SSimon J. Gerraty }
4800957b409SSimon J. Gerraty
4810957b409SSimon J. Gerraty esize_p = (size + 1) >> 1;
4820957b409SSimon J. Gerraty esize_q = size - esize_p;
4830957b409SSimon J. Gerraty sk->n_bitlen = size;
4840957b409SSimon J. Gerraty sk->p = kbuf_priv;
4850957b409SSimon J. Gerraty sk->plen = (esize_p + 7) >> 3;
4860957b409SSimon J. Gerraty sk->q = sk->p + sk->plen;
4870957b409SSimon J. Gerraty sk->qlen = (esize_q + 7) >> 3;
4880957b409SSimon J. Gerraty sk->dp = sk->q + sk->qlen;
4890957b409SSimon J. Gerraty sk->dplen = sk->plen;
4900957b409SSimon J. Gerraty sk->dq = sk->dp + sk->dplen;
4910957b409SSimon J. Gerraty sk->dqlen = sk->qlen;
4920957b409SSimon J. Gerraty sk->iq = sk->dq + sk->dqlen;
4930957b409SSimon J. Gerraty sk->iqlen = sk->plen;
4940957b409SSimon J. Gerraty
4950957b409SSimon J. Gerraty if (pk != NULL) {
4960957b409SSimon J. Gerraty pk->n = kbuf_pub;
4970957b409SSimon J. Gerraty pk->nlen = (size + 7) >> 3;
4980957b409SSimon J. Gerraty pk->e = pk->n + pk->nlen;
4990957b409SSimon J. Gerraty pk->elen = 4;
5000957b409SSimon J. Gerraty br_enc32be(pk->e, pubexp);
5010957b409SSimon J. Gerraty while (*pk->e == 0) {
5020957b409SSimon J. Gerraty pk->e ++;
5030957b409SSimon J. Gerraty pk->elen --;
5040957b409SSimon J. Gerraty }
5050957b409SSimon J. Gerraty }
5060957b409SSimon J. Gerraty
5070957b409SSimon J. Gerraty /*
5080957b409SSimon J. Gerraty * We now switch to encoded sizes.
5090957b409SSimon J. Gerraty *
5100957b409SSimon J. Gerraty * floor((x * 16913) / (2^19)) is equal to floor(x/31) for all
5110957b409SSimon J. Gerraty * integers x from 0 to 34966; the intermediate product fits on
5120957b409SSimon J. Gerraty * 30 bits, thus we can use MUL31().
5130957b409SSimon J. Gerraty */
5140957b409SSimon J. Gerraty esize_p += MUL31(esize_p, 16913) >> 19;
5150957b409SSimon J. Gerraty esize_q += MUL31(esize_q, 16913) >> 19;
5160957b409SSimon J. Gerraty plen = (esize_p + 31) >> 5;
5170957b409SSimon J. Gerraty qlen = (esize_q + 31) >> 5;
5180957b409SSimon J. Gerraty p = tmp.t32;
5190957b409SSimon J. Gerraty q = p + 1 + plen;
5200957b409SSimon J. Gerraty t = q + 1 + qlen;
5210957b409SSimon J. Gerraty tlen = ((sizeof tmp.t32) / sizeof(uint32_t)) - (2 + plen + qlen);
5220957b409SSimon J. Gerraty
5230957b409SSimon J. Gerraty /*
5240957b409SSimon J. Gerraty * When looking for primes p and q, we temporarily divide
5250957b409SSimon J. Gerraty * candidates by 2, in order to compute the inverse of the
5260957b409SSimon J. Gerraty * public exponent.
5270957b409SSimon J. Gerraty */
5280957b409SSimon J. Gerraty
5290957b409SSimon J. Gerraty for (;;) {
5300957b409SSimon J. Gerraty mkprime(rng, p, esize_p, pubexp, t, tlen, mp31);
5310957b409SSimon J. Gerraty br_i31_rshift(p, 1);
5320957b409SSimon J. Gerraty if (invert_pubexp(t, p, pubexp, t + 1 + plen)) {
5330957b409SSimon J. Gerraty br_i31_add(p, p, 1);
5340957b409SSimon J. Gerraty p[1] |= 1;
5350957b409SSimon J. Gerraty br_i31_encode(sk->p, sk->plen, p);
5360957b409SSimon J. Gerraty br_i31_encode(sk->dp, sk->dplen, t);
5370957b409SSimon J. Gerraty break;
5380957b409SSimon J. Gerraty }
5390957b409SSimon J. Gerraty }
5400957b409SSimon J. Gerraty
5410957b409SSimon J. Gerraty for (;;) {
5420957b409SSimon J. Gerraty mkprime(rng, q, esize_q, pubexp, t, tlen, mp31);
5430957b409SSimon J. Gerraty br_i31_rshift(q, 1);
5440957b409SSimon J. Gerraty if (invert_pubexp(t, q, pubexp, t + 1 + qlen)) {
5450957b409SSimon J. Gerraty br_i31_add(q, q, 1);
5460957b409SSimon J. Gerraty q[1] |= 1;
5470957b409SSimon J. Gerraty br_i31_encode(sk->q, sk->qlen, q);
5480957b409SSimon J. Gerraty br_i31_encode(sk->dq, sk->dqlen, t);
5490957b409SSimon J. Gerraty break;
5500957b409SSimon J. Gerraty }
5510957b409SSimon J. Gerraty }
5520957b409SSimon J. Gerraty
5530957b409SSimon J. Gerraty /*
5540957b409SSimon J. Gerraty * If p and q have the same size, then it is possible that q > p
5550957b409SSimon J. Gerraty * (when the target modulus size is odd, we generate p with a
5560957b409SSimon J. Gerraty * greater bit length than q). If q > p, we want to swap p and q
5570957b409SSimon J. Gerraty * (and also dp and dq) for two reasons:
5580957b409SSimon J. Gerraty * - The final step below (inversion of q modulo p) is easier if
5590957b409SSimon J. Gerraty * p > q.
5600957b409SSimon J. Gerraty * - While BearSSL's RSA code is perfectly happy with RSA keys such
5610957b409SSimon J. Gerraty * that p < q, some other implementations have restrictions and
5620957b409SSimon J. Gerraty * require p > q.
5630957b409SSimon J. Gerraty *
5640957b409SSimon J. Gerraty * Note that we can do a simple non-constant-time swap here,
5650957b409SSimon J. Gerraty * because the only information we leak here is that we insist on
5660957b409SSimon J. Gerraty * returning p and q such that p > q, which is not a secret.
5670957b409SSimon J. Gerraty */
5680957b409SSimon J. Gerraty if (esize_p == esize_q && br_i31_sub(p, q, 0) == 1) {
5690957b409SSimon J. Gerraty bufswap(p, q, (1 + plen) * sizeof *p);
5700957b409SSimon J. Gerraty bufswap(sk->p, sk->q, sk->plen);
5710957b409SSimon J. Gerraty bufswap(sk->dp, sk->dq, sk->dplen);
5720957b409SSimon J. Gerraty }
5730957b409SSimon J. Gerraty
5740957b409SSimon J. Gerraty /*
5750957b409SSimon J. Gerraty * We have produced p, q, dp and dq. We can now compute iq = 1/d mod p.
5760957b409SSimon J. Gerraty *
5770957b409SSimon J. Gerraty * We ensured that p >= q, so this is just a matter of updating the
5780957b409SSimon J. Gerraty * header word for q (and possibly adding an extra word).
5790957b409SSimon J. Gerraty *
5800957b409SSimon J. Gerraty * Theoretically, the call below may fail, in case we were
5810957b409SSimon J. Gerraty * extraordinarily unlucky, and p = q. Another failure case is if
5820957b409SSimon J. Gerraty * Miller-Rabin failed us _twice_, and p and q are non-prime and
5830957b409SSimon J. Gerraty * have a factor is common. We report the error mostly because it
5840957b409SSimon J. Gerraty * is cheap and we can, but in practice this never happens (or, at
5850957b409SSimon J. Gerraty * least, it happens way less often than hardware glitches).
5860957b409SSimon J. Gerraty */
5870957b409SSimon J. Gerraty q[0] = p[0];
5880957b409SSimon J. Gerraty if (plen > qlen) {
5890957b409SSimon J. Gerraty q[plen] = 0;
5900957b409SSimon J. Gerraty t ++;
5910957b409SSimon J. Gerraty tlen --;
5920957b409SSimon J. Gerraty }
5930957b409SSimon J. Gerraty br_i31_zero(t, p[0]);
5940957b409SSimon J. Gerraty t[1] = 1;
5950957b409SSimon J. Gerraty r = br_i31_moddiv(t, q, p, br_i31_ninv31(p[1]), t + 1 + plen);
5960957b409SSimon J. Gerraty br_i31_encode(sk->iq, sk->iqlen, t);
5970957b409SSimon J. Gerraty
5980957b409SSimon J. Gerraty /*
5990957b409SSimon J. Gerraty * Compute the public modulus too, if required.
6000957b409SSimon J. Gerraty */
6010957b409SSimon J. Gerraty if (pk != NULL) {
6020957b409SSimon J. Gerraty br_i31_zero(t, p[0]);
6030957b409SSimon J. Gerraty br_i31_mulacc(t, p, q);
6040957b409SSimon J. Gerraty br_i31_encode(pk->n, pk->nlen, t);
6050957b409SSimon J. Gerraty }
6060957b409SSimon J. Gerraty
6070957b409SSimon J. Gerraty return r;
6080957b409SSimon J. Gerraty }
609