xref: /freebsd/contrib/bearssl/test/test_math.c (revision 0fca6ea1d4eea4c934cfff25ac9ee8ad6fe95583)
1 /*
2  * Copyright (c) 2016 Thomas Pornin <pornin@bolet.org>
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining
5  * a copy of this software and associated documentation files (the
6  * "Software"), to deal in the Software without restriction, including
7  * without limitation the rights to use, copy, modify, merge, publish,
8  * distribute, sublicense, and/or sell copies of the Software, and to
9  * permit persons to whom the Software is furnished to do so, subject to
10  * the following conditions:
11  *
12  * The above copyright notice and this permission notice shall be
13  * included in all copies or substantial portions of the Software.
14  *
15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
16  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
18  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
19  * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
20  * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
21  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
24 
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <string.h>
28 #include <stdarg.h>
29 #include <time.h>
30 
31 #include <gmp.h>
32 
33 #include "bearssl.h"
34 #include "inner.h"
35 
36 /*
37  * Pointers to implementations.
38  */
39 typedef struct {
40 	uint32_t word_size;
41 	void (*zero)(uint32_t *x, uint32_t bit_len);
42 	void (*decode)(uint32_t *x, const void *src, size_t len);
43 	uint32_t (*decode_mod)(uint32_t *x,
44 		const void *src, size_t len, const uint32_t *m);
45 	void (*reduce)(uint32_t *x, const uint32_t *a, const uint32_t *m);
46 	void (*decode_reduce)(uint32_t *x,
47 		const void *src, size_t len, const uint32_t *m);
48 	void (*encode)(void *dst, size_t len, const uint32_t *x);
49 	uint32_t (*add)(uint32_t *a, const uint32_t *b, uint32_t ctl);
50 	uint32_t (*sub)(uint32_t *a, const uint32_t *b, uint32_t ctl);
51 	uint32_t (*ninv)(uint32_t x);
52 	void (*montymul)(uint32_t *d, const uint32_t *x, const uint32_t *y,
53 		const uint32_t *m, uint32_t m0i);
54 	void (*to_monty)(uint32_t *x, const uint32_t *m);
55 	void (*from_monty)(uint32_t *x, const uint32_t *m, uint32_t m0i);
56 	void (*modpow)(uint32_t *x, const unsigned char *e, size_t elen,
57 		const uint32_t *m, uint32_t m0i, uint32_t *t1, uint32_t *t2);
58 } int_impl;
59 
60 static const int_impl i31_impl = {
61 	31,
62 	&br_i31_zero,
63 	&br_i31_decode,
64 	&br_i31_decode_mod,
65 	&br_i31_reduce,
66 	&br_i31_decode_reduce,
67 	&br_i31_encode,
68 	&br_i31_add,
69 	&br_i31_sub,
70 	&br_i31_ninv31,
71 	&br_i31_montymul,
72 	&br_i31_to_monty,
73 	&br_i31_from_monty,
74 	&br_i31_modpow
75 };
76 static const int_impl i32_impl = {
77 	32,
78 	&br_i32_zero,
79 	&br_i32_decode,
80 	&br_i32_decode_mod,
81 	&br_i32_reduce,
82 	&br_i32_decode_reduce,
83 	&br_i32_encode,
84 	&br_i32_add,
85 	&br_i32_sub,
86 	&br_i32_ninv32,
87 	&br_i32_montymul,
88 	&br_i32_to_monty,
89 	&br_i32_from_monty,
90 	&br_i32_modpow
91 };
92 
93 static const int_impl *impl;
94 
95 static gmp_randstate_t RNG;
96 
97 /*
98  * Get a random prime of length 'size' bits. This function also guarantees
99  * that x-1 is not a multiple of 65537.
100  */
101 static void
102 rand_prime(mpz_t x, int size)
103 {
104 	for (;;) {
105 		mpz_urandomb(x, RNG, size - 1);
106 		mpz_setbit(x, 0);
107 		mpz_setbit(x, size - 1);
108 		if (mpz_probab_prime_p(x, 50)) {
109 			mpz_sub_ui(x, x, 1);
110 			if (mpz_divisible_ui_p(x, 65537)) {
111 				continue;
112 			}
113 			mpz_add_ui(x, x, 1);
114 			return;
115 		}
116 	}
117 }
118 
119 /*
120  * Print out a GMP integer (for debug).
121  */
122 static void
123 print_z(mpz_t z)
124 {
125 	unsigned char zb[1000];
126 	size_t zlen, k;
127 
128 	mpz_export(zb, &zlen, 1, 1, 0, 0, z);
129 	if (zlen == 0) {
130 		printf(" 00");
131 		return;
132 	}
133 	if ((zlen & 3) != 0) {
134 		k = 4 - (zlen & 3);
135 		memmove(zb + k, zb, zlen);
136 		memset(zb, 0, k);
137 		zlen += k;
138 	}
139 	for (k = 0; k < zlen; k += 4) {
140 		printf(" %02X%02X%02X%02X",
141 			zb[k], zb[k + 1], zb[k + 2], zb[k + 3]);
142 	}
143 }
144 
145 /*
146  * Print out an i31 or i32 integer (for debug).
147  */
148 static void
149 print_u(uint32_t *x)
150 {
151 	size_t k;
152 
153 	if (x[0] == 0) {
154 		printf(" 00000000 (0, 0)");
155 		return;
156 	}
157 	for (k = (x[0] + 31) >> 5; k > 0; k --) {
158 		printf(" %08lX", (unsigned long)x[k]);
159 	}
160 	printf(" (%u, %u)", (unsigned)(x[0] >> 5), (unsigned)(x[0] & 31));
161 }
162 
163 /*
164  * Check that an i31/i32 number and a GMP number are equal.
165  */
166 static void
167 check_eqz(uint32_t *x, mpz_t z)
168 {
169 	unsigned char xb[1000];
170 	unsigned char zb[1000];
171 	size_t xlen, zlen;
172 	int good;
173 
174 	xlen = ((x[0] + 31) & ~(uint32_t)31) >> 3;
175 	impl->encode(xb, xlen, x);
176 	mpz_export(zb, &zlen, 1, 1, 0, 0, z);
177 	good = 1;
178 	if (xlen < zlen) {
179 		good = 0;
180 	} else if (xlen > zlen) {
181 		size_t u;
182 
183 		for (u = xlen; u > zlen; u --) {
184 			if (xb[xlen - u] != 0) {
185 				good = 0;
186 				break;
187 			}
188 		}
189 	}
190 	good = good && memcmp(xb + xlen - zlen, zb, zlen) == 0;
191 	if (!good) {
192 		size_t u;
193 
194 		printf("Mismatch:\n");
195 		printf("  x = ");
196 		print_u(x);
197 		printf("\n");
198 		printf("  ex = ");
199 		for (u = 0; u < xlen; u ++) {
200 			printf("%02X", xb[u]);
201 		}
202 		printf("\n");
203 		printf("  z = ");
204 		print_z(z);
205 		printf("\n");
206 		exit(EXIT_FAILURE);
207 	}
208 }
209 
210 /* obsolete
211 static void
212 mp_to_br(uint32_t *mx, uint32_t x_bitlen, mpz_t x)
213 {
214 	uint32_t x_ebitlen;
215 	size_t xlen;
216 
217 	if (mpz_sizeinbase(x, 2) > x_bitlen) {
218 		abort();
219 	}
220 	x_ebitlen = ((x_bitlen / 31) << 5) + (x_bitlen % 31);
221 	br_i31_zero(mx, x_ebitlen);
222 	mpz_export(mx + 1, &xlen, -1, sizeof *mx, 0, 1, x);
223 }
224 */
225 
226 static void
227 test_modint(void)
228 {
229 	int i, j, k;
230 	mpz_t p, a, b, v, t1;
231 
232 	printf("Test modular integers: ");
233 	fflush(stdout);
234 
235 	gmp_randinit_mt(RNG);
236 	mpz_init(p);
237 	mpz_init(a);
238 	mpz_init(b);
239 	mpz_init(v);
240 	mpz_init(t1);
241 	mpz_set_ui(t1, (unsigned long)time(NULL));
242 	gmp_randseed(RNG, t1);
243 	for (k = 2; k <= 128; k ++) {
244 		for (i = 0; i < 10; i ++) {
245 			unsigned char ep[100], ea[100], eb[100], ev[100];
246 			size_t plen, alen, blen, vlen;
247 			uint32_t mp[40], ma[40], mb[40], mv[60], mx[100];
248 			uint32_t mt1[40], mt2[40], mt3[40];
249 			uint32_t ctl;
250 			uint32_t mp0i;
251 
252 			rand_prime(p, k);
253 			mpz_urandomm(a, RNG, p);
254 			mpz_urandomm(b, RNG, p);
255 			mpz_urandomb(v, RNG, k + 60);
256 			if (mpz_sgn(b) == 0) {
257 				mpz_set_ui(b, 1);
258 			}
259 			mpz_export(ep, &plen, 1, 1, 0, 0, p);
260 			mpz_export(ea, &alen, 1, 1, 0, 0, a);
261 			mpz_export(eb, &blen, 1, 1, 0, 0, b);
262 			mpz_export(ev, &vlen, 1, 1, 0, 0, v);
263 
264 			impl->decode(mp, ep, plen);
265 			if (impl->decode_mod(ma, ea, alen, mp) != 1) {
266 				printf("Decode error\n");
267 				printf("  ea = ");
268 				print_z(a);
269 				printf("\n");
270 				printf("  p = ");
271 				print_u(mp);
272 				printf("\n");
273 				exit(EXIT_FAILURE);
274 			}
275 			mp0i = impl->ninv(mp[1]);
276 			if (impl->decode_mod(mb, eb, blen, mp) != 1) {
277 				printf("Decode error\n");
278 				printf("  eb = ");
279 				print_z(b);
280 				printf("\n");
281 				printf("  p = ");
282 				print_u(mp);
283 				printf("\n");
284 				exit(EXIT_FAILURE);
285 			}
286 			impl->decode(mv, ev, vlen);
287 			check_eqz(mp, p);
288 			check_eqz(ma, a);
289 			check_eqz(mb, b);
290 			check_eqz(mv, v);
291 
292 			impl->decode_mod(ma, ea, alen, mp);
293 			impl->decode_mod(mb, eb, blen, mp);
294 			ctl = impl->add(ma, mb, 1);
295 			ctl |= impl->sub(ma, mp, 0) ^ (uint32_t)1;
296 			impl->sub(ma, mp, ctl);
297 			mpz_add(t1, a, b);
298 			mpz_mod(t1, t1, p);
299 			check_eqz(ma, t1);
300 
301 			impl->decode_mod(ma, ea, alen, mp);
302 			impl->decode_mod(mb, eb, blen, mp);
303 			impl->add(ma, mp, impl->sub(ma, mb, 1));
304 			mpz_sub(t1, a, b);
305 			mpz_mod(t1, t1, p);
306 			check_eqz(ma, t1);
307 
308 			impl->decode_reduce(ma, ev, vlen, mp);
309 			mpz_mod(t1, v, p);
310 			check_eqz(ma, t1);
311 
312 			impl->decode(mv, ev, vlen);
313 			impl->reduce(ma, mv, mp);
314 			mpz_mod(t1, v, p);
315 			check_eqz(ma, t1);
316 
317 			impl->decode_mod(ma, ea, alen, mp);
318 			impl->to_monty(ma, mp);
319 			mpz_mul_2exp(t1, a, ((k + impl->word_size - 1)
320 				/ impl->word_size) * impl->word_size);
321 			mpz_mod(t1, t1, p);
322 			check_eqz(ma, t1);
323 			impl->from_monty(ma, mp, mp0i);
324 			check_eqz(ma, a);
325 
326 			impl->decode_mod(ma, ea, alen, mp);
327 			impl->decode_mod(mb, eb, blen, mp);
328 			impl->to_monty(ma, mp);
329 			impl->montymul(mt1, ma, mb, mp, mp0i);
330 			mpz_mul(t1, a, b);
331 			mpz_mod(t1, t1, p);
332 			check_eqz(mt1, t1);
333 
334 			impl->decode_mod(ma, ea, alen, mp);
335 			impl->modpow(ma, ev, vlen, mp, mp0i, mt1, mt2);
336 			mpz_powm(t1, a, v, p);
337 			check_eqz(ma, t1);
338 
339 			/*
340 			br_modint_decode(ma, mp, ea, alen);
341 			br_modint_decode(mb, mp, eb, blen);
342 			if (!br_modint_div(ma, mb, mp, mt1, mt2, mt3)) {
343 				fprintf(stderr, "division failed\n");
344 				exit(EXIT_FAILURE);
345 			}
346 			mpz_sub_ui(t1, p, 2);
347 			mpz_powm(t1, b, t1, p);
348 			mpz_mul(t1, a, t1);
349 			mpz_mod(t1, t1, p);
350 			check_eqz(ma, t1);
351 
352 			br_modint_decode(ma, mp, ea, alen);
353 			br_modint_decode(mb, mp, eb, blen);
354 			for (j = 0; j <= (2 * k + 5); j ++) {
355 				br_int_add(mx, j, ma, mb);
356 				mpz_add(t1, a, b);
357 				mpz_tdiv_r_2exp(t1, t1, j);
358 				check_eqz(mx, t1);
359 
360 				br_int_mul(mx, j, ma, mb);
361 				mpz_mul(t1, a, b);
362 				mpz_tdiv_r_2exp(t1, t1, j);
363 				check_eqz(mx, t1);
364 			}
365 			*/
366 		}
367 		printf(".");
368 		fflush(stdout);
369 	}
370 	mpz_clear(p);
371 	mpz_clear(a);
372 	mpz_clear(b);
373 	mpz_clear(v);
374 	mpz_clear(t1);
375 
376 	printf(" done.\n");
377 	fflush(stdout);
378 }
379 
380 #if 0
381 static void
382 test_RSA_core(void)
383 {
384 	int i, j, k;
385 	mpz_t n, e, d, p, q, dp, dq, iq, t1, t2, phi;
386 
387 	printf("Test RSA core: ");
388 	fflush(stdout);
389 
390 	gmp_randinit_mt(RNG);
391 	mpz_init(n);
392 	mpz_init(e);
393 	mpz_init(d);
394 	mpz_init(p);
395 	mpz_init(q);
396 	mpz_init(dp);
397 	mpz_init(dq);
398 	mpz_init(iq);
399 	mpz_init(t1);
400 	mpz_init(t2);
401 	mpz_init(phi);
402 	mpz_set_ui(t1, (unsigned long)time(NULL));
403 	gmp_randseed(RNG, t1);
404 
405 	/*
406 	 * To test corner cases, we want to try RSA keys such that the
407 	 * lengths of both factors can be arbitrary modulo 2^32. Factors
408 	 * p and q need not be of the same length; p can be greater than
409 	 * q and q can be greater than p.
410 	 *
411 	 * To keep computation time reasonable, we use p and q factors of
412 	 * less than 128 bits; this is way too small for secure RSA,
413 	 * but enough to exercise all code paths (since we work only with
414 	 * 32-bit words).
415 	 */
416 	for (i = 64; i <= 96; i ++) {
417 		rand_prime(p, i);
418 		for (j = i - 33; j <= i + 33; j ++) {
419 			uint32_t mp[40], mq[40], mdp[40], mdq[40], miq[40];
420 
421 			/*
422 			 * Generate a RSA key pair, with p of length i bits,
423 			 * and q of length j bits.
424 			 */
425 			do {
426 				rand_prime(q, j);
427 			} while (mpz_cmp(p, q) == 0);
428 			mpz_mul(n, p, q);
429 			mpz_set_ui(e, 65537);
430 			mpz_sub_ui(t1, p, 1);
431 			mpz_sub_ui(t2, q, 1);
432 			mpz_mul(phi, t1, t2);
433 			mpz_invert(d, e, phi);
434 			mpz_mod(dp, d, t1);
435 			mpz_mod(dq, d, t2);
436 			mpz_invert(iq, q, p);
437 
438 			/*
439 			 * Convert the key pair elements to BearSSL arrays.
440 			 */
441 			mp_to_br(mp, mpz_sizeinbase(p, 2), p);
442 			mp_to_br(mq, mpz_sizeinbase(q, 2), q);
443 			mp_to_br(mdp, mpz_sizeinbase(dp, 2), dp);
444 			mp_to_br(mdq, mpz_sizeinbase(dq, 2), dq);
445 			mp_to_br(miq, mp[0], iq);
446 
447 			/*
448 			 * Compute and check ten public/private operations.
449 			 */
450 			for (k = 0; k < 10; k ++) {
451 				uint32_t mx[40];
452 
453 				mpz_urandomm(t1, RNG, n);
454 				mpz_powm(t2, t1, e, n);
455 				mp_to_br(mx, mpz_sizeinbase(n, 2), t2);
456 				br_rsa_private_core(mx, mp, mq, mdp, mdq, miq);
457 				check_eqz(mx, t1);
458 			}
459 		}
460 		printf(".");
461 		fflush(stdout);
462 	}
463 
464 	printf(" done.\n");
465 	fflush(stdout);
466 }
467 #endif
468 
469 int
470 main(void)
471 {
472 	printf("===== i32 ======\n");
473 	impl = &i32_impl;
474 	test_modint();
475 	printf("===== i31 ======\n");
476 	impl = &i31_impl;
477 	test_modint();
478 	/*
479 	test_RSA_core();
480 	*/
481 	return 0;
482 }
483