~apreiml/hare

bf745a78ede8903cce00446fbf3c5f6b42b13225 — Armin Preiml 7 months ago 8afe9b2 dev
crypto::rsa: keygen wip

Signed-off-by: Armin Preiml <apreiml@strohwolke.at>
1 files changed, 346 insertions(+), 0 deletions(-)

A crypto/rsa/keygen.ha
A crypto/rsa/keygen.ha => crypto/rsa/keygen.ha +346 -0
@@ 0,0 1,346 @@
// SPDX-License-Identifier: MPL-2.0
// (c) Hare authors <https://harelang.org>

// The following code was initially ported from BearSSL.
//
// Copyright (c) 2017 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.
use crypto::bigint::*;
use crypto::math::*;
use io;

use crypto::random; // DEBUG
use encoding::hex;
use fmt;
use strings;

// Make a random integer of the provided size. The size is encoded. The header
// word is untouched.
fn mkrand(rnd: *io::stream, x: []word, xsz: word) void = {
	const l = (xsz + 31) >> 5;
	io::readall(rnd, (x[1..l + 1]: *[*]u8)[..l * size(word)])!;

	for (let i = 1z; i < l; i += 1) {
		x[i] &= 0x7FFFFFFF;
	};

	const m = xsz & 31;
	if (m == 0) {
		x[l] &= 0x7FFFFFFF;
	} else {
		x[l] &= 0x7FFFFFFF >> (31 - m);
	};
};

// This is the big-endian unsigned representation of the product of small small
// primes from 13 to 1481.
const SMALL_PRIMES: [_]u8 = [
	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,
];

// 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.
fn trial_divisions(x: []word, t: []word) word = {
	const x0i = ninv(x[1]);
	const tstart = 1 + ((x[0] + 31) >> 5);
	let y = t[tstart..];

	fmt::println("=========")!;
	dump(x);
	fmt::println("small", hex::encodestr(SMALL_PRIMES))!;
	encodereduce(y, SMALL_PRIMES, x);
	dump(y);
	fmt::println("---------")!;
	return moddiv(y, y, x, x0i, t);
};

type modpow_opt = fn () void;

// 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.
fn miller_rabin(rng: *io::stream, x: []word, n: int, t: []word) u32 = {
	// 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).

	const tlen = len(t);
	// Compute (x-1)/2 (encoded).
	const xm1d2_len = ((x[0] - (x[0] >> 5)) + 7) >> 3;
	let xm1d2 = (t: *[*]u8)[..xm1d2_len];
	decode(xm1d2, x);

	let cc: uint = 0;
	for (let u = 0z; u < xm1d2_len; u += 1) {
		let w = xm1d2[u];
		xm1d2[u] = (w >> 1) | cc: u8;
		cc = w << 7;
	};

	// We used some words of the provided buffer for (x-1)/2.
	const xm1d2_len_u32 = (xm1d2_len + 3) >> 2;
	let t = t[xm1d2_len_u32..];

	const xlen = (x[0] + 31) >> 5;
	const asize: word = x[0] - 1 - eq0u32(x[0] & 31);
	const x0i = ninv(x[1]);
	for (n > 0; n -= 1) {
		//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.
		let 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).
		let t2 = t[1 + xlen..];
		if ((len(t2) & 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 = t2[1..];
		};

		modpow(a, xm1d2, x, x0i, t2);

		// 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).
		const eq1 = a[1] ^ 1;
		const eqm1 = a[1] ^ (x[1] - 1);
		for (let u: word = 2; u <= xlen; u += 1) {
			eq1 |= a[u];
			eqm1 |= a[u] ^ x[u];
		};

		if ((eq0u32(eq1) | eq0u32(eqm1)) == 0) {
			return 0;
		};
	};
	return 1;
};

fn dump(x: []word) void = {
	fmt::printf("[{}]", x[0]: u32)!;// , x[0]: u32, hex::encodestr((x[1..]: *[*]u8)[..(len(x)-1) * size(word)]))!;
	let dest: [4096]u8 = [0...];
	decode(dest, x);
	fmt::printf(strings::ltrim(hex::encodestr(dest), '0'))!;
	fmt::println()!;
};

// 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.
// 'esize' is the encoded size of resulting 'x'.
export fn mkprime(rng: *io::stream, x: []word, esize: word, pubexp: u32, t: []word) void =
	//uint32_t *x, uint32_t esize,
	//uint32_t pubexp, uint32_t *t, size_t tlen, br_i31_modpow_opt_type mp31)
{
	x[0] = esize;
	const sz = (esize + 31) >> 5;
	for (true) {
		// 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);
		dump(x);
		if ((esize & 31) == 0) {
			x[sz] |= 0x60000000;
		} else if ((esize & 31) == 1) {
			x[sz] |= 0x00000001;
			x[sz - 1] |= 0x40000000;
		} else {
			x[sz] |= 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
		let m3: word = 0;
		let m5: word = 0;
		let m7: word = 0;
		let m11: word = 0;
		let s7: int = 0;
		let s11: int = 0;
		for (let u: word = 0; u < sz; u += 1) {

			let w: word = x[1 + u];
			let w3: word = (w & 0xFFFF) + (w >> 16);   // max: 98302
			let w5: word = (w & 0xFFFF) + (w >> 16);   // max: 98302
			let w7: word = (w & 0x7FFF) + (w >> 15);   // max: 98302
			let w11: word = (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: word;
			m7 = (m7 & 0x1FF) + (m7 >> 9);     // max: 1280
			s7 += 1;
			if (s7 == 3) {
				s7 = 0;
			};

			m11 += w11 << s11: word;
			s11 += 1;
			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: 2
		m5 = (m5 & 0x0F) + (m5 >> 4);      // max: 31
		m5 -= 20 & -gtu32(m5, 19);
		m5 -= 10 & -gtu32(m5, 9);
		m5 -= 5 & -gtu32(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 & -gtu32(m11, 43);
		m11 -= 22 & -gtu32(m11, 21);
		m11 -= 11 & -gtu32(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) {
			fmt::println("  mod failed")!;
			continue;
		};
		if ((pubexp == 3 && m3 == 1)
				|| (pubexp == 5 && m5 == 1)
				|| (pubexp == 7 && m7 == 1)
				|| (pubexp == 11 && m11 == 1)) {
			fmt::println("  pubexp failed")!;
			continue;
		};

		// More trial divisions.
		if (trial_divisions(x, t) == 0) {
			fmt::println("  trial failed")!;
			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.
		const rounds: int = if (esize < 309) {
			yield 12;
		} else if (esize < 464) {
			yield 9;
		} else if (esize < 670) {
			yield 6;
		} else if (esize < 877) {
			yield 4;
		} else if (esize < 1341) {
			yield 3;
		} else {
			yield 2;
		};

		if (miller_rabin(rng, x, rounds, t) == 1) {
			return;
		};
		fmt::println("  miller failed")!;
	};
};