1/*
2 * Copyright 2002-2018 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
60static 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)
129static 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 */
231static 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
244void 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 */
252int 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. */
292int 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 */
390int 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 BNerr(BN_F_BN_GF2M_MOD, 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 */
410int 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 */
465int 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 = NULL;
471 bn_check_top(a);
472 bn_check_top(b);
473 bn_check_top(p);
474 if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
475 goto err;
476 ret = BN_GF2m_poly2arr(p, arr, max);
477 if (!ret || ret > max) {
478 BNerr(BN_F_BN_GF2M_MOD_MUL, BN_R_INVALID_LENGTH);
479 goto err;
480 }
481 ret = BN_GF2m_mod_mul_arr(r, a, b, arr, ctx);
482 bn_check_top(r);
483 err:
484 OPENSSL_free(arr);
485 return ret;
486}
487
488/* Square a, reduce the result mod p, and store it in a. r could be a. */
489int BN_GF2m_mod_sqr_arr(BIGNUM *r, const BIGNUM *a, const int p[],
490 BN_CTX *ctx)
491{
492 int i, ret = 0;
493 BIGNUM *s;
494
495 bn_check_top(a);
496 BN_CTX_start(ctx);
497 if ((s = BN_CTX_get(ctx)) == NULL)
498 goto err;
499 if (!bn_wexpand(s, 2 * a->top))
500 goto err;
501
502 for (i = a->top - 1; i >= 0; i--) {
503 s->d[2 * i + 1] = SQR1(a->d[i]);
504 s->d[2 * i] = SQR0(a->d[i]);
505 }
506
507 s->top = 2 * a->top;
508 bn_correct_top(s);
509 if (!BN_GF2m_mod_arr(r, s, p))
510 goto err;
511 bn_check_top(r);
512 ret = 1;
513 err:
514 BN_CTX_end(ctx);
515 return ret;
516}
517
518/*
519 * Square a, reduce the result mod p, and store it in a. r could be a. This
520 * function calls down to the BN_GF2m_mod_sqr_arr implementation; this
521 * wrapper function is only provided for convenience; for best performance,
522 * use the BN_GF2m_mod_sqr_arr function.
523 */
524int BN_GF2m_mod_sqr(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
525{
526 int ret = 0;
527 const int max = BN_num_bits(p) + 1;
528 int *arr = NULL;
529
530 bn_check_top(a);
531 bn_check_top(p);
532 if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
533 goto err;
534 ret = BN_GF2m_poly2arr(p, arr, max);
535 if (!ret || ret > max) {
536 BNerr(BN_F_BN_GF2M_MOD_SQR, BN_R_INVALID_LENGTH);
537 goto err;
538 }
539 ret = BN_GF2m_mod_sqr_arr(r, a, arr, ctx);
540 bn_check_top(r);
541 err:
542 OPENSSL_free(arr);
543 return ret;
544}
545
546/*
547 * Invert a, reduce modulo p, and store the result in r. r could be a. Uses
548 * Modified Almost Inverse Algorithm (Algorithm 10) from Hankerson, D.,
549 * Hernandez, J.L., and Menezes, A. "Software Implementation of Elliptic
550 * Curve Cryptography Over Binary Fields".
551 */
552static int BN_GF2m_mod_inv_vartime(BIGNUM *r, const BIGNUM *a,
553 const BIGNUM *p, BN_CTX *ctx)
554{
555 BIGNUM *b, *c = NULL, *u = NULL, *v = NULL, *tmp;
556 int ret = 0;
557
558 bn_check_top(a);
559 bn_check_top(p);
560
561 BN_CTX_start(ctx);
562
563 b = BN_CTX_get(ctx);
564 c = BN_CTX_get(ctx);
565 u = BN_CTX_get(ctx);
566 v = BN_CTX_get(ctx);
567 if (v == NULL)
568 goto err;
569
570 if (!BN_GF2m_mod(u, a, p))
571 goto err;
572 if (BN_is_zero(u))
573 goto err;
574
575 if (!BN_copy(v, p))
576 goto err;
577# if 0
578 if (!BN_one(b))
579 goto err;
580
581 while (1) {
582 while (!BN_is_odd(u)) {
583 if (BN_is_zero(u))
584 goto err;
585 if (!BN_rshift1(u, u))
586 goto err;
587 if (BN_is_odd(b)) {
588 if (!BN_GF2m_add(b, b, p))
589 goto err;
590 }
591 if (!BN_rshift1(b, b))
592 goto err;
593 }
594
595 if (BN_abs_is_word(u, 1))
596 break;
597
598 if (BN_num_bits(u) < BN_num_bits(v)) {
599 tmp = u;
600 u = v;
601 v = tmp;
602 tmp = b;
603 b = c;
604 c = tmp;
605 }
606
607 if (!BN_GF2m_add(u, u, v))
608 goto err;
609 if (!BN_GF2m_add(b, b, c))
610 goto err;
611 }
612# else
613 {
614 int i;
615 int ubits = BN_num_bits(u);
616 int vbits = BN_num_bits(v); /* v is copy of p */
617 int top = p->top;
618 BN_ULONG *udp, *bdp, *vdp, *cdp;
619
620 if (!bn_wexpand(u, top))
621 goto err;
622 udp = u->d;
623 for (i = u->top; i < top; i++)
624 udp[i] = 0;
625 u->top = top;
626 if (!bn_wexpand(b, top))
627 goto err;
628 bdp = b->d;
629 bdp[0] = 1;
630 for (i = 1; i < top; i++)
631 bdp[i] = 0;
632 b->top = top;
633 if (!bn_wexpand(c, top))
634 goto err;
635 cdp = c->d;
636 for (i = 0; i < top; i++)
637 cdp[i] = 0;
638 c->top = top;
639 vdp = v->d; /* It pays off to "cache" *->d pointers,
640 * because it allows optimizer to be more
641 * aggressive. But we don't have to "cache"
642 * p->d, because *p is declared 'const'... */
643 while (1) {
644 while (ubits && !(udp[0] & 1)) {
645 BN_ULONG u0, u1, b0, b1, mask;
646
647 u0 = udp[0];
648 b0 = bdp[0];
649 mask = (BN_ULONG)0 - (b0 & 1);
650 b0 ^= p->d[0] & mask;
651 for (i = 0; i < top - 1; i++) {
652 u1 = udp[i + 1];
653 udp[i] = ((u0 >> 1) | (u1 << (BN_BITS2 - 1))) & BN_MASK2;
654 u0 = u1;
655 b1 = bdp[i + 1] ^ (p->d[i + 1] & mask);
656 bdp[i] = ((b0 >> 1) | (b1 << (BN_BITS2 - 1))) & BN_MASK2;
657 b0 = b1;
658 }
659 udp[i] = u0 >> 1;
660 bdp[i] = b0 >> 1;
661 ubits--;
662 }
663
664 if (ubits <= BN_BITS2) {
665 if (udp[0] == 0) /* poly was reducible */
666 goto err;
667 if (udp[0] == 1)
668 break;
669 }
670
671 if (ubits < vbits) {
672 i = ubits;
673 ubits = vbits;
674 vbits = i;
675 tmp = u;
676 u = v;
677 v = tmp;
678 tmp = b;
679 b = c;
680 c = tmp;
681 udp = vdp;
682 vdp = v->d;
683 bdp = cdp;
684 cdp = c->d;
685 }
686 for (i = 0; i < top; i++) {
687 udp[i] ^= vdp[i];
688 bdp[i] ^= cdp[i];
689 }
690 if (ubits == vbits) {
691 BN_ULONG ul;
692 int utop = (ubits - 1) / BN_BITS2;
693
694 while ((ul = udp[utop]) == 0 && utop)
695 utop--;
696 ubits = utop * BN_BITS2 + BN_num_bits_word(ul);
697 }
698 }
699 bn_correct_top(b);
700 }
701# endif
702
703 if (!BN_copy(r, b))
704 goto err;
705 bn_check_top(r);
706 ret = 1;
707
708 err:
709# ifdef BN_DEBUG /* BN_CTX_end would complain about the
710 * expanded form */
711 bn_correct_top(c);
712 bn_correct_top(u);
713 bn_correct_top(v);
714# endif
715 BN_CTX_end(ctx);
716 return ret;
717}
718
719/*-
720 * Wrapper for BN_GF2m_mod_inv_vartime that blinds the input before calling.
721 * This is not constant time.
722 * But it does eliminate first order deduction on the input.
723 */
724int BN_GF2m_mod_inv(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
725{
726 BIGNUM *b = NULL;
727 int ret = 0;
728
729 BN_CTX_start(ctx);
730 if ((b = BN_CTX_get(ctx)) == NULL)
731 goto err;
732
733 /* generate blinding value */
734 do {
735 if (!BN_priv_rand_ex(b, BN_num_bits(p) - 1,
736 BN_RAND_TOP_ANY, BN_RAND_BOTTOM_ANY, ctx))
737 goto err;
738 } while (BN_is_zero(b));
739
740 /* r := a * b */
741 if (!BN_GF2m_mod_mul(r, a, b, p, ctx))
742 goto err;
743
744 /* r := 1/(a * b) */
745 if (!BN_GF2m_mod_inv_vartime(r, r, p, ctx))
746 goto err;
747
748 /* r := b/(a * b) = 1/a */
749 if (!BN_GF2m_mod_mul(r, r, b, p, ctx))
750 goto err;
751
752 ret = 1;
753
754 err:
755 BN_CTX_end(ctx);
756 return ret;
757}
758
759/*
760 * Invert xx, reduce modulo p, and store the result in r. r could be xx.
761 * This function calls down to the BN_GF2m_mod_inv implementation; this
762 * wrapper function is only provided for convenience; for best performance,
763 * use the BN_GF2m_mod_inv function.
764 */
765int BN_GF2m_mod_inv_arr(BIGNUM *r, const BIGNUM *xx, const int p[],
766 BN_CTX *ctx)
767{
768 BIGNUM *field;
769 int ret = 0;
770
771 bn_check_top(xx);
772 BN_CTX_start(ctx);
773 if ((field = BN_CTX_get(ctx)) == NULL)
774 goto err;
775 if (!BN_GF2m_arr2poly(p, field))
776 goto err;
777
778 ret = BN_GF2m_mod_inv(r, xx, field, ctx);
779 bn_check_top(r);
780
781 err:
782 BN_CTX_end(ctx);
783 return ret;
784}
785
786/*
787 * Divide y by x, reduce modulo p, and store the result in r. r could be x
788 * or y, x could equal y.
789 */
790int BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x,
791 const BIGNUM *p, BN_CTX *ctx)
792{
793 BIGNUM *xinv = NULL;
794 int ret = 0;
795
796 bn_check_top(y);
797 bn_check_top(x);
798 bn_check_top(p);
799
800 BN_CTX_start(ctx);
801 xinv = BN_CTX_get(ctx);
802 if (xinv == NULL)
803 goto err;
804
805 if (!BN_GF2m_mod_inv(xinv, x, p, ctx))
806 goto err;
807 if (!BN_GF2m_mod_mul(r, y, xinv, p, ctx))
808 goto err;
809 bn_check_top(r);
810 ret = 1;
811
812 err:
813 BN_CTX_end(ctx);
814 return ret;
815}
816
817/*
818 * Divide yy by xx, reduce modulo p, and store the result in r. r could be xx
819 * * or yy, xx could equal yy. This function calls down to the
820 * BN_GF2m_mod_div implementation; this wrapper function is only provided for
821 * convenience; for best performance, use the BN_GF2m_mod_div function.
822 */
823int BN_GF2m_mod_div_arr(BIGNUM *r, const BIGNUM *yy, const BIGNUM *xx,
824 const int p[], BN_CTX *ctx)
825{
826 BIGNUM *field;
827 int ret = 0;
828
829 bn_check_top(yy);
830 bn_check_top(xx);
831
832 BN_CTX_start(ctx);
833 if ((field = BN_CTX_get(ctx)) == NULL)
834 goto err;
835 if (!BN_GF2m_arr2poly(p, field))
836 goto err;
837
838 ret = BN_GF2m_mod_div(r, yy, xx, field, ctx);
839 bn_check_top(r);
840
841 err:
842 BN_CTX_end(ctx);
843 return ret;
844}
845
846/*
847 * Compute the bth power of a, reduce modulo p, and store the result in r. r
848 * could be a. Uses simple square-and-multiply algorithm A.5.1 from IEEE
849 * P1363.
850 */
851int BN_GF2m_mod_exp_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
852 const int p[], BN_CTX *ctx)
853{
854 int ret = 0, i, n;
855 BIGNUM *u;
856
857 bn_check_top(a);
858 bn_check_top(b);
859
860 if (BN_is_zero(b))
861 return BN_one(r);
862
863 if (BN_abs_is_word(b, 1))
864 return (BN_copy(r, a) != NULL);
865
866 BN_CTX_start(ctx);
867 if ((u = BN_CTX_get(ctx)) == NULL)
868 goto err;
869
870 if (!BN_GF2m_mod_arr(u, a, p))
871 goto err;
872
873 n = BN_num_bits(b) - 1;
874 for (i = n - 1; i >= 0; i--) {
875 if (!BN_GF2m_mod_sqr_arr(u, u, p, ctx))
876 goto err;
877 if (BN_is_bit_set(b, i)) {
878 if (!BN_GF2m_mod_mul_arr(u, u, a, p, ctx))
879 goto err;
880 }
881 }
882 if (!BN_copy(r, u))
883 goto err;
884 bn_check_top(r);
885 ret = 1;
886 err:
887 BN_CTX_end(ctx);
888 return ret;
889}
890
891/*
892 * Compute the bth power of a, reduce modulo p, and store the result in r. r
893 * could be a. This function calls down to the BN_GF2m_mod_exp_arr
894 * implementation; this wrapper function is only provided for convenience;
895 * for best performance, use the BN_GF2m_mod_exp_arr function.
896 */
897int BN_GF2m_mod_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
898 const BIGNUM *p, BN_CTX *ctx)
899{
900 int ret = 0;
901 const int max = BN_num_bits(p) + 1;
902 int *arr = NULL;
903 bn_check_top(a);
904 bn_check_top(b);
905 bn_check_top(p);
906 if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
907 goto err;
908 ret = BN_GF2m_poly2arr(p, arr, max);
909 if (!ret || ret > max) {
910 BNerr(BN_F_BN_GF2M_MOD_EXP, BN_R_INVALID_LENGTH);
911 goto err;
912 }
913 ret = BN_GF2m_mod_exp_arr(r, a, b, arr, ctx);
914 bn_check_top(r);
915 err:
916 OPENSSL_free(arr);
917 return ret;
918}
919
920/*
921 * Compute the square root of a, reduce modulo p, and store the result in r.
922 * r could be a. Uses exponentiation as in algorithm A.4.1 from IEEE P1363.
923 */
924int BN_GF2m_mod_sqrt_arr(BIGNUM *r, const BIGNUM *a, const int p[],
925 BN_CTX *ctx)
926{
927 int ret = 0;
928 BIGNUM *u;
929
930 bn_check_top(a);
931
932 if (p[0] == 0) {
933 /* reduction mod 1 => return 0 */
934 BN_zero(r);
935 return 1;
936 }
937
938 BN_CTX_start(ctx);
939 if ((u = BN_CTX_get(ctx)) == NULL)
940 goto err;
941
942 if (!BN_set_bit(u, p[0] - 1))
943 goto err;
944 ret = BN_GF2m_mod_exp_arr(r, a, u, p, ctx);
945 bn_check_top(r);
946
947 err:
948 BN_CTX_end(ctx);
949 return ret;
950}
951
952/*
953 * Compute the square root of a, reduce modulo p, and store the result in r.
954 * r could be a. This function calls down to the BN_GF2m_mod_sqrt_arr
955 * implementation; this wrapper function is only provided for convenience;
956 * for best performance, use the BN_GF2m_mod_sqrt_arr function.
957 */
958int BN_GF2m_mod_sqrt(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
959{
960 int ret = 0;
961 const int max = BN_num_bits(p) + 1;
962 int *arr = NULL;
963 bn_check_top(a);
964 bn_check_top(p);
965 if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
966 goto err;
967 ret = BN_GF2m_poly2arr(p, arr, max);
968 if (!ret || ret > max) {
969 BNerr(BN_F_BN_GF2M_MOD_SQRT, BN_R_INVALID_LENGTH);
970 goto err;
971 }
972 ret = BN_GF2m_mod_sqrt_arr(r, a, arr, ctx);
973 bn_check_top(r);
974 err:
975 OPENSSL_free(arr);
976 return ret;
977}
978
979/*
980 * Find r such that r^2 + r = a mod p. r could be a. If no r exists returns
981 * 0. Uses algorithms A.4.7 and A.4.6 from IEEE P1363.
982 */
983int BN_GF2m_mod_solve_quad_arr(BIGNUM *r, const BIGNUM *a_, const int p[],
984 BN_CTX *ctx)
985{
986 int ret = 0, count = 0, j;
987 BIGNUM *a, *z, *rho, *w, *w2, *tmp;
988
989 bn_check_top(a_);
990
991 if (p[0] == 0) {
992 /* reduction mod 1 => return 0 */
993 BN_zero(r);
994 return 1;
995 }
996
997 BN_CTX_start(ctx);
998 a = BN_CTX_get(ctx);
999 z = BN_CTX_get(ctx);
1000 w = BN_CTX_get(ctx);
1001 if (w == NULL)
1002 goto err;
1003
1004 if (!BN_GF2m_mod_arr(a, a_, p))
1005 goto err;
1006
1007 if (BN_is_zero(a)) {
1008 BN_zero(r);
1009 ret = 1;
1010 goto err;
1011 }
1012
1013 if (p[0] & 0x1) { /* m is odd */
1014 /* compute half-trace of a */
1015 if (!BN_copy(z, a))
1016 goto err;
1017 for (j = 1; j <= (p[0] - 1) / 2; j++) {
1018 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1019 goto err;
1020 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1021 goto err;
1022 if (!BN_GF2m_add(z, z, a))
1023 goto err;
1024 }
1025
1026 } else { /* m is even */
1027
1028 rho = BN_CTX_get(ctx);
1029 w2 = BN_CTX_get(ctx);
1030 tmp = BN_CTX_get(ctx);
1031 if (tmp == NULL)
1032 goto err;
1033 do {
1034 if (!BN_priv_rand_ex(rho, p[0], BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ANY,
1035 ctx))
1036 goto err;
1037 if (!BN_GF2m_mod_arr(rho, rho, p))
1038 goto err;
1039 BN_zero(z);
1040 if (!BN_copy(w, rho))
1041 goto err;
1042 for (j = 1; j <= p[0] - 1; j++) {
1043 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1044 goto err;
1045 if (!BN_GF2m_mod_sqr_arr(w2, w, p, ctx))
1046 goto err;
1047 if (!BN_GF2m_mod_mul_arr(tmp, w2, a, p, ctx))
1048 goto err;
1049 if (!BN_GF2m_add(z, z, tmp))
1050 goto err;
1051 if (!BN_GF2m_add(w, w2, rho))
1052 goto err;
1053 }
1054 count++;
1055 } while (BN_is_zero(w) && (count < MAX_ITERATIONS));
1056 if (BN_is_zero(w)) {
1057 BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_TOO_MANY_ITERATIONS);
1058 goto err;
1059 }
1060 }
1061
1062 if (!BN_GF2m_mod_sqr_arr(w, z, p, ctx))
1063 goto err;
1064 if (!BN_GF2m_add(w, z, w))
1065 goto err;
1066 if (BN_GF2m_cmp(w, a)) {
1067 BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_NO_SOLUTION);
1068 goto err;
1069 }
1070
1071 if (!BN_copy(r, z))
1072 goto err;
1073 bn_check_top(r);
1074
1075 ret = 1;
1076
1077 err:
1078 BN_CTX_end(ctx);
1079 return ret;
1080}
1081
1082/*
1083 * Find r such that r^2 + r = a mod p. r could be a. If no r exists returns
1084 * 0. This function calls down to the BN_GF2m_mod_solve_quad_arr
1085 * implementation; this wrapper function is only provided for convenience;
1086 * for best performance, use the BN_GF2m_mod_solve_quad_arr function.
1087 */
1088int BN_GF2m_mod_solve_quad(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
1089 BN_CTX *ctx)
1090{
1091 int ret = 0;
1092 const int max = BN_num_bits(p) + 1;
1093 int *arr = NULL;
1094 bn_check_top(a);
1095 bn_check_top(p);
1096 if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
1097 goto err;
1098 ret = BN_GF2m_poly2arr(p, arr, max);
1099 if (!ret || ret > max) {
1100 BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD, BN_R_INVALID_LENGTH);
1101 goto err;
1102 }
1103 ret = BN_GF2m_mod_solve_quad_arr(r, a, arr, ctx);
1104 bn_check_top(r);
1105 err:
1106 OPENSSL_free(arr);
1107 return ret;
1108}
1109
1110/*
1111 * Convert the bit-string representation of a polynomial ( \sum_{i=0}^n a_i *
1112 * x^i) into an array of integers corresponding to the bits with non-zero
1113 * coefficient. Array is terminated with -1. Up to max elements of the array
1114 * will be filled. Return value is total number of array elements that would
1115 * be filled if array was large enough.
1116 */
1117int BN_GF2m_poly2arr(const BIGNUM *a, int p[], int max)
1118{
1119 int i, j, k = 0;
1120 BN_ULONG mask;
1121
1122 if (BN_is_zero(a))
1123 return 0;
1124
1125 for (i = a->top - 1; i >= 0; i--) {
1126 if (!a->d[i])
1127 /* skip word if a->d[i] == 0 */
1128 continue;
1129 mask = BN_TBIT;
1130 for (j = BN_BITS2 - 1; j >= 0; j--) {
1131 if (a->d[i] & mask) {
1132 if (k < max)
1133 p[k] = BN_BITS2 * i + j;
1134 k++;
1135 }
1136 mask >>= 1;
1137 }
1138 }
1139
1140 if (k < max) {
1141 p[k] = -1;
1142 k++;
1143 }
1144
1145 return k;
1146}
1147
1148/*
1149 * Convert the coefficient array representation of a polynomial to a
1150 * bit-string. The array must be terminated by -1.
1151 */
1152int BN_GF2m_arr2poly(const int p[], BIGNUM *a)
1153{
1154 int i;
1155
1156 bn_check_top(a);
1157 BN_zero(a);
1158 for (i = 0; p[i] != -1; i++) {
1159 if (BN_set_bit(a, p[i]) == 0)
1160 return 0;
1161 }
1162 bn_check_top(a);
1163
1164 return 1;
1165}
1166
1167#endif
1168