1 /*
2 * Copyright 2002-2021 The OpenSSL Project Authors. All Rights Reserved.
3 * Copyright (c) 2002, Oracle and/or its affiliates. All rights reserved
4 *
5 * Licensed under the Apache License 2.0 (the "License"). You may not use
6 * this file except in compliance with the License. You can obtain a copy
7 * in the file LICENSE in the source distribution or at
8 * https://www.openssl.org/source/license.html
9 */
10
11 #include <assert.h>
12 #include <limits.h>
13 #include <stdio.h>
14 #include "internal/cryptlib.h"
15 #include "bn_local.h"
16
17 #ifndef OPENSSL_NO_EC2M
18
19 /*
20 * Maximum number of iterations before BN_GF2m_mod_solve_quad_arr should
21 * fail.
22 */
23 # define MAX_ITERATIONS 50
24
25 # define SQR_nibble(w) ((((w) & 8) << 3) \
26 | (((w) & 4) << 2) \
27 | (((w) & 2) << 1) \
28 | ((w) & 1))
29
30
31 /* Platform-specific macros to accelerate squaring. */
32 # if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
33 # define SQR1(w) \
34 SQR_nibble((w) >> 60) << 56 | SQR_nibble((w) >> 56) << 48 | \
35 SQR_nibble((w) >> 52) << 40 | SQR_nibble((w) >> 48) << 32 | \
36 SQR_nibble((w) >> 44) << 24 | SQR_nibble((w) >> 40) << 16 | \
37 SQR_nibble((w) >> 36) << 8 | SQR_nibble((w) >> 32)
38 # define SQR0(w) \
39 SQR_nibble((w) >> 28) << 56 | SQR_nibble((w) >> 24) << 48 | \
40 SQR_nibble((w) >> 20) << 40 | SQR_nibble((w) >> 16) << 32 | \
41 SQR_nibble((w) >> 12) << 24 | SQR_nibble((w) >> 8) << 16 | \
42 SQR_nibble((w) >> 4) << 8 | SQR_nibble((w) )
43 # endif
44 # ifdef THIRTY_TWO_BIT
45 # define SQR1(w) \
46 SQR_nibble((w) >> 28) << 24 | SQR_nibble((w) >> 24) << 16 | \
47 SQR_nibble((w) >> 20) << 8 | SQR_nibble((w) >> 16)
48 # define SQR0(w) \
49 SQR_nibble((w) >> 12) << 24 | SQR_nibble((w) >> 8) << 16 | \
50 SQR_nibble((w) >> 4) << 8 | SQR_nibble((w) )
51 # endif
52
53 # if !defined(OPENSSL_BN_ASM_GF2m)
54 /*
55 * Product of two polynomials a, b each with degree < BN_BITS2 - 1, result is
56 * a polynomial r with degree < 2 * BN_BITS - 1 The caller MUST ensure that
57 * the variables have the right amount of space allocated.
58 */
59 # ifdef THIRTY_TWO_BIT
bn_GF2m_mul_1x1(BN_ULONG * r1,BN_ULONG * r0,const BN_ULONG a,const BN_ULONG b)60 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
61 const BN_ULONG b)
62 {
63 register BN_ULONG h, l, s;
64 BN_ULONG tab[8], top2b = a >> 30;
65 register BN_ULONG a1, a2, a4;
66
67 a1 = a & (0x3FFFFFFF);
68 a2 = a1 << 1;
69 a4 = a2 << 1;
70
71 tab[0] = 0;
72 tab[1] = a1;
73 tab[2] = a2;
74 tab[3] = a1 ^ a2;
75 tab[4] = a4;
76 tab[5] = a1 ^ a4;
77 tab[6] = a2 ^ a4;
78 tab[7] = a1 ^ a2 ^ a4;
79
80 s = tab[b & 0x7];
81 l = s;
82 s = tab[b >> 3 & 0x7];
83 l ^= s << 3;
84 h = s >> 29;
85 s = tab[b >> 6 & 0x7];
86 l ^= s << 6;
87 h ^= s >> 26;
88 s = tab[b >> 9 & 0x7];
89 l ^= s << 9;
90 h ^= s >> 23;
91 s = tab[b >> 12 & 0x7];
92 l ^= s << 12;
93 h ^= s >> 20;
94 s = tab[b >> 15 & 0x7];
95 l ^= s << 15;
96 h ^= s >> 17;
97 s = tab[b >> 18 & 0x7];
98 l ^= s << 18;
99 h ^= s >> 14;
100 s = tab[b >> 21 & 0x7];
101 l ^= s << 21;
102 h ^= s >> 11;
103 s = tab[b >> 24 & 0x7];
104 l ^= s << 24;
105 h ^= s >> 8;
106 s = tab[b >> 27 & 0x7];
107 l ^= s << 27;
108 h ^= s >> 5;
109 s = tab[b >> 30];
110 l ^= s << 30;
111 h ^= s >> 2;
112
113 /* compensate for the top two bits of a */
114
115 if (top2b & 01) {
116 l ^= b << 30;
117 h ^= b >> 2;
118 }
119 if (top2b & 02) {
120 l ^= b << 31;
121 h ^= b >> 1;
122 }
123
124 *r1 = h;
125 *r0 = l;
126 }
127 # endif
128 # if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
bn_GF2m_mul_1x1(BN_ULONG * r1,BN_ULONG * r0,const BN_ULONG a,const BN_ULONG b)129 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
130 const BN_ULONG b)
131 {
132 register BN_ULONG h, l, s;
133 BN_ULONG tab[16], top3b = a >> 61;
134 register BN_ULONG a1, a2, a4, a8;
135
136 a1 = a & (0x1FFFFFFFFFFFFFFFULL);
137 a2 = a1 << 1;
138 a4 = a2 << 1;
139 a8 = a4 << 1;
140
141 tab[0] = 0;
142 tab[1] = a1;
143 tab[2] = a2;
144 tab[3] = a1 ^ a2;
145 tab[4] = a4;
146 tab[5] = a1 ^ a4;
147 tab[6] = a2 ^ a4;
148 tab[7] = a1 ^ a2 ^ a4;
149 tab[8] = a8;
150 tab[9] = a1 ^ a8;
151 tab[10] = a2 ^ a8;
152 tab[11] = a1 ^ a2 ^ a8;
153 tab[12] = a4 ^ a8;
154 tab[13] = a1 ^ a4 ^ a8;
155 tab[14] = a2 ^ a4 ^ a8;
156 tab[15] = a1 ^ a2 ^ a4 ^ a8;
157
158 s = tab[b & 0xF];
159 l = s;
160 s = tab[b >> 4 & 0xF];
161 l ^= s << 4;
162 h = s >> 60;
163 s = tab[b >> 8 & 0xF];
164 l ^= s << 8;
165 h ^= s >> 56;
166 s = tab[b >> 12 & 0xF];
167 l ^= s << 12;
168 h ^= s >> 52;
169 s = tab[b >> 16 & 0xF];
170 l ^= s << 16;
171 h ^= s >> 48;
172 s = tab[b >> 20 & 0xF];
173 l ^= s << 20;
174 h ^= s >> 44;
175 s = tab[b >> 24 & 0xF];
176 l ^= s << 24;
177 h ^= s >> 40;
178 s = tab[b >> 28 & 0xF];
179 l ^= s << 28;
180 h ^= s >> 36;
181 s = tab[b >> 32 & 0xF];
182 l ^= s << 32;
183 h ^= s >> 32;
184 s = tab[b >> 36 & 0xF];
185 l ^= s << 36;
186 h ^= s >> 28;
187 s = tab[b >> 40 & 0xF];
188 l ^= s << 40;
189 h ^= s >> 24;
190 s = tab[b >> 44 & 0xF];
191 l ^= s << 44;
192 h ^= s >> 20;
193 s = tab[b >> 48 & 0xF];
194 l ^= s << 48;
195 h ^= s >> 16;
196 s = tab[b >> 52 & 0xF];
197 l ^= s << 52;
198 h ^= s >> 12;
199 s = tab[b >> 56 & 0xF];
200 l ^= s << 56;
201 h ^= s >> 8;
202 s = tab[b >> 60];
203 l ^= s << 60;
204 h ^= s >> 4;
205
206 /* compensate for the top three bits of a */
207
208 if (top3b & 01) {
209 l ^= b << 61;
210 h ^= b >> 3;
211 }
212 if (top3b & 02) {
213 l ^= b << 62;
214 h ^= b >> 2;
215 }
216 if (top3b & 04) {
217 l ^= b << 63;
218 h ^= b >> 1;
219 }
220
221 *r1 = h;
222 *r0 = l;
223 }
224 # endif
225
226 /*
227 * Product of two polynomials a, b each with degree < 2 * BN_BITS2 - 1,
228 * result is a polynomial r with degree < 4 * BN_BITS2 - 1 The caller MUST
229 * ensure that the variables have the right amount of space allocated.
230 */
bn_GF2m_mul_2x2(BN_ULONG * r,const BN_ULONG a1,const BN_ULONG a0,const BN_ULONG b1,const BN_ULONG b0)231 static void bn_GF2m_mul_2x2(BN_ULONG *r, const BN_ULONG a1, const BN_ULONG a0,
232 const BN_ULONG b1, const BN_ULONG b0)
233 {
234 BN_ULONG m1, m0;
235 /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */
236 bn_GF2m_mul_1x1(r + 3, r + 2, a1, b1);
237 bn_GF2m_mul_1x1(r + 1, r, a0, b0);
238 bn_GF2m_mul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1);
239 /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */
240 r[2] ^= m1 ^ r[1] ^ r[3]; /* h0 ^= m1 ^ l1 ^ h1; */
241 r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */
242 }
243 # else
244 void bn_GF2m_mul_2x2(BN_ULONG *r, BN_ULONG a1, BN_ULONG a0, BN_ULONG b1,
245 BN_ULONG b0);
246 # endif
247
248 /*
249 * Add polynomials a and b and store result in r; r could be a or b, a and b
250 * could be equal; r is the bitwise XOR of a and b.
251 */
BN_GF2m_add(BIGNUM * r,const BIGNUM * a,const BIGNUM * b)252 int BN_GF2m_add(BIGNUM *r, const BIGNUM *a, const BIGNUM *b)
253 {
254 int i;
255 const BIGNUM *at, *bt;
256
257 bn_check_top(a);
258 bn_check_top(b);
259
260 if (a->top < b->top) {
261 at = b;
262 bt = a;
263 } else {
264 at = a;
265 bt = b;
266 }
267
268 if (bn_wexpand(r, at->top) == NULL)
269 return 0;
270
271 for (i = 0; i < bt->top; i++) {
272 r->d[i] = at->d[i] ^ bt->d[i];
273 }
274 for (; i < at->top; i++) {
275 r->d[i] = at->d[i];
276 }
277
278 r->top = at->top;
279 bn_correct_top(r);
280
281 return 1;
282 }
283
284 /*-
285 * Some functions allow for representation of the irreducible polynomials
286 * as an int[], say p. The irreducible f(t) is then of the form:
287 * t^p[0] + t^p[1] + ... + t^p[k]
288 * where m = p[0] > p[1] > ... > p[k] = 0.
289 */
290
291 /* Performs modular reduction of a and store result in r. r could be a. */
BN_GF2m_mod_arr(BIGNUM * r,const BIGNUM * a,const int p[])292 int BN_GF2m_mod_arr(BIGNUM *r, const BIGNUM *a, const int p[])
293 {
294 int j, k;
295 int n, dN, d0, d1;
296 BN_ULONG zz, *z;
297
298 bn_check_top(a);
299
300 if (p[0] == 0) {
301 /* reduction mod 1 => return 0 */
302 BN_zero(r);
303 return 1;
304 }
305
306 /*
307 * Since the algorithm does reduction in the r value, if a != r, copy the
308 * contents of a into r so we can do reduction in r.
309 */
310 if (a != r) {
311 if (!bn_wexpand(r, a->top))
312 return 0;
313 for (j = 0; j < a->top; j++) {
314 r->d[j] = a->d[j];
315 }
316 r->top = a->top;
317 }
318 z = r->d;
319
320 /* start reduction */
321 dN = p[0] / BN_BITS2;
322 for (j = r->top - 1; j > dN;) {
323 zz = z[j];
324 if (z[j] == 0) {
325 j--;
326 continue;
327 }
328 z[j] = 0;
329
330 for (k = 1; p[k] != 0; k++) {
331 /* reducing component t^p[k] */
332 n = p[0] - p[k];
333 d0 = n % BN_BITS2;
334 d1 = BN_BITS2 - d0;
335 n /= BN_BITS2;
336 z[j - n] ^= (zz >> d0);
337 if (d0)
338 z[j - n - 1] ^= (zz << d1);
339 }
340
341 /* reducing component t^0 */
342 n = dN;
343 d0 = p[0] % BN_BITS2;
344 d1 = BN_BITS2 - d0;
345 z[j - n] ^= (zz >> d0);
346 if (d0)
347 z[j - n - 1] ^= (zz << d1);
348 }
349
350 /* final round of reduction */
351 while (j == dN) {
352
353 d0 = p[0] % BN_BITS2;
354 zz = z[dN] >> d0;
355 if (zz == 0)
356 break;
357 d1 = BN_BITS2 - d0;
358
359 /* clear up the top d1 bits */
360 if (d0)
361 z[dN] = (z[dN] << d1) >> d1;
362 else
363 z[dN] = 0;
364 z[0] ^= zz; /* reduction t^0 component */
365
366 for (k = 1; p[k] != 0; k++) {
367 BN_ULONG tmp_ulong;
368
369 /* reducing component t^p[k] */
370 n = p[k] / BN_BITS2;
371 d0 = p[k] % BN_BITS2;
372 d1 = BN_BITS2 - d0;
373 z[n] ^= (zz << d0);
374 if (d0 && (tmp_ulong = zz >> d1))
375 z[n + 1] ^= tmp_ulong;
376 }
377
378 }
379
380 bn_correct_top(r);
381 return 1;
382 }
383
384 /*
385 * Performs modular reduction of a by p and store result in r. r could be a.
386 * This function calls down to the BN_GF2m_mod_arr implementation; this wrapper
387 * function is only provided for convenience; for best performance, use the
388 * BN_GF2m_mod_arr function.
389 */
BN_GF2m_mod(BIGNUM * r,const BIGNUM * a,const BIGNUM * p)390 int BN_GF2m_mod(BIGNUM *r, const BIGNUM *a, const BIGNUM *p)
391 {
392 int ret = 0;
393 int arr[6];
394 bn_check_top(a);
395 bn_check_top(p);
396 ret = BN_GF2m_poly2arr(p, arr, OSSL_NELEM(arr));
397 if (!ret || ret > (int)OSSL_NELEM(arr)) {
398 ERR_raise(ERR_LIB_BN, BN_R_INVALID_LENGTH);
399 return 0;
400 }
401 ret = BN_GF2m_mod_arr(r, a, arr);
402 bn_check_top(r);
403 return ret;
404 }
405
406 /*
407 * Compute the product of two polynomials a and b, reduce modulo p, and store
408 * the result in r. r could be a or b; a could be b.
409 */
BN_GF2m_mod_mul_arr(BIGNUM * r,const BIGNUM * a,const BIGNUM * b,const int p[],BN_CTX * ctx)410 int BN_GF2m_mod_mul_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
411 const int p[], BN_CTX *ctx)
412 {
413 int zlen, i, j, k, ret = 0;
414 BIGNUM *s;
415 BN_ULONG x1, x0, y1, y0, zz[4];
416
417 bn_check_top(a);
418 bn_check_top(b);
419
420 if (a == b) {
421 return BN_GF2m_mod_sqr_arr(r, a, p, ctx);
422 }
423
424 BN_CTX_start(ctx);
425 if ((s = BN_CTX_get(ctx)) == NULL)
426 goto err;
427
428 zlen = a->top + b->top + 4;
429 if (!bn_wexpand(s, zlen))
430 goto err;
431 s->top = zlen;
432
433 for (i = 0; i < zlen; i++)
434 s->d[i] = 0;
435
436 for (j = 0; j < b->top; j += 2) {
437 y0 = b->d[j];
438 y1 = ((j + 1) == b->top) ? 0 : b->d[j + 1];
439 for (i = 0; i < a->top; i += 2) {
440 x0 = a->d[i];
441 x1 = ((i + 1) == a->top) ? 0 : a->d[i + 1];
442 bn_GF2m_mul_2x2(zz, x1, x0, y1, y0);
443 for (k = 0; k < 4; k++)
444 s->d[i + j + k] ^= zz[k];
445 }
446 }
447
448 bn_correct_top(s);
449 if (BN_GF2m_mod_arr(r, s, p))
450 ret = 1;
451 bn_check_top(r);
452
453 err:
454 BN_CTX_end(ctx);
455 return ret;
456 }
457
458 /*
459 * Compute the product of two polynomials a and b, reduce modulo p, and store
460 * the result in r. r could be a or b; a could equal b. This function calls
461 * down to the BN_GF2m_mod_mul_arr implementation; this wrapper function is
462 * only provided for convenience; for best performance, use the
463 * BN_GF2m_mod_mul_arr function.
464 */
BN_GF2m_mod_mul(BIGNUM * r,const BIGNUM * a,const BIGNUM * b,const BIGNUM * p,BN_CTX * ctx)465 int BN_GF2m_mod_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
466 const BIGNUM *p, BN_CTX *ctx)
467 {
468 int ret = 0;
469 const int max = BN_num_bits(p) + 1;
470 int *arr;
471
472 bn_check_top(a);
473 bn_check_top(b);
474 bn_check_top(p);
475
476 arr = OPENSSL_malloc(sizeof(*arr) * max);
477 if (arr == NULL) {
478 ERR_raise(ERR_LIB_BN, ERR_R_MALLOC_FAILURE);
479 return 0;
480 }
481 ret = BN_GF2m_poly2arr(p, arr, max);
482 if (!ret || ret > max) {
483 ERR_raise(ERR_LIB_BN, BN_R_INVALID_LENGTH);
484 goto err;
485 }
486 ret = BN_GF2m_mod_mul_arr(r, a, b, arr, ctx);
487 bn_check_top(r);
488 err:
489 OPENSSL_free(arr);
490 return ret;
491 }
492
493 /* Square a, reduce the result mod p, and store it in a. r could be a. */
BN_GF2m_mod_sqr_arr(BIGNUM * r,const BIGNUM * a,const int p[],BN_CTX * ctx)494 int BN_GF2m_mod_sqr_arr(BIGNUM *r, const BIGNUM *a, const int p[],
495 BN_CTX *ctx)
496 {
497 int i, ret = 0;
498 BIGNUM *s;
499
500 bn_check_top(a);
501 BN_CTX_start(ctx);
502 if ((s = BN_CTX_get(ctx)) == NULL)
503 goto err;
504 if (!bn_wexpand(s, 2 * a->top))
505 goto err;
506
507 for (i = a->top - 1; i >= 0; i--) {
508 s->d[2 * i + 1] = SQR1(a->d[i]);
509 s->d[2 * i] = SQR0(a->d[i]);
510 }
511
512 s->top = 2 * a->top;
513 bn_correct_top(s);
514 if (!BN_GF2m_mod_arr(r, s, p))
515 goto err;
516 bn_check_top(r);
517 ret = 1;
518 err:
519 BN_CTX_end(ctx);
520 return ret;
521 }
522
523 /*
524 * Square a, reduce the result mod p, and store it in a. r could be a. This
525 * function calls down to the BN_GF2m_mod_sqr_arr implementation; this
526 * wrapper function is only provided for convenience; for best performance,
527 * use the BN_GF2m_mod_sqr_arr function.
528 */
BN_GF2m_mod_sqr(BIGNUM * r,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)529 int BN_GF2m_mod_sqr(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
530 {
531 int ret = 0;
532 const int max = BN_num_bits(p) + 1;
533 int *arr;
534
535 bn_check_top(a);
536 bn_check_top(p);
537
538 arr = OPENSSL_malloc(sizeof(*arr) * max);
539 if (arr == NULL) {
540 ERR_raise(ERR_LIB_BN, ERR_R_MALLOC_FAILURE);
541 return 0;
542 }
543 ret = BN_GF2m_poly2arr(p, arr, max);
544 if (!ret || ret > max) {
545 ERR_raise(ERR_LIB_BN, BN_R_INVALID_LENGTH);
546 goto err;
547 }
548 ret = BN_GF2m_mod_sqr_arr(r, a, arr, ctx);
549 bn_check_top(r);
550 err:
551 OPENSSL_free(arr);
552 return ret;
553 }
554
555 /*
556 * Invert a, reduce modulo p, and store the result in r. r could be a. Uses
557 * Modified Almost Inverse Algorithm (Algorithm 10) from Hankerson, D.,
558 * Hernandez, J.L., and Menezes, A. "Software Implementation of Elliptic
559 * Curve Cryptography Over Binary Fields".
560 */
BN_GF2m_mod_inv_vartime(BIGNUM * r,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)561 static int BN_GF2m_mod_inv_vartime(BIGNUM *r, const BIGNUM *a,
562 const BIGNUM *p, BN_CTX *ctx)
563 {
564 BIGNUM *b, *c = NULL, *u = NULL, *v = NULL, *tmp;
565 int ret = 0;
566
567 bn_check_top(a);
568 bn_check_top(p);
569
570 BN_CTX_start(ctx);
571
572 b = BN_CTX_get(ctx);
573 c = BN_CTX_get(ctx);
574 u = BN_CTX_get(ctx);
575 v = BN_CTX_get(ctx);
576 if (v == NULL)
577 goto err;
578
579 if (!BN_GF2m_mod(u, a, p))
580 goto err;
581 if (BN_is_zero(u))
582 goto err;
583
584 if (!BN_copy(v, p))
585 goto err;
586 # if 0
587 if (!BN_one(b))
588 goto err;
589
590 while (1) {
591 while (!BN_is_odd(u)) {
592 if (BN_is_zero(u))
593 goto err;
594 if (!BN_rshift1(u, u))
595 goto err;
596 if (BN_is_odd(b)) {
597 if (!BN_GF2m_add(b, b, p))
598 goto err;
599 }
600 if (!BN_rshift1(b, b))
601 goto err;
602 }
603
604 if (BN_abs_is_word(u, 1))
605 break;
606
607 if (BN_num_bits(u) < BN_num_bits(v)) {
608 tmp = u;
609 u = v;
610 v = tmp;
611 tmp = b;
612 b = c;
613 c = tmp;
614 }
615
616 if (!BN_GF2m_add(u, u, v))
617 goto err;
618 if (!BN_GF2m_add(b, b, c))
619 goto err;
620 }
621 # else
622 {
623 int i;
624 int ubits = BN_num_bits(u);
625 int vbits = BN_num_bits(v); /* v is copy of p */
626 int top = p->top;
627 BN_ULONG *udp, *bdp, *vdp, *cdp;
628
629 if (!bn_wexpand(u, top))
630 goto err;
631 udp = u->d;
632 for (i = u->top; i < top; i++)
633 udp[i] = 0;
634 u->top = top;
635 if (!bn_wexpand(b, top))
636 goto err;
637 bdp = b->d;
638 bdp[0] = 1;
639 for (i = 1; i < top; i++)
640 bdp[i] = 0;
641 b->top = top;
642 if (!bn_wexpand(c, top))
643 goto err;
644 cdp = c->d;
645 for (i = 0; i < top; i++)
646 cdp[i] = 0;
647 c->top = top;
648 vdp = v->d; /* It pays off to "cache" *->d pointers,
649 * because it allows optimizer to be more
650 * aggressive. But we don't have to "cache"
651 * p->d, because *p is declared 'const'... */
652 while (1) {
653 while (ubits && !(udp[0] & 1)) {
654 BN_ULONG u0, u1, b0, b1, mask;
655
656 u0 = udp[0];
657 b0 = bdp[0];
658 mask = (BN_ULONG)0 - (b0 & 1);
659 b0 ^= p->d[0] & mask;
660 for (i = 0; i < top - 1; i++) {
661 u1 = udp[i + 1];
662 udp[i] = ((u0 >> 1) | (u1 << (BN_BITS2 - 1))) & BN_MASK2;
663 u0 = u1;
664 b1 = bdp[i + 1] ^ (p->d[i + 1] & mask);
665 bdp[i] = ((b0 >> 1) | (b1 << (BN_BITS2 - 1))) & BN_MASK2;
666 b0 = b1;
667 }
668 udp[i] = u0 >> 1;
669 bdp[i] = b0 >> 1;
670 ubits--;
671 }
672
673 if (ubits <= BN_BITS2) {
674 if (udp[0] == 0) /* poly was reducible */
675 goto err;
676 if (udp[0] == 1)
677 break;
678 }
679
680 if (ubits < vbits) {
681 i = ubits;
682 ubits = vbits;
683 vbits = i;
684 tmp = u;
685 u = v;
686 v = tmp;
687 tmp = b;
688 b = c;
689 c = tmp;
690 udp = vdp;
691 vdp = v->d;
692 bdp = cdp;
693 cdp = c->d;
694 }
695 for (i = 0; i < top; i++) {
696 udp[i] ^= vdp[i];
697 bdp[i] ^= cdp[i];
698 }
699 if (ubits == vbits) {
700 BN_ULONG ul;
701 int utop = (ubits - 1) / BN_BITS2;
702
703 while ((ul = udp[utop]) == 0 && utop)
704 utop--;
705 ubits = utop * BN_BITS2 + BN_num_bits_word(ul);
706 }
707 }
708 bn_correct_top(b);
709 }
710 # endif
711
712 if (!BN_copy(r, b))
713 goto err;
714 bn_check_top(r);
715 ret = 1;
716
717 err:
718 # ifdef BN_DEBUG
719 /* BN_CTX_end would complain about the expanded form */
720 bn_correct_top(c);
721 bn_correct_top(u);
722 bn_correct_top(v);
723 # endif
724 BN_CTX_end(ctx);
725 return ret;
726 }
727
728 /*-
729 * Wrapper for BN_GF2m_mod_inv_vartime that blinds the input before calling.
730 * This is not constant time.
731 * But it does eliminate first order deduction on the input.
732 */
BN_GF2m_mod_inv(BIGNUM * r,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)733 int BN_GF2m_mod_inv(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
734 {
735 BIGNUM *b = NULL;
736 int ret = 0;
737 int numbits;
738
739 BN_CTX_start(ctx);
740 if ((b = BN_CTX_get(ctx)) == NULL)
741 goto err;
742
743 /* Fail on a non-sensical input p value */
744 numbits = BN_num_bits(p);
745 if (numbits <= 1)
746 goto err;
747
748 /* generate blinding value */
749 do {
750 if (!BN_priv_rand_ex(b, numbits - 1,
751 BN_RAND_TOP_ANY, BN_RAND_BOTTOM_ANY, 0, ctx))
752 goto err;
753 } while (BN_is_zero(b));
754
755 /* r := a * b */
756 if (!BN_GF2m_mod_mul(r, a, b, p, ctx))
757 goto err;
758
759 /* r := 1/(a * b) */
760 if (!BN_GF2m_mod_inv_vartime(r, r, p, ctx))
761 goto err;
762
763 /* r := b/(a * b) = 1/a */
764 if (!BN_GF2m_mod_mul(r, r, b, p, ctx))
765 goto err;
766
767 ret = 1;
768
769 err:
770 BN_CTX_end(ctx);
771 return ret;
772 }
773
774 /*
775 * Invert xx, reduce modulo p, and store the result in r. r could be xx.
776 * This function calls down to the BN_GF2m_mod_inv implementation; this
777 * wrapper function is only provided for convenience; for best performance,
778 * use the BN_GF2m_mod_inv function.
779 */
BN_GF2m_mod_inv_arr(BIGNUM * r,const BIGNUM * xx,const int p[],BN_CTX * ctx)780 int BN_GF2m_mod_inv_arr(BIGNUM *r, const BIGNUM *xx, const int p[],
781 BN_CTX *ctx)
782 {
783 BIGNUM *field;
784 int ret = 0;
785
786 bn_check_top(xx);
787 BN_CTX_start(ctx);
788 if ((field = BN_CTX_get(ctx)) == NULL)
789 goto err;
790 if (!BN_GF2m_arr2poly(p, field))
791 goto err;
792
793 ret = BN_GF2m_mod_inv(r, xx, field, ctx);
794 bn_check_top(r);
795
796 err:
797 BN_CTX_end(ctx);
798 return ret;
799 }
800
801 /*
802 * Divide y by x, reduce modulo p, and store the result in r. r could be x
803 * or y, x could equal y.
804 */
BN_GF2m_mod_div(BIGNUM * r,const BIGNUM * y,const BIGNUM * x,const BIGNUM * p,BN_CTX * ctx)805 int BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x,
806 const BIGNUM *p, BN_CTX *ctx)
807 {
808 BIGNUM *xinv = NULL;
809 int ret = 0;
810
811 bn_check_top(y);
812 bn_check_top(x);
813 bn_check_top(p);
814
815 BN_CTX_start(ctx);
816 xinv = BN_CTX_get(ctx);
817 if (xinv == NULL)
818 goto err;
819
820 if (!BN_GF2m_mod_inv(xinv, x, p, ctx))
821 goto err;
822 if (!BN_GF2m_mod_mul(r, y, xinv, p, ctx))
823 goto err;
824 bn_check_top(r);
825 ret = 1;
826
827 err:
828 BN_CTX_end(ctx);
829 return ret;
830 }
831
832 /*
833 * Divide yy by xx, reduce modulo p, and store the result in r. r could be xx
834 * * or yy, xx could equal yy. This function calls down to the
835 * BN_GF2m_mod_div implementation; this wrapper function is only provided for
836 * convenience; for best performance, use the BN_GF2m_mod_div function.
837 */
BN_GF2m_mod_div_arr(BIGNUM * r,const BIGNUM * yy,const BIGNUM * xx,const int p[],BN_CTX * ctx)838 int BN_GF2m_mod_div_arr(BIGNUM *r, const BIGNUM *yy, const BIGNUM *xx,
839 const int p[], BN_CTX *ctx)
840 {
841 BIGNUM *field;
842 int ret = 0;
843
844 bn_check_top(yy);
845 bn_check_top(xx);
846
847 BN_CTX_start(ctx);
848 if ((field = BN_CTX_get(ctx)) == NULL)
849 goto err;
850 if (!BN_GF2m_arr2poly(p, field))
851 goto err;
852
853 ret = BN_GF2m_mod_div(r, yy, xx, field, ctx);
854 bn_check_top(r);
855
856 err:
857 BN_CTX_end(ctx);
858 return ret;
859 }
860
861 /*
862 * Compute the bth power of a, reduce modulo p, and store the result in r. r
863 * could be a. Uses simple square-and-multiply algorithm A.5.1 from IEEE
864 * P1363.
865 */
BN_GF2m_mod_exp_arr(BIGNUM * r,const BIGNUM * a,const BIGNUM * b,const int p[],BN_CTX * ctx)866 int BN_GF2m_mod_exp_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
867 const int p[], BN_CTX *ctx)
868 {
869 int ret = 0, i, n;
870 BIGNUM *u;
871
872 bn_check_top(a);
873 bn_check_top(b);
874
875 if (BN_is_zero(b))
876 return BN_one(r);
877
878 if (BN_abs_is_word(b, 1))
879 return (BN_copy(r, a) != NULL);
880
881 BN_CTX_start(ctx);
882 if ((u = BN_CTX_get(ctx)) == NULL)
883 goto err;
884
885 if (!BN_GF2m_mod_arr(u, a, p))
886 goto err;
887
888 n = BN_num_bits(b) - 1;
889 for (i = n - 1; i >= 0; i--) {
890 if (!BN_GF2m_mod_sqr_arr(u, u, p, ctx))
891 goto err;
892 if (BN_is_bit_set(b, i)) {
893 if (!BN_GF2m_mod_mul_arr(u, u, a, p, ctx))
894 goto err;
895 }
896 }
897 if (!BN_copy(r, u))
898 goto err;
899 bn_check_top(r);
900 ret = 1;
901 err:
902 BN_CTX_end(ctx);
903 return ret;
904 }
905
906 /*
907 * Compute the bth power of a, reduce modulo p, and store the result in r. r
908 * could be a. This function calls down to the BN_GF2m_mod_exp_arr
909 * implementation; this wrapper function is only provided for convenience;
910 * for best performance, use the BN_GF2m_mod_exp_arr function.
911 */
BN_GF2m_mod_exp(BIGNUM * r,const BIGNUM * a,const BIGNUM * b,const BIGNUM * p,BN_CTX * ctx)912 int BN_GF2m_mod_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
913 const BIGNUM *p, BN_CTX *ctx)
914 {
915 int ret = 0;
916 const int max = BN_num_bits(p) + 1;
917 int *arr;
918
919 bn_check_top(a);
920 bn_check_top(b);
921 bn_check_top(p);
922
923 arr = OPENSSL_malloc(sizeof(*arr) * max);
924 if (arr == NULL) {
925 ERR_raise(ERR_LIB_BN, ERR_R_MALLOC_FAILURE);
926 return 0;
927 }
928 ret = BN_GF2m_poly2arr(p, arr, max);
929 if (!ret || ret > max) {
930 ERR_raise(ERR_LIB_BN, BN_R_INVALID_LENGTH);
931 goto err;
932 }
933 ret = BN_GF2m_mod_exp_arr(r, a, b, arr, ctx);
934 bn_check_top(r);
935 err:
936 OPENSSL_free(arr);
937 return ret;
938 }
939
940 /*
941 * Compute the square root of a, reduce modulo p, and store the result in r.
942 * r could be a. Uses exponentiation as in algorithm A.4.1 from IEEE P1363.
943 */
BN_GF2m_mod_sqrt_arr(BIGNUM * r,const BIGNUM * a,const int p[],BN_CTX * ctx)944 int BN_GF2m_mod_sqrt_arr(BIGNUM *r, const BIGNUM *a, const int p[],
945 BN_CTX *ctx)
946 {
947 int ret = 0;
948 BIGNUM *u;
949
950 bn_check_top(a);
951
952 if (p[0] == 0) {
953 /* reduction mod 1 => return 0 */
954 BN_zero(r);
955 return 1;
956 }
957
958 BN_CTX_start(ctx);
959 if ((u = BN_CTX_get(ctx)) == NULL)
960 goto err;
961
962 if (!BN_set_bit(u, p[0] - 1))
963 goto err;
964 ret = BN_GF2m_mod_exp_arr(r, a, u, p, ctx);
965 bn_check_top(r);
966
967 err:
968 BN_CTX_end(ctx);
969 return ret;
970 }
971
972 /*
973 * Compute the square root of a, reduce modulo p, and store the result in r.
974 * r could be a. This function calls down to the BN_GF2m_mod_sqrt_arr
975 * implementation; this wrapper function is only provided for convenience;
976 * for best performance, use the BN_GF2m_mod_sqrt_arr function.
977 */
BN_GF2m_mod_sqrt(BIGNUM * r,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)978 int BN_GF2m_mod_sqrt(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
979 {
980 int ret = 0;
981 const int max = BN_num_bits(p) + 1;
982 int *arr;
983
984 bn_check_top(a);
985 bn_check_top(p);
986
987 arr = OPENSSL_malloc(sizeof(*arr) * max);
988 if (arr == NULL) {
989 ERR_raise(ERR_LIB_BN, ERR_R_MALLOC_FAILURE);
990 return 0;
991 }
992 ret = BN_GF2m_poly2arr(p, arr, max);
993 if (!ret || ret > max) {
994 ERR_raise(ERR_LIB_BN, BN_R_INVALID_LENGTH);
995 goto err;
996 }
997 ret = BN_GF2m_mod_sqrt_arr(r, a, arr, ctx);
998 bn_check_top(r);
999 err:
1000 OPENSSL_free(arr);
1001 return ret;
1002 }
1003
1004 /*
1005 * Find r such that r^2 + r = a mod p. r could be a. If no r exists returns
1006 * 0. Uses algorithms A.4.7 and A.4.6 from IEEE P1363.
1007 */
BN_GF2m_mod_solve_quad_arr(BIGNUM * r,const BIGNUM * a_,const int p[],BN_CTX * ctx)1008 int BN_GF2m_mod_solve_quad_arr(BIGNUM *r, const BIGNUM *a_, const int p[],
1009 BN_CTX *ctx)
1010 {
1011 int ret = 0, count = 0, j;
1012 BIGNUM *a, *z, *rho, *w, *w2, *tmp;
1013
1014 bn_check_top(a_);
1015
1016 if (p[0] == 0) {
1017 /* reduction mod 1 => return 0 */
1018 BN_zero(r);
1019 return 1;
1020 }
1021
1022 BN_CTX_start(ctx);
1023 a = BN_CTX_get(ctx);
1024 z = BN_CTX_get(ctx);
1025 w = BN_CTX_get(ctx);
1026 if (w == NULL)
1027 goto err;
1028
1029 if (!BN_GF2m_mod_arr(a, a_, p))
1030 goto err;
1031
1032 if (BN_is_zero(a)) {
1033 BN_zero(r);
1034 ret = 1;
1035 goto err;
1036 }
1037
1038 if (p[0] & 0x1) { /* m is odd */
1039 /* compute half-trace of a */
1040 if (!BN_copy(z, a))
1041 goto err;
1042 for (j = 1; j <= (p[0] - 1) / 2; j++) {
1043 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1044 goto err;
1045 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1046 goto err;
1047 if (!BN_GF2m_add(z, z, a))
1048 goto err;
1049 }
1050
1051 } else { /* m is even */
1052
1053 rho = BN_CTX_get(ctx);
1054 w2 = BN_CTX_get(ctx);
1055 tmp = BN_CTX_get(ctx);
1056 if (tmp == NULL)
1057 goto err;
1058 do {
1059 if (!BN_priv_rand_ex(rho, p[0], BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ANY,
1060 0, ctx))
1061 goto err;
1062 if (!BN_GF2m_mod_arr(rho, rho, p))
1063 goto err;
1064 BN_zero(z);
1065 if (!BN_copy(w, rho))
1066 goto err;
1067 for (j = 1; j <= p[0] - 1; j++) {
1068 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1069 goto err;
1070 if (!BN_GF2m_mod_sqr_arr(w2, w, p, ctx))
1071 goto err;
1072 if (!BN_GF2m_mod_mul_arr(tmp, w2, a, p, ctx))
1073 goto err;
1074 if (!BN_GF2m_add(z, z, tmp))
1075 goto err;
1076 if (!BN_GF2m_add(w, w2, rho))
1077 goto err;
1078 }
1079 count++;
1080 } while (BN_is_zero(w) && (count < MAX_ITERATIONS));
1081 if (BN_is_zero(w)) {
1082 ERR_raise(ERR_LIB_BN, BN_R_TOO_MANY_ITERATIONS);
1083 goto err;
1084 }
1085 }
1086
1087 if (!BN_GF2m_mod_sqr_arr(w, z, p, ctx))
1088 goto err;
1089 if (!BN_GF2m_add(w, z, w))
1090 goto err;
1091 if (BN_GF2m_cmp(w, a)) {
1092 ERR_raise(ERR_LIB_BN, BN_R_NO_SOLUTION);
1093 goto err;
1094 }
1095
1096 if (!BN_copy(r, z))
1097 goto err;
1098 bn_check_top(r);
1099
1100 ret = 1;
1101
1102 err:
1103 BN_CTX_end(ctx);
1104 return ret;
1105 }
1106
1107 /*
1108 * Find r such that r^2 + r = a mod p. r could be a. If no r exists returns
1109 * 0. This function calls down to the BN_GF2m_mod_solve_quad_arr
1110 * implementation; this wrapper function is only provided for convenience;
1111 * for best performance, use the BN_GF2m_mod_solve_quad_arr function.
1112 */
BN_GF2m_mod_solve_quad(BIGNUM * r,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)1113 int BN_GF2m_mod_solve_quad(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
1114 BN_CTX *ctx)
1115 {
1116 int ret = 0;
1117 const int max = BN_num_bits(p) + 1;
1118 int *arr;
1119
1120 bn_check_top(a);
1121 bn_check_top(p);
1122
1123 arr = OPENSSL_malloc(sizeof(*arr) * max);
1124 if (arr == NULL) {
1125 ERR_raise(ERR_LIB_BN, ERR_R_MALLOC_FAILURE);
1126 goto err;
1127 }
1128 ret = BN_GF2m_poly2arr(p, arr, max);
1129 if (!ret || ret > max) {
1130 ERR_raise(ERR_LIB_BN, BN_R_INVALID_LENGTH);
1131 goto err;
1132 }
1133 ret = BN_GF2m_mod_solve_quad_arr(r, a, arr, ctx);
1134 bn_check_top(r);
1135 err:
1136 OPENSSL_free(arr);
1137 return ret;
1138 }
1139
1140 /*
1141 * Convert the bit-string representation of a polynomial ( \sum_{i=0}^n a_i *
1142 * x^i) into an array of integers corresponding to the bits with non-zero
1143 * coefficient. Array is terminated with -1. Up to max elements of the array
1144 * will be filled. Return value is total number of array elements that would
1145 * be filled if array was large enough.
1146 */
BN_GF2m_poly2arr(const BIGNUM * a,int p[],int max)1147 int BN_GF2m_poly2arr(const BIGNUM *a, int p[], int max)
1148 {
1149 int i, j, k = 0;
1150 BN_ULONG mask;
1151
1152 if (BN_is_zero(a))
1153 return 0;
1154
1155 for (i = a->top - 1; i >= 0; i--) {
1156 if (!a->d[i])
1157 /* skip word if a->d[i] == 0 */
1158 continue;
1159 mask = BN_TBIT;
1160 for (j = BN_BITS2 - 1; j >= 0; j--) {
1161 if (a->d[i] & mask) {
1162 if (k < max)
1163 p[k] = BN_BITS2 * i + j;
1164 k++;
1165 }
1166 mask >>= 1;
1167 }
1168 }
1169
1170 if (k < max) {
1171 p[k] = -1;
1172 k++;
1173 }
1174
1175 return k;
1176 }
1177
1178 /*
1179 * Convert the coefficient array representation of a polynomial to a
1180 * bit-string. The array must be terminated by -1.
1181 */
BN_GF2m_arr2poly(const int p[],BIGNUM * a)1182 int BN_GF2m_arr2poly(const int p[], BIGNUM *a)
1183 {
1184 int i;
1185
1186 bn_check_top(a);
1187 BN_zero(a);
1188 for (i = 0; p[i] != -1; i++) {
1189 if (BN_set_bit(a, p[i]) == 0)
1190 return 0;
1191 }
1192 bn_check_top(a);
1193
1194 return 1;
1195 }
1196
1197 #endif
1198