Added RSA key generation code (i15, i31, i62).

This commit is contained in:
Thomas Pornin 2018-07-31 23:00:26 +02:00
parent a0054ad821
commit 491a45337d
15 changed files with 2865 additions and 13 deletions

View file

@ -48,10 +48,22 @@
* already set their root keys to RSA-4096, so we should be able to
* process such keys.
*
* This value MUST be a multiple of 64.
* This value MUST be a multiple of 64. This value MUST NOT exceed 47666
* (some computations in RSA key generation rely on the factor size being
* no more than 23833 bits). RSA key sizes beyond 3072 bits don't make a
* lot of sense anyway.
*/
#define BR_MAX_RSA_SIZE 4096
/*
* Minimum size for a RSA modulus (in bits); this value is used only to
* filter out invalid parameters for key pair generation. Normally,
* applications should not use RSA keys smaller than 2048 bits; but some
* specific cases might need shorter keys, for legacy or research
* purposes.
*/
#define BR_MIN_RSA_SIZE 512
/*
* Maximum size for a RSA factor (in bits). This is for RSA private-key
* operations. Default is to support factors up to a bit more than half
@ -1474,6 +1486,23 @@ uint32_t br_i31_modpow_opt(uint32_t *x, const unsigned char *e, size_t elen,
*/
void br_i31_mulacc(uint32_t *d, const uint32_t *a, const uint32_t *b);
/*
* Compute x/y mod m, result in x. Values x and y must be between 0 and
* m-1, and have the same announced bit length as m. Modulus m must be
* odd. The "m0i" parameter is equal to -1/m mod 2^31. The array 't'
* must point to a temporary area that can hold at least three integers
* of the size of m.
*
* m may not overlap x and y. x and y may overlap each other (this can
* be useful to test whether a value is invertible modulo m). t must be
* disjoint from all other arrays.
*
* Returned value is 1 on success, 0 otherwise. Success is attained if
* y is invertible modulo m.
*/
uint32_t br_i31_moddiv(uint32_t *x, const uint32_t *y,
const uint32_t *m, uint32_t m0i, uint32_t *t);
/* ==================================================================== */
/*
@ -1528,9 +1557,36 @@ void br_i15_reduce(uint16_t *x, const uint16_t *a, const uint16_t *m);
void br_i15_mulacc(uint16_t *d, const uint16_t *a, const uint16_t *b);
uint32_t br_i15_moddiv(uint16_t *x, const uint16_t *y,
const uint16_t *m, uint16_t m0i, uint16_t *t);
/*
* Variant of br_i31_modpow_opt() that internally uses 64x64->128
* multiplications. It expects the same parameters as br_i31_modpow_opt(),
* except that the temporaries should be 64-bit integers, not 32-bit
* integers.
*/
uint32_t br_i62_modpow_opt(uint32_t *x31, const unsigned char *e, size_t elen,
const uint32_t *m31, uint32_t m0i31, uint64_t *tmp, size_t twlen);
/*
* Type for a function with the same API as br_i31_modpow_opt() (some
* implementations of this type may have stricter alignment requirements
* on the temporaries).
*/
typedef uint32_t (*br_i31_modpow_opt_type)(uint32_t *x,
const unsigned char *e, size_t elen,
const uint32_t *m, uint32_t m0i, uint32_t *tmp, size_t twlen);
/*
* Wrapper for br_i62_modpow_opt() that uses the same type as
* br_i31_modpow_opt(); however, it requires its 'tmp' argument to the
* 64-bit aligned.
*/
uint32_t br_i62_modpow_opt_as_i31(uint32_t *x,
const unsigned char *e, size_t elen,
const uint32_t *m, uint32_t m0i, uint32_t *tmp, size_t twlen);
/* ==================================================================== */
static inline size_t
@ -1912,6 +1968,15 @@ uint32_t br_rsa_oaep_unpad(const br_hash_class *dig,
void br_mgf1_xor(void *data, size_t len,
const br_hash_class *dig, const void *seed, size_t seed_len);
/*
* Inner function for RSA key generation; used by the "i31" and "i62"
* implementations.
*/
uint32_t br_rsa_i31_keygen_inner(const br_prng_class **rng,
br_rsa_private_key *sk, unsigned char *kbuf_priv,
br_rsa_public_key *pk, unsigned char *kbuf_pub,
unsigned size, uint32_t pubexp, br_i31_modpow_opt_type mp31);
/* ==================================================================== */
/*
* Elliptic curves.

465
src/int/i15_moddiv.c Normal file
View file

@ -0,0 +1,465 @@
/*
* Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include "inner.h"
/*
* In this file, we handle big integers with a custom format, i.e.
* without the usual one-word header. Value is split into 15-bit words,
* each stored in a 16-bit slot (top bit is zero) in little-endian
* order. The length (in words) is provided explicitly. In some cases,
* the value can be negative (using two's complement representation). In
* some cases, the top word is allowed to have a 16th bit.
*/
/*
* Negate big integer conditionally. The value consists of 'len' words,
* with 15 bits in each word (the top bit of each word should be 0,
* except possibly for the last word). If 'ctl' is 1, the negation is
* computed; otherwise, if 'ctl' is 0, then the value is unchanged.
*/
static void
cond_negate(uint16_t *a, size_t len, uint32_t ctl)
{
size_t k;
uint32_t cc, xm;
cc = ctl;
xm = 0x7FFF & -ctl;
for (k = 0; k < len; k ++) {
uint32_t aw;
aw = a[k];
aw = (aw ^ xm) + cc;
a[k] = aw & 0x7FFF;
cc = (aw >> 15) & 1;
}
}
/*
* Finish modular reduction. Rules on input parameters:
*
* if neg = 1, then -m <= a < 0
* if neg = 0, then 0 <= a < 2*m
*
* If neg = 0, then the top word of a[] may use 16 bits.
*
* Also, modulus m must be odd.
*/
static void
finish_mod(uint16_t *a, size_t len, const uint16_t *m, uint32_t neg)
{
size_t k;
uint32_t cc, xm, ym;
/*
* First pass: compare a (assumed nonnegative) with m.
*/
cc = 0;
for (k = 0; k < len; k ++) {
uint32_t aw, mw;
aw = a[k];
mw = m[k];
cc = (aw - mw - cc) >> 31;
}
/*
* At this point:
* if neg = 1, then we must add m (regardless of cc)
* if neg = 0 and cc = 0, then we must subtract m
* if neg = 0 and cc = 1, then we must do nothing
*/
xm = 0x7FFF & -neg;
ym = -(neg | (1 - cc));
cc = neg;
for (k = 0; k < len; k ++) {
uint32_t aw, mw;
aw = a[k];
mw = (m[k] ^ xm) & ym;
aw = aw - mw - cc;
a[k] = aw & 0x7FFF;
cc = aw >> 31;
}
}
/*
* Compute:
* a <- (a*pa+b*pb)/(2^15)
* b <- (a*qa+b*qb)/(2^15)
* The division is assumed to be exact (i.e. the low word is dropped).
* If the final a is negative, then it is negated. Similarly for b.
* Returned value is the combination of two bits:
* bit 0: 1 if a had to be negated, 0 otherwise
* bit 1: 1 if b had to be negated, 0 otherwise
*
* Factors pa, pb, qa and qb must be at most 2^15 in absolute value.
* Source integers a and b must be nonnegative; top word is not allowed
* to contain an extra 16th bit.
*/
static uint32_t
co_reduce(uint16_t *a, uint16_t *b, size_t len,
int32_t pa, int32_t pb, int32_t qa, int32_t qb)
{
size_t k;
int32_t cca, ccb;
uint32_t nega, negb;
cca = 0;
ccb = 0;
for (k = 0; k < len; k ++) {
uint32_t wa, wb, za, zb;
uint16_t tta, ttb;
/*
* Since:
* |pa| <= 2^15
* |pb| <= 2^15
* 0 <= wa <= 2^15 - 1
* 0 <= wb <= 2^15 - 1
* |cca| <= 2^16 - 1
* Then:
* |za| <= (2^15-1)*(2^16) + (2^16-1) = 2^31 - 1
*
* Thus, the new value of cca is such that |cca| <= 2^16 - 1.
* The same applies to ccb.
*/
wa = a[k];
wb = b[k];
za = wa * (uint32_t)pa + wb * (uint32_t)pb + (uint32_t)cca;
zb = wa * (uint32_t)qa + wb * (uint32_t)qb + (uint32_t)ccb;
if (k > 0) {
a[k - 1] = za & 0x7FFF;
b[k - 1] = zb & 0x7FFF;
}
tta = za >> 15;
ttb = zb >> 15;
cca = *(int16_t *)&tta;
ccb = *(int16_t *)&ttb;
}
a[len - 1] = (uint16_t)cca;
b[len - 1] = (uint16_t)ccb;
nega = (uint32_t)cca >> 31;
negb = (uint32_t)ccb >> 31;
cond_negate(a, len, nega);
cond_negate(b, len, negb);
return nega | (negb << 1);
}
/*
* Compute:
* a <- (a*pa+b*pb)/(2^15) mod m
* b <- (a*qa+b*qb)/(2^15) mod m
*
* m0i is equal to -1/m[0] mod 2^15.
*
* Factors pa, pb, qa and qb must be at most 2^15 in absolute value.
* Source integers a and b must be nonnegative; top word is not allowed
* to contain an extra 16th bit.
*/
static void
co_reduce_mod(uint16_t *a, uint16_t *b, size_t len,
int32_t pa, int32_t pb, int32_t qa, int32_t qb,
const uint16_t *m, uint16_t m0i)
{
size_t k;
int32_t cca, ccb, fa, fb;
cca = 0;
ccb = 0;
fa = ((a[0] * (uint32_t)pa + b[0] * (uint32_t)pb) * m0i) & 0x7FFF;
fb = ((a[0] * (uint32_t)qa + b[0] * (uint32_t)qb) * m0i) & 0x7FFF;
for (k = 0; k < len; k ++) {
uint32_t wa, wb, za, zb;
uint32_t tta, ttb;
/*
* In this loop, carries 'cca' and 'ccb' always fit on
* 17 bits (in absolute value).
*/
wa = a[k];
wb = b[k];
za = wa * (uint32_t)pa + wb * (uint32_t)pb
+ m[k] * (uint32_t)fa + (uint32_t)cca;
zb = wa * (uint32_t)qa + wb * (uint32_t)qb
+ m[k] * (uint32_t)fb + (uint32_t)ccb;
if (k > 0) {
a[k - 1] = za & 0x7FFF;
b[k - 1] = zb & 0x7FFF;
}
/*
* The XOR-and-sub construction below does an arithmetic
* right shift in a portable way (technically, right-shifting
* a negative signed value is implementation-defined in C).
*/
#define M ((uint32_t)1 << 16)
tta = za >> 15;
ttb = zb >> 15;
tta = (tta ^ M) - M;
ttb = (ttb ^ M) - M;
cca = *(int32_t *)&tta;
ccb = *(int32_t *)&ttb;
#undef M
}
a[len - 1] = (uint32_t)cca;
b[len - 1] = (uint32_t)ccb;
/*
* At this point:
* -m <= a < 2*m
* -m <= b < 2*m
* (this is a case of Montgomery reduction)
* The top word of 'a' and 'b' may have a 16-th bit set.
* We may have to add or subtract the modulus.
*/
finish_mod(a, len, m, (uint32_t)cca >> 31);
finish_mod(b, len, m, (uint32_t)ccb >> 31);
}
/* see inner.h */
uint32_t
br_i15_moddiv(uint16_t *x, const uint16_t *y, const uint16_t *m, uint16_t m0i,
uint16_t *t)
{
/*
* Algorithm is an extended binary GCD. We maintain four values
* a, b, u and v, with the following invariants:
*
* a * x = y * u mod m
* b * x = y * v mod m
*
* Starting values are:
*
* a = y
* b = m
* u = x
* v = 0
*
* The formal definition of the algorithm is a sequence of steps:
*
* - If a is even, then a <- a/2 and u <- u/2 mod m.
* - Otherwise, if b is even, then b <- b/2 and v <- v/2 mod m.
* - Otherwise, if a > b, then a <- (a-b)/2 and u <- (u-v)/2 mod m.
* - Otherwise, b <- (b-a)/2 and v <- (v-u)/2 mod m.
*
* Algorithm stops when a = b. At that point, they both are equal
* to GCD(y,m); the modular division succeeds if that value is 1.
* The result of the modular division is then u (or v: both are
* equal at that point).
*
* Each step makes either a or b shrink by at least one bit; hence,
* if m has bit length k bits, then 2k-2 steps are sufficient.
*
*
* Though complexity is quadratic in the size of m, the bit-by-bit
* processing is not very efficient. We can speed up processing by
* remarking that the decisions are taken based only on observation
* of the top and low bits of a and b.
*
* In the loop below, at each iteration, we use the two top words
* of a and b, and the low words of a and b, to compute reduction
* parameters pa, pb, qa and qb such that the new values for a
* and b are:
*
* a' = (a*pa + b*pb) / (2^15)
* b' = (a*qa + b*qb) / (2^15)
*
* the division being exact.
*
* Since the choices are based on the top words, they may be slightly
* off, requiring an optional correction: if a' < 0, then we replace
* pa with -pa, and pb with -pb. The total length of a and b is
* thus reduced by at least 14 bits at each iteration.
*
* The stopping conditions are still the same, though: when a
* and b become equal, they must be both odd (since m is odd,
* the GCD cannot be even), therefore the next operation is a
* subtraction, and one of the values becomes 0. At that point,
* nothing else happens, i.e. one value is stuck at 0, and the
* other one is the GCD.
*/
size_t len, k;
uint16_t *a, *b, *u, *v;
uint32_t num, r;
len = (m[0] + 15) >> 4;
a = t;
b = a + len;
u = x + 1;
v = b + len;
memcpy(a, y + 1, len * sizeof *y);
memcpy(b, m + 1, len * sizeof *m);
memset(v, 0, len * sizeof *v);
/*
* Loop below ensures that a and b are reduced by some bits each,
* for a total of at least 14 bits.
*/
for (num = ((m[0] - (m[0] >> 4)) << 1) + 14; num >= 14; num -= 14) {
size_t j;
uint32_t c0, c1;
uint32_t a0, a1, b0, b1;
uint32_t a_hi, b_hi, a_lo, b_lo;
int32_t pa, pb, qa, qb;
int i;
/*
* Extract top words of a and b. If j is the highest
* index >= 1 such that a[j] != 0 or b[j] != 0, then we want
* (a[j] << 15) + a[j - 1], and (b[j] << 15) + b[j - 1].
* If a and b are down to one word each, then we use a[0]
* and b[0].
*/
c0 = (uint32_t)-1;
c1 = (uint32_t)-1;
a0 = 0;
a1 = 0;
b0 = 0;
b1 = 0;
j = len;
while (j -- > 0) {
uint32_t aw, bw;
aw = a[j];
bw = b[j];
a0 ^= (a0 ^ aw) & c0;
a1 ^= (a1 ^ aw) & c1;
b0 ^= (b0 ^ bw) & c0;
b1 ^= (b1 ^ bw) & c1;
c1 = c0;
c0 &= (((aw | bw) + 0xFFFF) >> 16) - (uint32_t)1;
}
/*
* If c1 = 0, then we grabbed two words for a and b.
* If c1 != 0 but c0 = 0, then we grabbed one word. It
* is not possible that c1 != 0 and c0 != 0, because that
* would mean that both integers are zero.
*/
a1 |= a0 & c1;
a0 &= ~c1;
b1 |= b0 & c1;
b0 &= ~c1;
a_hi = (a0 << 15) + a1;
b_hi = (b0 << 15) + b1;
a_lo = a[0];
b_lo = b[0];
/*
* Compute reduction factors:
*
* a' = a*pa + b*pb
* b' = a*qa + b*qb
*
* such that a' and b' are both multiple of 2^15, but are
* only marginally larger than a and b.
*/
pa = 1;
pb = 0;
qa = 0;
qb = 1;
for (i = 0; i < 15; i ++) {
/*
* At each iteration:
*
* a <- (a-b)/2 if: a is odd, b is odd, a_hi > b_hi
* b <- (b-a)/2 if: a is odd, b is odd, a_hi <= b_hi
* a <- a/2 if: a is even
* b <- b/2 if: a is odd, b is even
*
* We multiply a_lo and b_lo by 2 at each
* iteration, thus a division by 2 really is a
* non-multiplication by 2.
*/
uint32_t r, oa, ob, cAB, cBA, cA;
/*
* cAB = 1 if b must be subtracted from a
* cBA = 1 if a must be subtracted from b
* cA = 1 if a is divided by 2, 0 otherwise
*
* Rules:
*
* cAB and cBA cannot be both 1.
* if a is not divided by 2, b is.
*/
r = GT(a_hi, b_hi);
oa = (a_lo >> i) & 1;
ob = (b_lo >> i) & 1;
cAB = oa & ob & r;
cBA = oa & ob & NOT(r);
cA = cAB | NOT(oa);
/*
* Conditional subtractions.
*/
a_lo -= b_lo & -cAB;
a_hi -= b_hi & -cAB;
pa -= qa & -(int32_t)cAB;
pb -= qb & -(int32_t)cAB;
b_lo -= a_lo & -cBA;
b_hi -= a_hi & -cBA;
qa -= pa & -(int32_t)cBA;
qb -= pb & -(int32_t)cBA;
/*
* Shifting.
*/
a_lo += a_lo & (cA - 1);
pa += pa & ((int32_t)cA - 1);
pb += pb & ((int32_t)cA - 1);
a_hi ^= (a_hi ^ (a_hi >> 1)) & -cA;
b_lo += b_lo & -cA;
qa += qa & -(int32_t)cA;
qb += qb & -(int32_t)cA;
b_hi ^= (b_hi ^ (b_hi >> 1)) & (cA - 1);
}
/*
* Replace a and b with new values a' and b'.
*/
r = co_reduce(a, b, len, pa, pb, qa, qb);
pa -= pa * ((r & 1) << 1);
pb -= pb * ((r & 1) << 1);
qa -= qa * (r & 2);
qb -= qb * (r & 2);
co_reduce_mod(u, v, len, pa, pb, qa, qb, m + 1, m0i);
}
/*
* Now one of the arrays should be 0, and the other contains
* the GCD. If a is 0, then u is 0 as well, and v contains
* the division result.
* Result is correct if and only if GCD is 1.
*/
r = (a[0] | b[0]) ^ 1;
u[0] |= v[0];
for (k = 1; k < len; k ++) {
r |= a[k] | b[k];
u[k] |= v[k];
}
return EQ0(r);
}

488
src/int/i31_moddiv.c Normal file
View file

@ -0,0 +1,488 @@
/*
* Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include "inner.h"
/*
* In this file, we handle big integers with a custom format, i.e.
* without the usual one-word header. Value is split into 31-bit words,
* each stored in a 32-bit slot (top bit is zero) in little-endian
* order. The length (in words) is provided explicitly. In some cases,
* the value can be negative (using two's complement representation). In
* some cases, the top word is allowed to have a 32th bit.
*/
/*
* Negate big integer conditionally. The value consists of 'len' words,
* with 31 bits in each word (the top bit of each word should be 0,
* except possibly for the last word). If 'ctl' is 1, the negation is
* computed; otherwise, if 'ctl' is 0, then the value is unchanged.
*/
static void
cond_negate(uint32_t *a, size_t len, uint32_t ctl)
{
size_t k;
uint32_t cc, xm;
cc = ctl;
xm = -ctl >> 1;
for (k = 0; k < len; k ++) {
uint32_t aw;
aw = a[k];
aw = (aw ^ xm) + cc;
a[k] = aw & 0x7FFFFFFF;
cc = aw >> 31;
}
}
/*
* Finish modular reduction. Rules on input parameters:
*
* if neg = 1, then -m <= a < 0
* if neg = 0, then 0 <= a < 2*m
*
* If neg = 0, then the top word of a[] may use 32 bits.
*
* Also, modulus m must be odd.
*/
static void
finish_mod(uint32_t *a, size_t len, const uint32_t *m, uint32_t neg)
{
size_t k;
uint32_t cc, xm, ym;
/*
* First pass: compare a (assumed nonnegative) with m.
* Note that if the final word uses the top extra bit, then
* subtracting m must yield a value less than 2^31, since we
* assumed that a < 2*m.
*/
cc = 0;
for (k = 0; k < len; k ++) {
uint32_t aw, mw;
aw = a[k];
mw = m[k];
cc = (aw - mw - cc) >> 31;
}
/*
* At this point:
* if neg = 1, then we must add m (regardless of cc)
* if neg = 0 and cc = 0, then we must subtract m
* if neg = 0 and cc = 1, then we must do nothing
*/
xm = -neg >> 1;
ym = -(neg | (1 - cc));
cc = neg;
for (k = 0; k < len; k ++) {
uint32_t aw, mw;
aw = a[k];
mw = (m[k] ^ xm) & ym;
aw = aw - mw - cc;
a[k] = aw & 0x7FFFFFFF;
cc = aw >> 31;
}
}
/*
* Compute:
* a <- (a*pa+b*pb)/(2^31)
* b <- (a*qa+b*qb)/(2^31)
* The division is assumed to be exact (i.e. the low word is dropped).
* If the final a is negative, then it is negated. Similarly for b.
* Returned value is the combination of two bits:
* bit 0: 1 if a had to be negated, 0 otherwise
* bit 1: 1 if b had to be negated, 0 otherwise
*
* Factors pa, pb, qa and qb must be at most 2^31 in absolute value.
* Source integers a and b must be nonnegative; top word is not allowed
* to contain an extra 32th bit.
*/
static uint32_t
co_reduce(uint32_t *a, uint32_t *b, size_t len,
int64_t pa, int64_t pb, int64_t qa, int64_t qb)
{
size_t k;
int64_t cca, ccb;
uint32_t nega, negb;
cca = 0;
ccb = 0;
for (k = 0; k < len; k ++) {
uint32_t wa, wb;
uint64_t za, zb;
uint64_t tta, ttb;
/*
* Since:
* |pa| <= 2^31
* |pb| <= 2^31
* 0 <= wa <= 2^31 - 1
* 0 <= wb <= 2^31 - 1
* |cca| <= 2^32 - 1
* Then:
* |za| <= (2^31-1)*(2^32) + (2^32-1) = 2^63 - 1
*
* Thus, the new value of cca is such that |cca| <= 2^32 - 1.
* The same applies to ccb.
*/
wa = a[k];
wb = b[k];
za = wa * (uint64_t)pa + wb * (uint64_t)pb + (uint64_t)cca;
zb = wa * (uint64_t)qa + wb * (uint64_t)qb + (uint64_t)ccb;
if (k > 0) {
a[k - 1] = za & 0x7FFFFFFF;
b[k - 1] = zb & 0x7FFFFFFF;
}
/*
* For the new values of cca and ccb, we need a signed
* right-shift; since, in C, right-shifting a signed
* negative value is implementation-defined, we use a
* custom portable sign extension expression.
*/
#define M ((uint64_t)1 << 32)
tta = za >> 31;
ttb = zb >> 31;
tta = (tta ^ M) - M;
ttb = (ttb ^ M) - M;
cca = *(int64_t *)&tta;
ccb = *(int64_t *)&ttb;
#undef M
}
a[len - 1] = (uint32_t)cca;
b[len - 1] = (uint32_t)ccb;
nega = (uint32_t)((uint64_t)cca >> 63);
negb = (uint32_t)((uint64_t)ccb >> 63);
cond_negate(a, len, nega);
cond_negate(b, len, negb);
return nega | (negb << 1);
}
/*
* Compute:
* a <- (a*pa+b*pb)/(2^31) mod m
* b <- (a*qa+b*qb)/(2^31) mod m
*
* m0i is equal to -1/m[0] mod 2^31.
*
* Factors pa, pb, qa and qb must be at most 2^31 in absolute value.
* Source integers a and b must be nonnegative; top word is not allowed
* to contain an extra 32th bit.
*/
static void
co_reduce_mod(uint32_t *a, uint32_t *b, size_t len,
int64_t pa, int64_t pb, int64_t qa, int64_t qb,
const uint32_t *m, uint32_t m0i)
{
size_t k;
int64_t cca, ccb;
uint32_t fa, fb;
cca = 0;
ccb = 0;
fa = ((a[0] * (uint32_t)pa + b[0] * (uint32_t)pb) * m0i) & 0x7FFFFFFF;
fb = ((a[0] * (uint32_t)qa + b[0] * (uint32_t)qb) * m0i) & 0x7FFFFFFF;
for (k = 0; k < len; k ++) {
uint32_t wa, wb;
uint64_t za, zb;
uint64_t tta, ttb;
/*
* In this loop, carries 'cca' and 'ccb' always fit on
* 33 bits (in absolute value).
*/
wa = a[k];
wb = b[k];
za = wa * (uint64_t)pa + wb * (uint64_t)pb
+ m[k] * (uint64_t)fa + (uint64_t)cca;
zb = wa * (uint64_t)qa + wb * (uint64_t)qb
+ m[k] * (uint64_t)fb + (uint64_t)ccb;
if (k > 0) {
a[k - 1] = (uint32_t)za & 0x7FFFFFFF;
b[k - 1] = (uint32_t)zb & 0x7FFFFFFF;
}
#define M ((uint64_t)1 << 32)
tta = za >> 31;
ttb = zb >> 31;
tta = (tta ^ M) - M;
ttb = (ttb ^ M) - M;
cca = *(int64_t *)&tta;
ccb = *(int64_t *)&ttb;
#undef M
}
a[len - 1] = (uint32_t)cca;
b[len - 1] = (uint32_t)ccb;
/*
* At this point:
* -m <= a < 2*m
* -m <= b < 2*m
* (this is a case of Montgomery reduction)
* The top word of 'a' and 'b' may have a 32-th bit set.
* We may have to add or subtract the modulus.
*/
finish_mod(a, len, m, (uint32_t)((uint64_t)cca >> 63));
finish_mod(b, len, m, (uint32_t)((uint64_t)ccb >> 63));
}
/* see inner.h */
uint32_t
br_i31_moddiv(uint32_t *x, const uint32_t *y, const uint32_t *m, uint32_t m0i,
uint32_t *t)
{
/*
* Algorithm is an extended binary GCD. We maintain four values
* a, b, u and v, with the following invariants:
*
* a * x = y * u mod m
* b * x = y * v mod m
*
* Starting values are:
*
* a = y
* b = m
* u = x
* v = 0
*
* The formal definition of the algorithm is a sequence of steps:
*
* - If a is even, then a <- a/2 and u <- u/2 mod m.
* - Otherwise, if b is even, then b <- b/2 and v <- v/2 mod m.
* - Otherwise, if a > b, then a <- (a-b)/2 and u <- (u-v)/2 mod m.
* - Otherwise, b <- (b-a)/2 and v <- (v-u)/2 mod m.
*
* Algorithm stops when a = b. At that point, they both are equal
* to GCD(y,m); the modular division succeeds if that value is 1.
* The result of the modular division is then u (or v: both are
* equal at that point).
*
* Each step makes either a or b shrink by at least one bit; hence,
* if m has bit length k bits, then 2k-2 steps are sufficient.
*
*
* Though complexity is quadratic in the size of m, the bit-by-bit
* processing is not very efficient. We can speed up processing by
* remarking that the decisions are taken based only on observation
* of the top and low bits of a and b.
*
* In the loop below, at each iteration, we use the two top words
* of a and b, and the low words of a and b, to compute reduction
* parameters pa, pb, qa and qb such that the new values for a
* and b are:
*
* a' = (a*pa + b*pb) / (2^31)
* b' = (a*qa + b*qb) / (2^31)
*
* the division being exact.
*
* Since the choices are based on the top words, they may be slightly
* off, requiring an optional correction: if a' < 0, then we replace
* pa with -pa, and pb with -pb. The total length of a and b is
* thus reduced by at least 30 bits at each iteration.
*
* The stopping conditions are still the same, though: when a
* and b become equal, they must be both odd (since m is odd,
* the GCD cannot be even), therefore the next operation is a
* subtraction, and one of the values becomes 0. At that point,
* nothing else happens, i.e. one value is stuck at 0, and the
* other one is the GCD.
*/
size_t len, k;
uint32_t *a, *b, *u, *v;
uint32_t num, r;
len = (m[0] + 31) >> 5;
a = t;
b = a + len;
u = x + 1;
v = b + len;
memcpy(a, y + 1, len * sizeof *y);
memcpy(b, m + 1, len * sizeof *m);
memset(v, 0, len * sizeof *v);
/*
* Loop below ensures that a and b are reduced by some bits each,
* for a total of at least 30 bits.
*/
for (num = ((m[0] - (m[0] >> 5)) << 1) + 30; num >= 30; num -= 30) {
size_t j;
uint32_t c0, c1;
uint32_t a0, a1, b0, b1;
uint64_t a_hi, b_hi;
uint32_t a_lo, b_lo;
int64_t pa, pb, qa, qb;
int i;
/*
* Extract top words of a and b. If j is the highest
* index >= 1 such that a[j] != 0 or b[j] != 0, then we want
* (a[j] << 31) + a[j - 1], and (b[j] << 31) + b[j - 1].
* If a and b are down to one word each, then we use a[0]
* and b[0].
*/
c0 = (uint32_t)-1;
c1 = (uint32_t)-1;
a0 = 0;
a1 = 0;
b0 = 0;
b1 = 0;
j = len;
while (j -- > 0) {
uint32_t aw, bw;
aw = a[j];
bw = b[j];
a0 ^= (a0 ^ aw) & c0;
a1 ^= (a1 ^ aw) & c1;
b0 ^= (b0 ^ bw) & c0;
b1 ^= (b1 ^ bw) & c1;
c1 = c0;
c0 &= (((aw | bw) + 0x7FFFFFFF) >> 31) - (uint32_t)1;
}
/*
* If c1 = 0, then we grabbed two words for a and b.
* If c1 != 0 but c0 = 0, then we grabbed one word. It
* is not possible that c1 != 0 and c0 != 0, because that
* would mean that both integers are zero.
*/
a1 |= a0 & c1;
a0 &= ~c1;
b1 |= b0 & c1;
b0 &= ~c1;
a_hi = ((uint64_t)a0 << 31) + a1;
b_hi = ((uint64_t)b0 << 31) + b1;
a_lo = a[0];
b_lo = b[0];
/*
* Compute reduction factors:
*
* a' = a*pa + b*pb
* b' = a*qa + b*qb
*
* such that a' and b' are both multiple of 2^31, but are
* only marginally larger than a and b.
*/
pa = 1;
pb = 0;
qa = 0;
qb = 1;
for (i = 0; i < 31; i ++) {
/*
* At each iteration:
*
* a <- (a-b)/2 if: a is odd, b is odd, a_hi > b_hi
* b <- (b-a)/2 if: a is odd, b is odd, a_hi <= b_hi
* a <- a/2 if: a is even
* b <- b/2 if: a is odd, b is even
*
* We multiply a_lo and b_lo by 2 at each
* iteration, thus a division by 2 really is a
* non-multiplication by 2.
*/
uint32_t r, oa, ob, cAB, cBA, cA;
uint64_t rz;
/*
* r = GT(a_hi, b_hi)
* But the GT() function works on uint32_t operands,
* so we inline a 64-bit version here.
*/
rz = b_hi - a_hi;
r = (uint32_t)((rz ^ ((a_hi ^ b_hi)
& (a_hi ^ rz))) >> 63);
/*
* cAB = 1 if b must be subtracted from a
* cBA = 1 if a must be subtracted from b
* cA = 1 if a is divided by 2, 0 otherwise
*
* Rules:
*
* cAB and cBA cannot be both 1.
* if a is not divided by 2, b is.
*/
oa = (a_lo >> i) & 1;
ob = (b_lo >> i) & 1;
cAB = oa & ob & r;
cBA = oa & ob & NOT(r);
cA = cAB | NOT(oa);
/*
* Conditional subtractions.
*/
a_lo -= b_lo & -cAB;
a_hi -= b_hi & -(uint64_t)cAB;
pa -= qa & -(int64_t)cAB;
pb -= qb & -(int64_t)cAB;
b_lo -= a_lo & -cBA;
b_hi -= a_hi & -(uint64_t)cBA;
qa -= pa & -(int64_t)cBA;
qb -= pb & -(int64_t)cBA;
/*
* Shifting.
*/
a_lo += a_lo & (cA - 1);
pa += pa & ((int64_t)cA - 1);
pb += pb & ((int64_t)cA - 1);
a_hi ^= (a_hi ^ (a_hi >> 1)) & -(uint64_t)cA;
b_lo += b_lo & -cA;
qa += qa & -(int64_t)cA;
qb += qb & -(int64_t)cA;
b_hi ^= (b_hi ^ (b_hi >> 1)) & ((uint64_t)cA - 1);
}
/*
* Replace a and b with new values a' and b'.
*/
r = co_reduce(a, b, len, pa, pb, qa, qb);
pa -= pa * ((r & 1) << 1);
pb -= pb * ((r & 1) << 1);
qa -= qa * (r & 2);
qb -= qb * (r & 2);
co_reduce_mod(u, v, len, pa, pb, qa, qb, m + 1, m0i);
}
/*
* Now one of the arrays should be 0, and the other contains
* the GCD. If a is 0, then u is 0 as well, and v contains
* the division result.
* Result is correct if and only if GCD is 1.
*/
r = (a[0] | b[0]) ^ 1;
u[0] |= v[0];
for (k = 1; k < len; k ++) {
r |= a[k] | b[k];
u[k] |= v[k];
}
return EQ0(r);
}

View file

@ -244,7 +244,7 @@ br_i62_modpow_opt(uint32_t *x31, const unsigned char *e, size_t elen,
mw62num = (mw31num + 1) >> 1;
/*
* In order to apply this function, we must have enough room tp
* In order to apply this function, we must have enough room to
* copy the operand and modulus into the temporary array, along
* with at least two temporaries. If there is not enough room,
* switch to br_i31_modpow(). We also use br_i31_modpow() if the
@ -477,3 +477,17 @@ br_i62_modpow_opt(uint32_t *x31, const unsigned char *e, size_t elen,
}
#endif
/* see inner.h */
uint32_t
br_i62_modpow_opt_as_i31(uint32_t *x31, const unsigned char *e, size_t elen,
const uint32_t *m31, uint32_t m0i31, uint32_t *tmp, size_t twlen)
{
/*
* As documented, this function expects the 'tmp' argument to be
* 64-bit aligned. This is OK since this function is internal (it
* is not part of BearSSL's public API).
*/
return br_i62_modpow_opt(x31, e, elen, m31, m0i31,
(uint64_t *)tmp, twlen >> 1);
}

View file

@ -0,0 +1,38 @@
/*
* Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include "inner.h"
/* see bearssl_rsa.h */
br_rsa_keygen
br_rsa_keygen_get_default(void)
{
#if BR_INT128 || BR_UMUL128
return &br_rsa_i62_keygen;
#elif BR_LOMUL
return &br_rsa_i15_keygen;
#else
return &br_rsa_i31_keygen;
#endif
}

583
src/rsa/rsa_i15_keygen.c Normal file
View file

@ -0,0 +1,583 @@
/*
* Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include "inner.h"
/*
* Make a random integer of the provided size. The size is encoded.
* The header word is untouched.
*/
static void
mkrand(const br_prng_class **rng, uint16_t *x, uint32_t esize)
{
size_t u, len;
unsigned m;
len = (esize + 15) >> 4;
(*rng)->generate(rng, x + 1, len * sizeof(uint16_t));
for (u = 1; u < len; u ++) {
x[u] &= 0x7FFF;
}
m = esize & 15;
if (m == 0) {
x[len] &= 0x7FFF;
} else {
x[len] &= 0x7FFF >> (15 - m);
}
}
/*
* This is the big-endian unsigned representation of the product of
* all small primes from 13 to 1481.
*/
static const unsigned char SMALL_PRIMES[] = {
0x2E, 0xAB, 0x92, 0xD1, 0x8B, 0x12, 0x47, 0x31, 0x54, 0x0A,
0x99, 0x5D, 0x25, 0x5E, 0xE2, 0x14, 0x96, 0x29, 0x1E, 0xB7,
0x78, 0x70, 0xCC, 0x1F, 0xA5, 0xAB, 0x8D, 0x72, 0x11, 0x37,
0xFB, 0xD8, 0x1E, 0x3F, 0x5B, 0x34, 0x30, 0x17, 0x8B, 0xE5,
0x26, 0x28, 0x23, 0xA1, 0x8A, 0xA4, 0x29, 0xEA, 0xFD, 0x9E,
0x39, 0x60, 0x8A, 0xF3, 0xB5, 0xA6, 0xEB, 0x3F, 0x02, 0xB6,
0x16, 0xC3, 0x96, 0x9D, 0x38, 0xB0, 0x7D, 0x82, 0x87, 0x0C,
0xF7, 0xBE, 0x24, 0xE5, 0x5F, 0x41, 0x04, 0x79, 0x76, 0x40,
0xE7, 0x00, 0x22, 0x7E, 0xB5, 0x85, 0x7F, 0x8D, 0x01, 0x50,
0xE9, 0xD3, 0x29, 0x42, 0x08, 0xB3, 0x51, 0x40, 0x7B, 0xD7,
0x8D, 0xCC, 0x10, 0x01, 0x64, 0x59, 0x28, 0xB6, 0x53, 0xF3,
0x50, 0x4E, 0xB1, 0xF2, 0x58, 0xCD, 0x6E, 0xF5, 0x56, 0x3E,
0x66, 0x2F, 0xD7, 0x07, 0x7F, 0x52, 0x4C, 0x13, 0x24, 0xDC,
0x8E, 0x8D, 0xCC, 0xED, 0x77, 0xC4, 0x21, 0xD2, 0xFD, 0x08,
0xEA, 0xD7, 0xC0, 0x5C, 0x13, 0x82, 0x81, 0x31, 0x2F, 0x2B,
0x08, 0xE4, 0x80, 0x04, 0x7A, 0x0C, 0x8A, 0x3C, 0xDC, 0x22,
0xE4, 0x5A, 0x7A, 0xB0, 0x12, 0x5E, 0x4A, 0x76, 0x94, 0x77,
0xC2, 0x0E, 0x92, 0xBA, 0x8A, 0xA0, 0x1F, 0x14, 0x51, 0x1E,
0x66, 0x6C, 0x38, 0x03, 0x6C, 0xC7, 0x4A, 0x4B, 0x70, 0x80,
0xAF, 0xCA, 0x84, 0x51, 0xD8, 0xD2, 0x26, 0x49, 0xF5, 0xA8,
0x5E, 0x35, 0x4B, 0xAC, 0xCE, 0x29, 0x92, 0x33, 0xB7, 0xA2,
0x69, 0x7D, 0x0C, 0xE0, 0x9C, 0xDB, 0x04, 0xD6, 0xB4, 0xBC,
0x39, 0xD7, 0x7F, 0x9E, 0x9D, 0x78, 0x38, 0x7F, 0x51, 0x54,
0x50, 0x8B, 0x9E, 0x9C, 0x03, 0x6C, 0xF5, 0x9D, 0x2C, 0x74,
0x57, 0xF0, 0x27, 0x2A, 0xC3, 0x47, 0xCA, 0xB9, 0xD7, 0x5C,
0xFF, 0xC2, 0xAC, 0x65, 0x4E, 0xBD
};
/*
* We need temporary values for at least 7 integers of the same size
* as a factor (including header word); more space helps with performance
* (in modular exponentiations), but we much prefer to remain under
* 2 kilobytes in total, to save stack space. The macro TEMPS below
* exceeds 1024 (which is a count in 16-bit words) when BR_MAX_RSA_SIZE
* is greater than 4350 (default value is 4096, so the 2-kB limit is
* maintained unless BR_MAX_RSA_SIZE was modified).
*/
#define MAX(x, y) ((x) > (y) ? (x) : (y))
#define TEMPS MAX(1024, 7 * ((((BR_MAX_RSA_SIZE + 1) >> 1) + 29) / 15))
/*
* Perform trial division on a candidate prime. This computes
* y = SMALL_PRIMES mod x, then tries to compute y/y mod x. The
* br_i15_moddiv() function will report an error if y is not invertible
* modulo x. Returned value is 1 on success (none of the small primes
* divides x), 0 on error (a non-trivial GCD is obtained).
*
* This function assumes that x is odd.
*/
static uint32_t
trial_divisions(const uint16_t *x, uint16_t *t)
{
uint16_t *y;
uint16_t x0i;
y = t;
t += 1 + ((x[0] + 15) >> 4);
x0i = br_i15_ninv15(x[1]);
br_i15_decode_reduce(y, SMALL_PRIMES, sizeof SMALL_PRIMES, x);
return br_i15_moddiv(y, y, x, x0i, t);
}
/*
* Perform n rounds of Miller-Rabin on the candidate prime x. This
* function assumes that x = 3 mod 4.
*
* Returned value is 1 on success (all rounds completed successfully),
* 0 otherwise.
*/
static uint32_t
miller_rabin(const br_prng_class **rng, const uint16_t *x, int n,
uint16_t *t, size_t tlen)
{
/*
* Since x = 3 mod 4, the Miller-Rabin test is simple:
* - get a random base a (such that 1 < a < x-1)
* - compute z = a^((x-1)/2) mod x
* - if z != 1 and z != x-1, the number x is composite
*
* We generate bases 'a' randomly with a size which is
* one bit less than x, which ensures that a < x-1. It
* is not useful to verify that a > 1 because the probability
* that we get a value a equal to 0 or 1 is much smaller
* than the probability of our Miller-Rabin tests not to
* detect a composite, which is already quite smaller than the
* probability of the hardware misbehaving and return a
* composite integer because of some glitch (e.g. bad RAM
* or ill-timed cosmic ray).
*/
unsigned char *xm1d2;
size_t xlen, xm1d2_len, xm1d2_len_u16, u;
uint32_t asize;
unsigned cc;
uint16_t x0i;
/*
* Compute (x-1)/2 (encoded).
*/
xm1d2 = (unsigned char *)t;
xm1d2_len = ((x[0] - (x[0] >> 4)) + 7) >> 3;
br_i15_encode(xm1d2, xm1d2_len, x);
cc = 0;
for (u = 0; u < xm1d2_len; u ++) {
unsigned w;
w = xm1d2[u];
xm1d2[u] = (unsigned char)((w >> 1) | cc);
cc = w << 7;
}
/*
* We used some words of the provided buffer for (x-1)/2.
*/
xm1d2_len_u16 = (xm1d2_len + 1) >> 1;
t += xm1d2_len_u16;
tlen -= xm1d2_len_u16;
xlen = (x[0] + 15) >> 4;
asize = x[0] - 1 - EQ0(x[0] & 15);
x0i = br_i15_ninv15(x[1]);
while (n -- > 0) {
uint16_t *a;
uint32_t eq1, eqm1;
/*
* Generate a random base. We don't need the base to be
* really uniform modulo x, so we just get a random
* number which is one bit shorter than x.
*/
a = t;
a[0] = x[0];
a[xlen] = 0;
mkrand(rng, a, asize);
/*
* Compute a^((x-1)/2) mod x. We assume here that the
* function will not fail (the temporary array is large
* enough).
*/
br_i15_modpow_opt(a, xm1d2, xm1d2_len,
x, x0i, t + 1 + xlen, tlen - 1 - xlen);
/*
* We must obtain either 1 or x-1. Note that x is odd,
* hence x-1 differs from x only in its low word (no
* carry).
*/
eq1 = a[1] ^ 1;
eqm1 = a[1] ^ (x[1] - 1);
for (u = 2; u <= xlen; u ++) {
eq1 |= a[u];
eqm1 |= a[u] ^ x[u];
}
if ((EQ0(eq1) | EQ0(eqm1)) == 0) {
return 0;
}
}
return 1;
}
/*
* Create a random prime of the provided size. 'size' is the _encoded_
* bit length. The two top bits and the two bottom bits are set to 1.
*/
static void
mkprime(const br_prng_class **rng, uint16_t *x, uint32_t esize,
uint32_t pubexp, uint16_t *t, size_t tlen)
{
size_t len;
x[0] = esize;
len = (esize + 15) >> 4;
for (;;) {
size_t u;
uint32_t m3, m5, m7, m11;
int rounds;
/*
* Generate random bits. We force the two top bits and the
* two bottom bits to 1.
*/
mkrand(rng, x, esize);
if ((esize & 15) == 0) {
x[len] |= 0x6000;
} else if ((esize & 15) == 1) {
x[len] |= 0x0001;
x[len - 1] |= 0x4000;
} else {
x[len] |= 0x0003 << ((esize & 15) - 2);
}
x[1] |= 0x0003;
/*
* Trial division with low primes (3, 5, 7 and 11). We
* use the following properties:
*
* 2^2 = 1 mod 3
* 2^4 = 1 mod 5
* 2^3 = 1 mod 7
* 2^10 = 1 mod 11
*/
m3 = 0;
m5 = 0;
m7 = 0;
m11 = 0;
for (u = 0; u < len; u ++) {
uint32_t w;
w = x[1 + u];
m3 += w << (u & 1);
m3 = (m3 & 0xFF) + (m3 >> 8);
m5 += w << ((4 - u) & 3);
m5 = (m5 & 0xFF) + (m5 >> 8);
m7 += w;
m7 = (m7 & 0x1FF) + (m7 >> 9);
m11 += w << (5 & -(u & 1));
m11 = (m11 & 0x3FF) + (m11 >> 10);
}
/*
* Maximum values of m* at this point:
* m3: 511
* m5: 2310
* m7: 510
* m11: 2047
* We use the same properties to make further reductions.
*/
m3 = (m3 & 0x0F) + (m3 >> 4); /* max: 46 */
m3 = (m3 & 0x0F) + (m3 >> 4); /* max: 16 */
m3 = ((m3 * 43) >> 5) & 3;
m5 = (m5 & 0xFF) + (m5 >> 8); /* max: 263 */
m5 = (m5 & 0x0F) + (m5 >> 4); /* max: 30 */
m5 = (m5 & 0x0F) + (m5 >> 4); /* max: 15 */
m5 -= 10 & -GT(m5, 9);
m5 -= 5 & -GT(m5, 4);
m7 = (m7 & 0x3F) + (m7 >> 6); /* max: 69 */
m7 = (m7 & 7) + (m7 >> 3); /* max: 14 */
m7 = ((m7 * 147) >> 7) & 7;
/*
* 2^5 = 32 = -1 mod 11.
*/
m11 = (m11 & 0x1F) + 66 - (m11 >> 5); /* max: 97 */
m11 -= 88 & -GT(m11, 87);
m11 -= 44 & -GT(m11, 43);
m11 -= 22 & -GT(m11, 21);
m11 -= 11 & -GT(m11, 10);
/*
* If any of these modulo is 0, then the candidate is
* not prime. Also, if pubexp is 3, 5, 7 or 11, and the
* corresponding modulus is 1, then the candidate must
* be rejected, because we need e to be invertible
* modulo p-1. We can use simple comparisons here
* because they won't leak information on a candidate
* that we keep, only on one that we reject (and is thus
* not secret).
*/
if (m3 == 0 || m5 == 0 || m7 == 0 || m11 == 0) {
continue;
}
if ((pubexp == 3 && m3 == 1)
|| (pubexp == 5 && m5 == 5)
|| (pubexp == 7 && m5 == 7)
|| (pubexp == 11 && m5 == 11))
{
continue;
}
/*
* More trial divisions.
*/
if (!trial_divisions(x, t)) {
continue;
}
/*
* Miller-Rabin algorithm. Since we selected a random
* integer, not a maliciously crafted integer, we can use
* relatively few rounds to lower the risk of a false
* positive (i.e. declaring prime a non-prime) under
* 2^(-80). It is not useful to lower the probability much
* below that, since that would be substantially below
* the probability of the hardware misbehaving. Sufficient
* numbers of rounds are extracted from the Handbook of
* Applied Cryptography, note 4.49 (page 149).
*
* Since we work on the encoded size (esize), we need to
* compare with encoded thresholds.
*/
if (esize < 320) {
rounds = 12;
} else if (esize < 480) {
rounds = 9;
} else if (esize < 693) {
rounds = 6;
} else if (esize < 906) {
rounds = 4;
} else if (esize < 1386) {
rounds = 3;
} else {
rounds = 2;
}
if (miller_rabin(rng, x, rounds, t, tlen)) {
return;
}
}
}
/*
* Let p be a prime (p > 2^33, p = 3 mod 4). Let m = (p-1)/2, provided
* as parameter (with announced bit length equal to that of p). This
* function computes d = 1/e mod p-1 (for an odd integer e). Returned
* value is 1 on success, 0 on error (an error is reported if e is not
* invertible modulo p-1).
*
* The temporary buffer (t) must have room for at least 4 integers of
* the size of p.
*/
static uint32_t
invert_pubexp(uint16_t *d, const uint16_t *m, uint32_t e, uint16_t *t)
{
uint16_t *f;
uint32_t r;
f = t;
t += 1 + ((m[0] + 15) >> 4);
/*
* Compute d = 1/e mod m. Since p = 3 mod 4, m is odd.
*/
br_i15_zero(d, m[0]);
d[1] = 1;
br_i15_zero(f, m[0]);
f[1] = e & 0x7FFF;
f[2] = (e >> 15) & 0x7FFF;
f[3] = e >> 30;
r = br_i15_moddiv(d, f, m, br_i15_ninv15(m[1]), t);
/*
* We really want d = 1/e mod p-1, with p = 2m. By the CRT,
* the result is either the d we got, or d + m.
*
* Let's write e*d = 1 + k*m, for some integer k. Integers e
* and m are odd. If d is odd, then e*d is odd, which implies
* that k must be even; in that case, e*d = 1 + (k/2)*2m, and
* thus d is already fine. Conversely, if d is even, then k
* is odd, and we must add m to d in order to get the correct
* result.
*/
br_i15_add(d, m, (uint32_t)(1 - (d[1] & 1)));
return r;
}
/*
* Swap two buffers in RAM. They must be disjoint.
*/
static void
bufswap(void *b1, void *b2, size_t len)
{
size_t u;
unsigned char *buf1, *buf2;
buf1 = b1;
buf2 = b2;
for (u = 0; u < len; u ++) {
unsigned w;
w = buf1[u];
buf1[u] = buf2[u];
buf2[u] = w;
}
}
/* see bearssl_rsa.h */
uint32_t
br_rsa_i15_keygen(const br_prng_class **rng,
br_rsa_private_key *sk, unsigned char *kbuf_priv,
br_rsa_public_key *pk, unsigned char *kbuf_pub,
unsigned size, uint32_t pubexp)
{
uint32_t esize_p, esize_q;
size_t plen, qlen, tlen;
uint16_t *p, *q, *t;
uint16_t tmp[TEMPS];
uint32_t r;
if (size < BR_MIN_RSA_SIZE || size > BR_MAX_RSA_SIZE) {
return 0;
}
if (pubexp == 0) {
pubexp = 3;
} else if (pubexp == 1 || (pubexp & 1) == 0) {
return 0;
}
esize_p = (size + 1) >> 1;
esize_q = size - esize_p;
sk->n_bitlen = size;
sk->p = kbuf_priv;
sk->plen = (esize_p + 7) >> 3;
sk->q = sk->p + sk->plen;
sk->qlen = (esize_q + 7) >> 3;
sk->dp = sk->q + sk->qlen;
sk->dplen = sk->plen;
sk->dq = sk->dp + sk->dplen;
sk->dqlen = sk->qlen;
sk->iq = sk->dq + sk->dqlen;
sk->iqlen = sk->plen;
if (pk != NULL) {
pk->n = kbuf_pub;
pk->nlen = (size + 7) >> 3;
pk->e = pk->n + pk->nlen;
pk->elen = 4;
br_enc32be(pk->e, pubexp);
while (*pk->e == 0) {
pk->e ++;
pk->elen --;
}
}
/*
* We now switch to encoded sizes.
*
* floor((x * 17477) / (2^18)) is equal to floor(x/15) for all
* integers x from 0 to 23833.
*/
esize_p += MUL15(esize_p, 17477) >> 18;
esize_q += MUL15(esize_q, 17477) >> 18;
plen = (esize_p + 15) >> 4;
qlen = (esize_q + 15) >> 4;
p = tmp;
q = p + 1 + plen;
t = q + 1 + qlen;
tlen = ((sizeof tmp) / sizeof(uint16_t)) - (2 + plen + qlen);
/*
* When looking for primes p and q, we temporarily divide
* candidates by 2, in order to compute the inverse of the
* public exponent.
*/
for (;;) {
mkprime(rng, p, esize_p, pubexp, t, tlen);
br_i15_rshift(p, 1);
if (invert_pubexp(t, p, pubexp, t + 1 + plen)) {
br_i15_add(p, p, 1);
p[1] |= 1;
br_i15_encode(sk->p, sk->plen, p);
br_i15_encode(sk->dp, sk->dplen, t);
break;
}
}
for (;;) {
mkprime(rng, q, esize_q, pubexp, t, tlen);
br_i15_rshift(q, 1);
if (invert_pubexp(t, q, pubexp, t + 1 + qlen)) {
br_i15_add(q, q, 1);
q[1] |= 1;
br_i15_encode(sk->q, sk->qlen, q);
br_i15_encode(sk->dq, sk->dqlen, t);
break;
}
}
/*
* If p and q have the same size, then it is possible that q > p
* (when the target modulus size is odd, we generate p with a
* greater bit length than q). If q > p, we want to swap p and q
* (and also dp and dq) for two reasons:
* - The final step below (inversion of q modulo p) is easier if
* p > q.
* - While BearSSL's RSA code is perfectly happy with RSA keys such
* that p < q, some other implementations have restrictions and
* require p > q.
*
* Note that we can do a simple non-constant-time swap here,
* because the only information we leak here is that we insist on
* returning p and q such that p > q, which is not a secret.
*/
if (esize_p == esize_q && br_i15_sub(p, q, 0) == 1) {
bufswap(p, q, (1 + plen) * sizeof *p);
bufswap(sk->p, sk->q, sk->plen);
bufswap(sk->dp, sk->dq, sk->dplen);
}
/*
* We have produced p, q, dp and dq. We can now compute iq = 1/d mod p.
*
* We ensured that p >= q, so this is just a matter of updating the
* header word for q (and possibly adding an extra word).
*
* Theoretically, the call below may fail, in case we were
* extraordinarily unlucky, and p = q. Another failure case is if
* Miller-Rabin failed us _twice_, and p and q are non-prime and
* have a factor is common. We report the error mostly because it
* is cheap and we can, but in practice this never happens (or, at
* least, it happens way less often than hardware glitches).
*/
q[0] = p[0];
if (plen > qlen) {
q[plen] = 0;
t ++;
tlen --;
}
br_i15_zero(t, p[0]);
t[1] = 1;
r = br_i15_moddiv(t, q, p, br_i15_ninv15(p[1]), t + 1 + plen);
br_i15_encode(sk->iq, sk->iqlen, t);
/*
* Compute the public modulus too, if required.
*/
if (pk != NULL) {
br_i15_zero(t, p[0]);
br_i15_mulacc(t, p, q);
br_i15_encode(pk->n, pk->nlen, t);
}
return r;
}

37
src/rsa/rsa_i31_keygen.c Normal file
View file

@ -0,0 +1,37 @@
/*
* Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include "inner.h"
/* see bearssl_rsa.h */
uint32_t
br_rsa_i31_keygen(const br_prng_class **rng,
br_rsa_private_key *sk, unsigned char *kbuf_priv,
br_rsa_public_key *pk, unsigned char *kbuf_pub,
unsigned size, uint32_t pubexp)
{
return br_rsa_i31_keygen_inner(rng,
sk, kbuf_priv, pk, kbuf_pub, size, pubexp,
&br_i31_modpow_opt);
}

View file

@ -0,0 +1,608 @@
/*
* Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include "inner.h"
/*
* Make a random integer of the provided size. The size is encoded.
* The header word is untouched.
*/
static void
mkrand(const br_prng_class **rng, uint32_t *x, uint32_t esize)
{
size_t u, len;
unsigned m;
len = (esize + 31) >> 5;
(*rng)->generate(rng, x + 1, len * sizeof(uint32_t));
for (u = 1; u < len; u ++) {
x[u] &= 0x7FFFFFFF;
}
m = esize & 31;
if (m == 0) {
x[len] &= 0x7FFFFFFF;
} else {
x[len] &= 0x7FFFFFFF >> (31 - m);
}
}
/*
* This is the big-endian unsigned representation of the product of
* all small primes from 13 to 1481.
*/
static const unsigned char SMALL_PRIMES[] = {
0x2E, 0xAB, 0x92, 0xD1, 0x8B, 0x12, 0x47, 0x31, 0x54, 0x0A,
0x99, 0x5D, 0x25, 0x5E, 0xE2, 0x14, 0x96, 0x29, 0x1E, 0xB7,
0x78, 0x70, 0xCC, 0x1F, 0xA5, 0xAB, 0x8D, 0x72, 0x11, 0x37,
0xFB, 0xD8, 0x1E, 0x3F, 0x5B, 0x34, 0x30, 0x17, 0x8B, 0xE5,
0x26, 0x28, 0x23, 0xA1, 0x8A, 0xA4, 0x29, 0xEA, 0xFD, 0x9E,
0x39, 0x60, 0x8A, 0xF3, 0xB5, 0xA6, 0xEB, 0x3F, 0x02, 0xB6,
0x16, 0xC3, 0x96, 0x9D, 0x38, 0xB0, 0x7D, 0x82, 0x87, 0x0C,
0xF7, 0xBE, 0x24, 0xE5, 0x5F, 0x41, 0x04, 0x79, 0x76, 0x40,
0xE7, 0x00, 0x22, 0x7E, 0xB5, 0x85, 0x7F, 0x8D, 0x01, 0x50,
0xE9, 0xD3, 0x29, 0x42, 0x08, 0xB3, 0x51, 0x40, 0x7B, 0xD7,
0x8D, 0xCC, 0x10, 0x01, 0x64, 0x59, 0x28, 0xB6, 0x53, 0xF3,
0x50, 0x4E, 0xB1, 0xF2, 0x58, 0xCD, 0x6E, 0xF5, 0x56, 0x3E,
0x66, 0x2F, 0xD7, 0x07, 0x7F, 0x52, 0x4C, 0x13, 0x24, 0xDC,
0x8E, 0x8D, 0xCC, 0xED, 0x77, 0xC4, 0x21, 0xD2, 0xFD, 0x08,
0xEA, 0xD7, 0xC0, 0x5C, 0x13, 0x82, 0x81, 0x31, 0x2F, 0x2B,
0x08, 0xE4, 0x80, 0x04, 0x7A, 0x0C, 0x8A, 0x3C, 0xDC, 0x22,
0xE4, 0x5A, 0x7A, 0xB0, 0x12, 0x5E, 0x4A, 0x76, 0x94, 0x77,
0xC2, 0x0E, 0x92, 0xBA, 0x8A, 0xA0, 0x1F, 0x14, 0x51, 0x1E,
0x66, 0x6C, 0x38, 0x03, 0x6C, 0xC7, 0x4A, 0x4B, 0x70, 0x80,
0xAF, 0xCA, 0x84, 0x51, 0xD8, 0xD2, 0x26, 0x49, 0xF5, 0xA8,
0x5E, 0x35, 0x4B, 0xAC, 0xCE, 0x29, 0x92, 0x33, 0xB7, 0xA2,
0x69, 0x7D, 0x0C, 0xE0, 0x9C, 0xDB, 0x04, 0xD6, 0xB4, 0xBC,
0x39, 0xD7, 0x7F, 0x9E, 0x9D, 0x78, 0x38, 0x7F, 0x51, 0x54,
0x50, 0x8B, 0x9E, 0x9C, 0x03, 0x6C, 0xF5, 0x9D, 0x2C, 0x74,
0x57, 0xF0, 0x27, 0x2A, 0xC3, 0x47, 0xCA, 0xB9, 0xD7, 0x5C,
0xFF, 0xC2, 0xAC, 0x65, 0x4E, 0xBD
};
/*
* We need temporary values for at least 7 integers of the same size
* as a factor (including header word); more space helps with performance
* (in modular exponentiations), but we much prefer to remain under
* 2 kilobytes in total, to save stack space. The macro TEMPS below
* exceeds 512 (which is a count in 32-bit words) when BR_MAX_RSA_SIZE
* is greater than 4464 (default value is 4096, so the 2-kB limit is
* maintained unless BR_MAX_RSA_SIZE was modified).
*/
#define MAX(x, y) ((x) > (y) ? (x) : (y))
#define ROUND2(x) ((((x) + 1) >> 1) << 1)
#define TEMPS MAX(512, ROUND2(7 * ((((BR_MAX_RSA_SIZE + 1) >> 1) + 61) / 31)))
/*
* Perform trial division on a candidate prime. This computes
* y = SMALL_PRIMES mod x, then tries to compute y/y mod x. The
* br_i31_moddiv() function will report an error if y is not invertible
* modulo x. Returned value is 1 on success (none of the small primes
* divides x), 0 on error (a non-trivial GCD is obtained).
*
* This function assumes that x is odd.
*/
static uint32_t
trial_divisions(const uint32_t *x, uint32_t *t)
{
uint32_t *y;
uint32_t x0i;
y = t;
t += 1 + ((x[0] + 31) >> 5);
x0i = br_i31_ninv31(x[1]);
br_i31_decode_reduce(y, SMALL_PRIMES, sizeof SMALL_PRIMES, x);
return br_i31_moddiv(y, y, x, x0i, t);
}
/*
* Perform n rounds of Miller-Rabin on the candidate prime x. This
* function assumes that x = 3 mod 4.
*
* Returned value is 1 on success (all rounds completed successfully),
* 0 otherwise.
*/
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)
{
/*
* Since x = 3 mod 4, the Miller-Rabin test is simple:
* - get a random base a (such that 1 < a < x-1)
* - compute z = a^((x-1)/2) mod x
* - if z != 1 and z != x-1, the number x is composite
*
* We generate bases 'a' randomly with a size which is
* one bit less than x, which ensures that a < x-1. It
* is not useful to verify that a > 1 because the probability
* that we get a value a equal to 0 or 1 is much smaller
* than the probability of our Miller-Rabin tests not to
* detect a composite, which is already quite smaller than the
* probability of the hardware misbehaving and return a
* composite integer because of some glitch (e.g. bad RAM
* or ill-timed cosmic ray).
*/
unsigned char *xm1d2;
size_t xlen, xm1d2_len, xm1d2_len_u32, u;
uint32_t asize;
unsigned cc;
uint32_t x0i;
/*
* Compute (x-1)/2 (encoded).
*/
xm1d2 = (unsigned char *)t;
xm1d2_len = ((x[0] - (x[0] >> 5)) + 7) >> 3;
br_i31_encode(xm1d2, xm1d2_len, x);
cc = 0;
for (u = 0; u < xm1d2_len; u ++) {
unsigned w;
w = xm1d2[u];
xm1d2[u] = (unsigned char)((w >> 1) | cc);
cc = w << 7;
}
/*
* We used some words of the provided buffer for (x-1)/2.
*/
xm1d2_len_u32 = (xm1d2_len + 3) >> 2;
t += xm1d2_len_u32;
tlen -= xm1d2_len_u32;
xlen = (x[0] + 31) >> 5;
asize = x[0] - 1 - EQ0(x[0] & 31);
x0i = br_i31_ninv31(x[1]);
while (n -- > 0) {
uint32_t *a, *t2;
uint32_t eq1, eqm1;
size_t t2len;
/*
* Generate a random base. We don't need the base to be
* really uniform modulo x, so we just get a random
* number which is one bit shorter than x.
*/
a = t;
a[0] = x[0];
a[xlen] = 0;
mkrand(rng, a, asize);
/*
* Compute a^((x-1)/2) mod x. We assume here that the
* function will not fail (the temporary array is large
* enough).
*/
t2 = t + 1 + xlen;
t2len = tlen - 1 - xlen;
if ((t2len & 1) != 0) {
/*
* Since the source array is 64-bit aligned and
* has an even number of elements (TEMPS), we
* can use the parity of the remaining length to
* detect and adjust alignment.
*/
t2 ++;
t2len --;
}
mp31(a, xm1d2, xm1d2_len, x, x0i, t2, t2len);
/*
* We must obtain either 1 or x-1. Note that x is odd,
* hence x-1 differs from x only in its low word (no
* carry).
*/
eq1 = a[1] ^ 1;
eqm1 = a[1] ^ (x[1] - 1);
for (u = 2; u <= xlen; u ++) {
eq1 |= a[u];
eqm1 |= a[u] ^ x[u];
}
if ((EQ0(eq1) | EQ0(eqm1)) == 0) {
return 0;
}
}
return 1;
}
/*
* Create a random prime of the provided size. 'size' is the _encoded_
* bit length. The two top bits and the two bottom bits are set to 1.
*/
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)
{
size_t len;
x[0] = esize;
len = (esize + 31) >> 5;
for (;;) {
size_t u;
uint32_t m3, m5, m7, m11;
int rounds, s7, s11;
/*
* Generate random bits. We force the two top bits and the
* two bottom bits to 1.
*/
mkrand(rng, x, esize);
if ((esize & 31) == 0) {
x[len] |= 0x60000000;
} else if ((esize & 31) == 1) {
x[len] |= 0x00000001;
x[len - 1] |= 0x40000000;
} else {
x[len] |= 0x00000003 << ((esize & 31) - 2);
}
x[1] |= 0x00000003;
/*
* Trial division with low primes (3, 5, 7 and 11). We
* use the following properties:
*
* 2^2 = 1 mod 3
* 2^4 = 1 mod 5
* 2^3 = 1 mod 7
* 2^10 = 1 mod 11
*/
m3 = 0;
m5 = 0;
m7 = 0;
m11 = 0;
s7 = 0;
s11 = 0;
for (u = 0; u < len; u ++) {
uint32_t w, w3, w5, w7, w11;
w = x[1 + u];
w3 = (w & 0xFFFF) + (w >> 16); /* max: 98302 */
w5 = (w & 0xFFFF) + (w >> 16); /* max: 98302 */
w7 = (w & 0x7FFF) + (w >> 15); /* max: 98302 */
w11 = (w & 0xFFFFF) + (w >> 20); /* max: 1050622 */
m3 += w3 << (u & 1);
m3 = (m3 & 0xFF) + (m3 >> 8); /* max: 1025 */
m5 += w5 << ((4 - u) & 3);
m5 = (m5 & 0xFFF) + (m5 >> 12); /* max: 4479 */
m7 += w7 << s7;
m7 = (m7 & 0x1FF) + (m7 >> 9); /* max: 1280 */
if (++ s7 == 3) {
s7 = 0;
}
m11 += w11 << s11;
if (++ s11 == 10) {
s11 = 0;
}
m11 = (m11 & 0x3FF) + (m11 >> 10); /* max: 526847 */
}
m3 = (m3 & 0x3F) + (m3 >> 6); /* max: 78 */
m3 = (m3 & 0x0F) + (m3 >> 4); /* max: 18 */
m3 = ((m3 * 43) >> 5) & 3;
m5 = (m5 & 0xFF) + (m5 >> 8); /* max: 271 */
m5 = (m5 & 0x0F) + (m5 >> 4); /* max: 31 */
m5 -= 20 & -GT(m5, 19);
m5 -= 10 & -GT(m5, 9);
m5 -= 5 & -GT(m5, 4);
m7 = (m7 & 0x3F) + (m7 >> 6); /* max: 82 */
m7 = (m7 & 0x07) + (m7 >> 3); /* max: 16 */
m7 = ((m7 * 147) >> 7) & 7;
/*
* 2^5 = 32 = -1 mod 11.
*/
m11 = (m11 & 0x3FF) + (m11 >> 10); /* max: 1536 */
m11 = (m11 & 0x3FF) + (m11 >> 10); /* max: 1023 */
m11 = (m11 & 0x1F) + 33 - (m11 >> 5); /* max: 64 */
m11 -= 44 & -GT(m11, 43);
m11 -= 22 & -GT(m11, 21);
m11 -= 11 & -GT(m11, 10);
/*
* If any of these modulo is 0, then the candidate is
* not prime. Also, if pubexp is 3, 5, 7 or 11, and the
* corresponding modulus is 1, then the candidate must
* be rejected, because we need e to be invertible
* modulo p-1. We can use simple comparisons here
* because they won't leak information on a candidate
* that we keep, only on one that we reject (and is thus
* not secret).
*/
if (m3 == 0 || m5 == 0 || m7 == 0 || m11 == 0) {
continue;
}
if ((pubexp == 3 && m3 == 1)
|| (pubexp == 5 && m5 == 5)
|| (pubexp == 7 && m5 == 7)
|| (pubexp == 11 && m5 == 11))
{
continue;
}
/*
* More trial divisions.
*/
if (!trial_divisions(x, t)) {
continue;
}
/*
* Miller-Rabin algorithm. Since we selected a random
* integer, not a maliciously crafted integer, we can use
* relatively few rounds to lower the risk of a false
* positive (i.e. declaring prime a non-prime) under
* 2^(-80). It is not useful to lower the probability much
* below that, since that would be substantially below
* the probability of the hardware misbehaving. Sufficient
* numbers of rounds are extracted from the Handbook of
* Applied Cryptography, note 4.49 (page 149).
*
* Since we work on the encoded size (esize), we need to
* compare with encoded thresholds.
*/
if (esize < 309) {
rounds = 12;
} else if (esize < 464) {
rounds = 9;
} else if (esize < 670) {
rounds = 6;
} else if (esize < 877) {
rounds = 4;
} else if (esize < 1341) {
rounds = 3;
} else {
rounds = 2;
}
if (miller_rabin(rng, x, rounds, t, tlen, mp31)) {
return;
}
}
}
/*
* Let p be a prime (p > 2^33, p = 3 mod 4). Let m = (p-1)/2, provided
* as parameter (with announced bit length equal to that of p). This
* function computes d = 1/e mod p-1 (for an odd integer e). Returned
* value is 1 on success, 0 on error (an error is reported if e is not
* invertible modulo p-1).
*
* The temporary buffer (t) must have room for at least 4 integers of
* the size of p.
*/
static uint32_t
invert_pubexp(uint32_t *d, const uint32_t *m, uint32_t e, uint32_t *t)
{
uint32_t *f;
uint32_t r;
f = t;
t += 1 + ((m[0] + 31) >> 5);
/*
* Compute d = 1/e mod m. Since p = 3 mod 4, m is odd.
*/
br_i31_zero(d, m[0]);
d[1] = 1;
br_i31_zero(f, m[0]);
f[1] = e & 0x7FFFFFFF;
f[2] = e >> 31;
r = br_i31_moddiv(d, f, m, br_i31_ninv31(m[1]), t);
/*
* We really want d = 1/e mod p-1, with p = 2m. By the CRT,
* the result is either the d we got, or d + m.
*
* Let's write e*d = 1 + k*m, for some integer k. Integers e
* and m are odd. If d is odd, then e*d is odd, which implies
* that k must be even; in that case, e*d = 1 + (k/2)*2m, and
* thus d is already fine. Conversely, if d is even, then k
* is odd, and we must add m to d in order to get the correct
* result.
*/
br_i31_add(d, m, (uint32_t)(1 - (d[1] & 1)));
return r;
}
/*
* Swap two buffers in RAM. They must be disjoint.
*/
static void
bufswap(void *b1, void *b2, size_t len)
{
size_t u;
unsigned char *buf1, *buf2;
buf1 = b1;
buf2 = b2;
for (u = 0; u < len; u ++) {
unsigned w;
w = buf1[u];
buf1[u] = buf2[u];
buf2[u] = w;
}
}
/* see inner.h */
uint32_t
br_rsa_i31_keygen_inner(const br_prng_class **rng,
br_rsa_private_key *sk, unsigned char *kbuf_priv,
br_rsa_public_key *pk, unsigned char *kbuf_pub,
unsigned size, uint32_t pubexp, br_i31_modpow_opt_type mp31)
{
uint32_t esize_p, esize_q;
size_t plen, qlen, tlen;
uint32_t *p, *q, *t;
union {
uint32_t t32[TEMPS];
uint64_t t64[TEMPS >> 1]; /* for 64-bit alignment */
} tmp;
uint32_t r;
if (size < BR_MIN_RSA_SIZE || size > BR_MAX_RSA_SIZE) {
return 0;
}
if (pubexp == 0) {
pubexp = 3;
} else if (pubexp == 1 || (pubexp & 1) == 0) {
return 0;
}
esize_p = (size + 1) >> 1;
esize_q = size - esize_p;
sk->n_bitlen = size;
sk->p = kbuf_priv;
sk->plen = (esize_p + 7) >> 3;
sk->q = sk->p + sk->plen;
sk->qlen = (esize_q + 7) >> 3;
sk->dp = sk->q + sk->qlen;
sk->dplen = sk->plen;
sk->dq = sk->dp + sk->dplen;
sk->dqlen = sk->qlen;
sk->iq = sk->dq + sk->dqlen;
sk->iqlen = sk->plen;
if (pk != NULL) {
pk->n = kbuf_pub;
pk->nlen = (size + 7) >> 3;
pk->e = pk->n + pk->nlen;
pk->elen = 4;
br_enc32be(pk->e, pubexp);
while (*pk->e == 0) {
pk->e ++;
pk->elen --;
}
}
/*
* We now switch to encoded sizes.
*
* floor((x * 16913) / (2^19)) is equal to floor(x/31) for all
* integers x from 0 to 34966; the intermediate product fits on
* 30 bits, thus we can use MUL31().
*/
esize_p += MUL31(esize_p, 16913) >> 19;
esize_q += MUL31(esize_q, 16913) >> 19;
plen = (esize_p + 31) >> 5;
qlen = (esize_q + 31) >> 5;
p = tmp.t32;
q = p + 1 + plen;
t = q + 1 + qlen;
tlen = ((sizeof tmp.t32) / sizeof(uint32_t)) - (2 + plen + qlen);
/*
* When looking for primes p and q, we temporarily divide
* candidates by 2, in order to compute the inverse of the
* public exponent.
*/
for (;;) {
mkprime(rng, p, esize_p, pubexp, t, tlen, mp31);
br_i31_rshift(p, 1);
if (invert_pubexp(t, p, pubexp, t + 1 + plen)) {
br_i31_add(p, p, 1);
p[1] |= 1;
br_i31_encode(sk->p, sk->plen, p);
br_i31_encode(sk->dp, sk->dplen, t);
break;
}
}
for (;;) {
mkprime(rng, q, esize_q, pubexp, t, tlen, mp31);
br_i31_rshift(q, 1);
if (invert_pubexp(t, q, pubexp, t + 1 + qlen)) {
br_i31_add(q, q, 1);
q[1] |= 1;
br_i31_encode(sk->q, sk->qlen, q);
br_i31_encode(sk->dq, sk->dqlen, t);
break;
}
}
/*
* If p and q have the same size, then it is possible that q > p
* (when the target modulus size is odd, we generate p with a
* greater bit length than q). If q > p, we want to swap p and q
* (and also dp and dq) for two reasons:
* - The final step below (inversion of q modulo p) is easier if
* p > q.
* - While BearSSL's RSA code is perfectly happy with RSA keys such
* that p < q, some other implementations have restrictions and
* require p > q.
*
* Note that we can do a simple non-constant-time swap here,
* because the only information we leak here is that we insist on
* returning p and q such that p > q, which is not a secret.
*/
if (esize_p == esize_q && br_i31_sub(p, q, 0) == 1) {
bufswap(p, q, (1 + plen) * sizeof *p);
bufswap(sk->p, sk->q, sk->plen);
bufswap(sk->dp, sk->dq, sk->dplen);
}
/*
* We have produced p, q, dp and dq. We can now compute iq = 1/d mod p.
*
* We ensured that p >= q, so this is just a matter of updating the
* header word for q (and possibly adding an extra word).
*
* Theoretically, the call below may fail, in case we were
* extraordinarily unlucky, and p = q. Another failure case is if
* Miller-Rabin failed us _twice_, and p and q are non-prime and
* have a factor is common. We report the error mostly because it
* is cheap and we can, but in practice this never happens (or, at
* least, it happens way less often than hardware glitches).
*/
q[0] = p[0];
if (plen > qlen) {
q[plen] = 0;
t ++;
tlen --;
}
br_i31_zero(t, p[0]);
t[1] = 1;
r = br_i31_moddiv(t, q, p, br_i31_ninv31(p[1]), t + 1 + plen);
br_i31_encode(sk->iq, sk->iqlen, t);
/*
* Compute the public modulus too, if required.
*/
if (pk != NULL) {
br_i31_zero(t, p[0]);
br_i31_mulacc(t, p, q);
br_i31_encode(pk->n, pk->nlen, t);
}
return r;
}

57
src/rsa/rsa_i62_keygen.c Normal file
View file

@ -0,0 +1,57 @@
/*
* Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include "inner.h"
#if BR_INT128 || BR_UMUL128
/* see bearssl_rsa.h */
uint32_t
br_rsa_i62_keygen(const br_prng_class **rng,
br_rsa_private_key *sk, unsigned char *kbuf_priv,
br_rsa_public_key *pk, unsigned char *kbuf_pub,
unsigned size, uint32_t pubexp)
{
return br_rsa_i31_keygen_inner(rng,
sk, kbuf_priv, pk, kbuf_pub, size, pubexp,
&br_i62_modpow_opt_as_i31);
}
/* see bearssl_rsa.h */
br_rsa_keygen
br_rsa_i62_keygen_get()
{
return &br_rsa_i62_keygen;
}
#else
/* see bearssl_rsa.h */
br_rsa_keygen
br_rsa_i62_keygen_get()
{
return 0;
}
#endif