1 /* $OpenBSD: bn_gf2m.c,v 1.23 2017/01/29 17:49:22 beck Exp $ */ 2 /* ==================================================================== 3 * Copyright 2002 Sun Microsystems, Inc. ALL RIGHTS RESERVED. 4 * 5 * The Elliptic Curve Public-Key Crypto Library (ECC Code) included 6 * herein is developed by SUN MICROSYSTEMS, INC., and is contributed 7 * to the OpenSSL project. 8 * 9 * The ECC Code is licensed pursuant to the OpenSSL open source 10 * license provided below. 11 * 12 * In addition, Sun covenants to all licensees who provide a reciprocal 13 * covenant with respect to their own patents if any, not to sue under 14 * current and future patent claims necessarily infringed by the making, 15 * using, practicing, selling, offering for sale and/or otherwise 16 * disposing of the ECC Code as delivered hereunder (or portions thereof), 17 * provided that such covenant shall not apply: 18 * 1) for code that a licensee deletes from the ECC Code; 19 * 2) separates from the ECC Code; or 20 * 3) for infringements caused by: 21 * i) the modification of the ECC Code or 22 * ii) the combination of the ECC Code with other software or 23 * devices where such combination causes the infringement. 24 * 25 * The software is originally written by Sheueling Chang Shantz and 26 * Douglas Stebila of Sun Microsystems Laboratories. 27 * 28 */ 29 30 /* NOTE: This file is licensed pursuant to the OpenSSL license below 31 * and may be modified; but after modifications, the above covenant 32 * may no longer apply! In such cases, the corresponding paragraph 33 * ["In addition, Sun covenants ... causes the infringement."] and 34 * this note can be edited out; but please keep the Sun copyright 35 * notice and attribution. */ 36 37 /* ==================================================================== 38 * Copyright (c) 1998-2002 The OpenSSL Project. All rights reserved. 39 * 40 * Redistribution and use in source and binary forms, with or without 41 * modification, are permitted provided that the following conditions 42 * are met: 43 * 44 * 1. Redistributions of source code must retain the above copyright 45 * notice, this list of conditions and the following disclaimer. 46 * 47 * 2. Redistributions in binary form must reproduce the above copyright 48 * notice, this list of conditions and the following disclaimer in 49 * the documentation and/or other materials provided with the 50 * distribution. 51 * 52 * 3. All advertising materials mentioning features or use of this 53 * software must display the following acknowledgment: 54 * "This product includes software developed by the OpenSSL Project 55 * for use in the OpenSSL Toolkit. (http://www.openssl.org/)" 56 * 57 * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to 58 * endorse or promote products derived from this software without 59 * prior written permission. For written permission, please contact 60 * openssl-core@openssl.org. 61 * 62 * 5. Products derived from this software may not be called "OpenSSL" 63 * nor may "OpenSSL" appear in their names without prior written 64 * permission of the OpenSSL Project. 65 * 66 * 6. Redistributions of any form whatsoever must retain the following 67 * acknowledgment: 68 * "This product includes software developed by the OpenSSL Project 69 * for use in the OpenSSL Toolkit (http://www.openssl.org/)" 70 * 71 * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY 72 * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 73 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 74 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE OpenSSL PROJECT OR 75 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 76 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 77 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 78 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 79 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 80 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 81 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED 82 * OF THE POSSIBILITY OF SUCH DAMAGE. 83 * ==================================================================== 84 * 85 * This product includes cryptographic software written by Eric Young 86 * (eay@cryptsoft.com). This product includes software written by Tim 87 * Hudson (tjh@cryptsoft.com). 88 * 89 */ 90 91 #include <limits.h> 92 #include <stdio.h> 93 94 #include <openssl/opensslconf.h> 95 96 #include <openssl/err.h> 97 98 #include "bn_lcl.h" 99 100 #ifndef OPENSSL_NO_EC2M 101 102 /* Maximum number of iterations before BN_GF2m_mod_solve_quad_arr should fail. */ 103 #define MAX_ITERATIONS 50 104 105 static const BN_ULONG SQR_tb[16] = 106 { 0, 1, 4, 5, 16, 17, 20, 21, 107 64, 65, 68, 69, 80, 81, 84, 85 }; 108 /* Platform-specific macros to accelerate squaring. */ 109 #ifdef _LP64 110 #define SQR1(w) \ 111 SQR_tb[(w) >> 60 & 0xF] << 56 | SQR_tb[(w) >> 56 & 0xF] << 48 | \ 112 SQR_tb[(w) >> 52 & 0xF] << 40 | SQR_tb[(w) >> 48 & 0xF] << 32 | \ 113 SQR_tb[(w) >> 44 & 0xF] << 24 | SQR_tb[(w) >> 40 & 0xF] << 16 | \ 114 SQR_tb[(w) >> 36 & 0xF] << 8 | SQR_tb[(w) >> 32 & 0xF] 115 #define SQR0(w) \ 116 SQR_tb[(w) >> 28 & 0xF] << 56 | SQR_tb[(w) >> 24 & 0xF] << 48 | \ 117 SQR_tb[(w) >> 20 & 0xF] << 40 | SQR_tb[(w) >> 16 & 0xF] << 32 | \ 118 SQR_tb[(w) >> 12 & 0xF] << 24 | SQR_tb[(w) >> 8 & 0xF] << 16 | \ 119 SQR_tb[(w) >> 4 & 0xF] << 8 | SQR_tb[(w) & 0xF] 120 #else 121 #define SQR1(w) \ 122 SQR_tb[(w) >> 28 & 0xF] << 24 | SQR_tb[(w) >> 24 & 0xF] << 16 | \ 123 SQR_tb[(w) >> 20 & 0xF] << 8 | SQR_tb[(w) >> 16 & 0xF] 124 #define SQR0(w) \ 125 SQR_tb[(w) >> 12 & 0xF] << 24 | SQR_tb[(w) >> 8 & 0xF] << 16 | \ 126 SQR_tb[(w) >> 4 & 0xF] << 8 | SQR_tb[(w) & 0xF] 127 #endif 128 129 #if !defined(OPENSSL_BN_ASM_GF2m) 130 /* Product of two polynomials a, b each with degree < BN_BITS2 - 1, 131 * result is a polynomial r with degree < 2 * BN_BITS - 1 132 * The caller MUST ensure that the variables have the right amount 133 * of space allocated. 134 */ 135 static void 136 bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a, const BN_ULONG b) 137 { 138 #ifndef _LP64 139 BN_ULONG h, l, s; 140 BN_ULONG tab[8], top2b = a >> 30; 141 BN_ULONG a1, a2, a4; 142 143 a1 = a & (0x3FFFFFFF); 144 a2 = a1 << 1; 145 a4 = a2 << 1; 146 147 tab[0] = 0; 148 tab[1] = a1; 149 tab[2] = a2; 150 tab[3] = a1 ^ a2; 151 tab[4] = a4; 152 tab[5] = a1 ^ a4; 153 tab[6] = a2 ^ a4; 154 tab[7] = a1 ^ a2 ^ a4; 155 156 s = tab[b & 0x7]; 157 l = s; 158 s = tab[b >> 3 & 0x7]; 159 l ^= s << 3; 160 h = s >> 29; 161 s = tab[b >> 6 & 0x7]; 162 l ^= s << 6; 163 h ^= s >> 26; 164 s = tab[b >> 9 & 0x7]; 165 l ^= s << 9; 166 h ^= s >> 23; 167 s = tab[b >> 12 & 0x7]; 168 l ^= s << 12; 169 h ^= s >> 20; 170 s = tab[b >> 15 & 0x7]; 171 l ^= s << 15; 172 h ^= s >> 17; 173 s = tab[b >> 18 & 0x7]; 174 l ^= s << 18; 175 h ^= s >> 14; 176 s = tab[b >> 21 & 0x7]; 177 l ^= s << 21; 178 h ^= s >> 11; 179 s = tab[b >> 24 & 0x7]; 180 l ^= s << 24; 181 h ^= s >> 8; 182 s = tab[b >> 27 & 0x7]; 183 l ^= s << 27; 184 h ^= s >> 5; 185 s = tab[b >> 30]; 186 l ^= s << 30; 187 h ^= s >> 2; 188 189 /* compensate for the top two bits of a */ 190 if (top2b & 01) { 191 l ^= b << 30; 192 h ^= b >> 2; 193 } 194 if (top2b & 02) { 195 l ^= b << 31; 196 h ^= b >> 1; 197 } 198 199 *r1 = h; 200 *r0 = l; 201 #else 202 BN_ULONG h, l, s; 203 BN_ULONG tab[16], top3b = a >> 61; 204 BN_ULONG a1, a2, a4, a8; 205 206 a1 = a & (0x1FFFFFFFFFFFFFFFULL); 207 a2 = a1 << 1; 208 a4 = a2 << 1; 209 a8 = a4 << 1; 210 211 tab[0] = 0; 212 tab[1] = a1; 213 tab[2] = a2; 214 tab[3] = a1 ^ a2; 215 tab[4] = a4; 216 tab[5] = a1 ^ a4; 217 tab[6] = a2 ^ a4; 218 tab[7] = a1 ^ a2 ^ a4; 219 tab[8] = a8; 220 tab[9] = a1 ^ a8; 221 tab[10] = a2 ^ a8; 222 tab[11] = a1 ^ a2 ^ a8; 223 tab[12] = a4 ^ a8; 224 tab[13] = a1 ^ a4 ^ a8; 225 tab[14] = a2 ^ a4 ^ a8; 226 tab[15] = a1 ^ a2 ^ a4 ^ a8; 227 228 s = tab[b & 0xF]; 229 l = s; 230 s = tab[b >> 4 & 0xF]; 231 l ^= s << 4; 232 h = s >> 60; 233 s = tab[b >> 8 & 0xF]; 234 l ^= s << 8; 235 h ^= s >> 56; 236 s = tab[b >> 12 & 0xF]; 237 l ^= s << 12; 238 h ^= s >> 52; 239 s = tab[b >> 16 & 0xF]; 240 l ^= s << 16; 241 h ^= s >> 48; 242 s = tab[b >> 20 & 0xF]; 243 l ^= s << 20; 244 h ^= s >> 44; 245 s = tab[b >> 24 & 0xF]; 246 l ^= s << 24; 247 h ^= s >> 40; 248 s = tab[b >> 28 & 0xF]; 249 l ^= s << 28; 250 h ^= s >> 36; 251 s = tab[b >> 32 & 0xF]; 252 l ^= s << 32; 253 h ^= s >> 32; 254 s = tab[b >> 36 & 0xF]; 255 l ^= s << 36; 256 h ^= s >> 28; 257 s = tab[b >> 40 & 0xF]; 258 l ^= s << 40; 259 h ^= s >> 24; 260 s = tab[b >> 44 & 0xF]; 261 l ^= s << 44; 262 h ^= s >> 20; 263 s = tab[b >> 48 & 0xF]; 264 l ^= s << 48; 265 h ^= s >> 16; 266 s = tab[b >> 52 & 0xF]; 267 l ^= s << 52; 268 h ^= s >> 12; 269 s = tab[b >> 56 & 0xF]; 270 l ^= s << 56; 271 h ^= s >> 8; 272 s = tab[b >> 60]; 273 l ^= s << 60; 274 h ^= s >> 4; 275 276 /* compensate for the top three bits of a */ 277 if (top3b & 01) { 278 l ^= b << 61; 279 h ^= b >> 3; 280 } 281 if (top3b & 02) { 282 l ^= b << 62; 283 h ^= b >> 2; 284 } 285 if (top3b & 04) { 286 l ^= b << 63; 287 h ^= b >> 1; 288 } 289 290 *r1 = h; 291 *r0 = l; 292 #endif 293 } 294 295 /* Product of two polynomials a, b each with degree < 2 * BN_BITS2 - 1, 296 * result is a polynomial r with degree < 4 * BN_BITS2 - 1 297 * The caller MUST ensure that the variables have the right amount 298 * of space allocated. 299 */ 300 static void 301 bn_GF2m_mul_2x2(BN_ULONG *r, const BN_ULONG a1, const BN_ULONG a0, 302 const BN_ULONG b1, const BN_ULONG b0) 303 { 304 BN_ULONG m1, m0; 305 306 /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */ 307 bn_GF2m_mul_1x1(r + 3, r + 2, a1, b1); 308 bn_GF2m_mul_1x1(r + 1, r, a0, b0); 309 bn_GF2m_mul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1); 310 /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */ 311 r[2] ^= m1 ^ r[1] ^ r[3]; /* h0 ^= m1 ^ l1 ^ h1; */ 312 r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */ 313 } 314 #else 315 void bn_GF2m_mul_2x2(BN_ULONG *r, BN_ULONG a1, BN_ULONG a0, BN_ULONG b1, 316 BN_ULONG b0); 317 #endif 318 319 /* Add polynomials a and b and store result in r; r could be a or b, a and b 320 * could be equal; r is the bitwise XOR of a and b. 321 */ 322 int 323 BN_GF2m_add(BIGNUM *r, const BIGNUM *a, const BIGNUM *b) 324 { 325 int i; 326 const BIGNUM *at, *bt; 327 328 bn_check_top(a); 329 bn_check_top(b); 330 331 if (a->top < b->top) { 332 at = b; 333 bt = a; 334 } else { 335 at = a; 336 bt = b; 337 } 338 339 if (bn_wexpand(r, at->top) == NULL) 340 return 0; 341 342 for (i = 0; i < bt->top; i++) { 343 r->d[i] = at->d[i] ^ bt->d[i]; 344 } 345 for (; i < at->top; i++) { 346 r->d[i] = at->d[i]; 347 } 348 349 r->top = at->top; 350 bn_correct_top(r); 351 352 return 1; 353 } 354 355 356 /* Some functions allow for representation of the irreducible polynomials 357 * as an int[], say p. The irreducible f(t) is then of the form: 358 * t^p[0] + t^p[1] + ... + t^p[k] 359 * where m = p[0] > p[1] > ... > p[k] = 0. 360 */ 361 362 363 /* Performs modular reduction of a and store result in r. r could be a. */ 364 int 365 BN_GF2m_mod_arr(BIGNUM *r, const BIGNUM *a, const int p[]) 366 { 367 int j, k; 368 int n, dN, d0, d1; 369 BN_ULONG zz, *z; 370 371 bn_check_top(a); 372 373 if (!p[0]) { 374 /* reduction mod 1 => return 0 */ 375 BN_zero(r); 376 return 1; 377 } 378 379 /* Since the algorithm does reduction in the r value, if a != r, copy 380 * the contents of a into r so we can do reduction in r. 381 */ 382 if (a != r) { 383 if (!bn_wexpand(r, a->top)) 384 return 0; 385 for (j = 0; j < a->top; j++) { 386 r->d[j] = a->d[j]; 387 } 388 r->top = a->top; 389 } 390 z = r->d; 391 392 /* start reduction */ 393 dN = p[0] / BN_BITS2; 394 for (j = r->top - 1; j > dN; ) { 395 zz = z[j]; 396 if (z[j] == 0) { 397 j--; 398 continue; 399 } 400 z[j] = 0; 401 402 for (k = 1; p[k] != 0; k++) { 403 /* reducing component t^p[k] */ 404 n = p[0] - p[k]; 405 d0 = n % BN_BITS2; 406 d1 = BN_BITS2 - d0; 407 n /= BN_BITS2; 408 z[j - n] ^= (zz >> d0); 409 if (d0) 410 z[j - n - 1] ^= (zz << d1); 411 } 412 413 /* reducing component t^0 */ 414 n = dN; 415 d0 = p[0] % BN_BITS2; 416 d1 = BN_BITS2 - d0; 417 z[j - n] ^= (zz >> d0); 418 if (d0) 419 z[j - n - 1] ^= (zz << d1); 420 } 421 422 /* final round of reduction */ 423 while (j == dN) { 424 425 d0 = p[0] % BN_BITS2; 426 zz = z[dN] >> d0; 427 if (zz == 0) 428 break; 429 d1 = BN_BITS2 - d0; 430 431 /* clear up the top d1 bits */ 432 if (d0) 433 z[dN] = (z[dN] << d1) >> d1; 434 else 435 z[dN] = 0; 436 z[0] ^= zz; /* reduction t^0 component */ 437 438 for (k = 1; p[k] != 0; k++) { 439 BN_ULONG tmp_ulong; 440 441 /* reducing component t^p[k]*/ 442 n = p[k] / BN_BITS2; 443 d0 = p[k] % BN_BITS2; 444 d1 = BN_BITS2 - d0; 445 z[n] ^= (zz << d0); 446 if (d0 && (tmp_ulong = zz >> d1)) 447 z[n + 1] ^= tmp_ulong; 448 } 449 450 451 } 452 453 bn_correct_top(r); 454 return 1; 455 } 456 457 /* Performs modular reduction of a by p and store result in r. r could be a. 458 * 459 * This function calls down to the BN_GF2m_mod_arr implementation; this wrapper 460 * function is only provided for convenience; for best performance, use the 461 * BN_GF2m_mod_arr function. 462 */ 463 int 464 BN_GF2m_mod(BIGNUM *r, const BIGNUM *a, const BIGNUM *p) 465 { 466 int ret = 0; 467 int arr[6]; 468 469 bn_check_top(a); 470 bn_check_top(p); 471 ret = BN_GF2m_poly2arr(p, arr, sizeof(arr) / sizeof(arr[0])); 472 if (!ret || ret > (int)(sizeof(arr) / sizeof(arr[0]))) { 473 BNerror(BN_R_INVALID_LENGTH); 474 return 0; 475 } 476 ret = BN_GF2m_mod_arr(r, a, arr); 477 bn_check_top(r); 478 return ret; 479 } 480 481 482 /* Compute the product of two polynomials a and b, reduce modulo p, and store 483 * the result in r. r could be a or b; a could be b. 484 */ 485 int 486 BN_GF2m_mod_mul_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, const int p[], 487 BN_CTX *ctx) 488 { 489 int zlen, i, j, k, ret = 0; 490 BIGNUM *s; 491 BN_ULONG x1, x0, y1, y0, zz[4]; 492 493 bn_check_top(a); 494 bn_check_top(b); 495 496 if (a == b) { 497 return BN_GF2m_mod_sqr_arr(r, a, p, ctx); 498 } 499 500 BN_CTX_start(ctx); 501 if ((s = BN_CTX_get(ctx)) == NULL) 502 goto err; 503 504 zlen = a->top + b->top + 4; 505 if (!bn_wexpand(s, zlen)) 506 goto err; 507 s->top = zlen; 508 509 for (i = 0; i < zlen; i++) 510 s->d[i] = 0; 511 512 for (j = 0; j < b->top; j += 2) { 513 y0 = b->d[j]; 514 y1 = ((j + 1) == b->top) ? 0 : b->d[j + 1]; 515 for (i = 0; i < a->top; i += 2) { 516 x0 = a->d[i]; 517 x1 = ((i + 1) == a->top) ? 0 : a->d[i + 1]; 518 bn_GF2m_mul_2x2(zz, x1, x0, y1, y0); 519 for (k = 0; k < 4; k++) 520 s->d[i + j + k] ^= zz[k]; 521 } 522 } 523 524 bn_correct_top(s); 525 if (BN_GF2m_mod_arr(r, s, p)) 526 ret = 1; 527 bn_check_top(r); 528 529 err: 530 BN_CTX_end(ctx); 531 return ret; 532 } 533 534 /* Compute the product of two polynomials a and b, reduce modulo p, and store 535 * the result in r. r could be a or b; a could equal b. 536 * 537 * This function calls down to the BN_GF2m_mod_mul_arr implementation; this wrapper 538 * function is only provided for convenience; for best performance, use the 539 * BN_GF2m_mod_mul_arr function. 540 */ 541 int 542 BN_GF2m_mod_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, const BIGNUM *p, 543 BN_CTX *ctx) 544 { 545 int ret = 0; 546 const int max = BN_num_bits(p) + 1; 547 int *arr = NULL; 548 549 bn_check_top(a); 550 bn_check_top(b); 551 bn_check_top(p); 552 if ((arr = reallocarray(NULL, max, sizeof(int))) == NULL) 553 goto err; 554 ret = BN_GF2m_poly2arr(p, arr, max); 555 if (!ret || ret > max) { 556 BNerror(BN_R_INVALID_LENGTH); 557 goto err; 558 } 559 ret = BN_GF2m_mod_mul_arr(r, a, b, arr, ctx); 560 bn_check_top(r); 561 562 err: 563 free(arr); 564 return ret; 565 } 566 567 568 /* Square a, reduce the result mod p, and store it in a. r could be a. */ 569 int 570 BN_GF2m_mod_sqr_arr(BIGNUM *r, const BIGNUM *a, const int p[], BN_CTX *ctx) 571 { 572 int i, ret = 0; 573 BIGNUM *s; 574 575 bn_check_top(a); 576 BN_CTX_start(ctx); 577 if ((s = BN_CTX_get(ctx)) == NULL) 578 goto err; 579 if (!bn_wexpand(s, 2 * a->top)) 580 goto err; 581 582 for (i = a->top - 1; i >= 0; i--) { 583 s->d[2 * i + 1] = SQR1(a->d[i]); 584 s->d[2 * i] = SQR0(a->d[i]); 585 } 586 587 s->top = 2 * a->top; 588 bn_correct_top(s); 589 if (!BN_GF2m_mod_arr(r, s, p)) 590 goto err; 591 bn_check_top(r); 592 ret = 1; 593 594 err: 595 BN_CTX_end(ctx); 596 return ret; 597 } 598 599 /* Square a, reduce the result mod p, and store it in a. r could be a. 600 * 601 * This function calls down to the BN_GF2m_mod_sqr_arr implementation; this wrapper 602 * function is only provided for convenience; for best performance, use the 603 * BN_GF2m_mod_sqr_arr function. 604 */ 605 int 606 BN_GF2m_mod_sqr(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 607 { 608 int ret = 0; 609 const int max = BN_num_bits(p) + 1; 610 int *arr = NULL; 611 612 bn_check_top(a); 613 bn_check_top(p); 614 if ((arr = reallocarray(NULL, max, sizeof(int))) == NULL) 615 goto err; 616 ret = BN_GF2m_poly2arr(p, arr, max); 617 if (!ret || ret > max) { 618 BNerror(BN_R_INVALID_LENGTH); 619 goto err; 620 } 621 ret = BN_GF2m_mod_sqr_arr(r, a, arr, ctx); 622 bn_check_top(r); 623 624 err: 625 free(arr); 626 return ret; 627 } 628 629 630 /* Invert a, reduce modulo p, and store the result in r. r could be a. 631 * Uses Modified Almost Inverse Algorithm (Algorithm 10) from 632 * Hankerson, D., Hernandez, J.L., and Menezes, A. "Software Implementation 633 * of Elliptic Curve Cryptography Over Binary Fields". 634 */ 635 int 636 BN_GF2m_mod_inv(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 637 { 638 BIGNUM *b, *c = NULL, *u = NULL, *v = NULL, *tmp; 639 int ret = 0; 640 641 bn_check_top(a); 642 bn_check_top(p); 643 644 BN_CTX_start(ctx); 645 646 if ((b = BN_CTX_get(ctx)) == NULL) 647 goto err; 648 if ((c = BN_CTX_get(ctx)) == NULL) 649 goto err; 650 if ((u = BN_CTX_get(ctx)) == NULL) 651 goto err; 652 if ((v = BN_CTX_get(ctx)) == NULL) 653 goto err; 654 655 if (!BN_GF2m_mod(u, a, p)) 656 goto err; 657 if (BN_is_zero(u)) 658 goto err; 659 660 if (!BN_copy(v, p)) 661 goto err; 662 #if 0 663 if (!BN_one(b)) 664 goto err; 665 666 while (1) { 667 while (!BN_is_odd(u)) { 668 if (BN_is_zero(u)) 669 goto err; 670 if (!BN_rshift1(u, u)) 671 goto err; 672 if (BN_is_odd(b)) { 673 if (!BN_GF2m_add(b, b, p)) 674 goto err; 675 } 676 if (!BN_rshift1(b, b)) 677 goto err; 678 } 679 680 if (BN_abs_is_word(u, 1)) 681 break; 682 683 if (BN_num_bits(u) < BN_num_bits(v)) { 684 tmp = u; 685 u = v; 686 v = tmp; 687 tmp = b; 688 b = c; 689 c = tmp; 690 } 691 692 if (!BN_GF2m_add(u, u, v)) 693 goto err; 694 if (!BN_GF2m_add(b, b, c)) 695 goto err; 696 } 697 #else 698 { 699 int i, ubits = BN_num_bits(u), 700 vbits = BN_num_bits(v), /* v is copy of p */ 701 top = p->top; 702 BN_ULONG *udp, *bdp, *vdp, *cdp; 703 704 if (!bn_wexpand(u, top)) 705 goto err; 706 udp = u->d; 707 for (i = u->top; i < top; i++) 708 udp[i] = 0; 709 u->top = top; 710 if (!bn_wexpand(b, top)) 711 goto err; 712 bdp = b->d; 713 bdp[0] = 1; 714 for (i = 1; i < top; i++) 715 bdp[i] = 0; 716 b->top = top; 717 if (!bn_wexpand(c, top)) 718 goto err; 719 cdp = c->d; 720 for (i = 0; i < top; i++) 721 cdp[i] = 0; 722 c->top = top; 723 vdp = v->d; /* It pays off to "cache" *->d pointers, because 724 * it allows optimizer to be more aggressive. 725 * But we don't have to "cache" p->d, because *p 726 * is declared 'const'... */ 727 while (1) { 728 while (ubits && !(udp[0]&1)) { 729 BN_ULONG u0, u1, b0, b1, mask; 730 731 u0 = udp[0]; 732 b0 = bdp[0]; 733 mask = (BN_ULONG)0 - (b0 & 1); 734 b0 ^= p->d[0] & mask; 735 for (i = 0; i < top - 1; i++) { 736 u1 = udp[i + 1]; 737 udp[i] = ((u0 >> 1) | 738 (u1 << (BN_BITS2 - 1))) & BN_MASK2; 739 u0 = u1; 740 b1 = bdp[i + 1] ^ (p->d[i + 1] & mask); 741 bdp[i] = ((b0 >> 1) | 742 (b1 << (BN_BITS2 - 1))) & BN_MASK2; 743 b0 = b1; 744 } 745 udp[i] = u0 >> 1; 746 bdp[i] = b0 >> 1; 747 ubits--; 748 } 749 750 if (ubits <= BN_BITS2) { 751 /* See if poly was reducible. */ 752 if (udp[0] == 0) 753 goto err; 754 if (udp[0] == 1) 755 break; 756 } 757 758 if (ubits < vbits) { 759 i = ubits; 760 ubits = vbits; 761 vbits = i; 762 tmp = u; 763 u = v; 764 v = tmp; 765 tmp = b; 766 b = c; 767 c = tmp; 768 udp = vdp; 769 vdp = v->d; 770 bdp = cdp; 771 cdp = c->d; 772 } 773 for (i = 0; i < top; i++) { 774 udp[i] ^= vdp[i]; 775 bdp[i] ^= cdp[i]; 776 } 777 if (ubits == vbits) { 778 BN_ULONG ul; 779 int utop = (ubits - 1) / BN_BITS2; 780 781 while ((ul = udp[utop]) == 0 && utop) 782 utop--; 783 ubits = utop*BN_BITS2 + BN_num_bits_word(ul); 784 } 785 } 786 bn_correct_top(b); 787 } 788 #endif 789 790 if (!BN_copy(r, b)) 791 goto err; 792 bn_check_top(r); 793 ret = 1; 794 795 err: 796 #ifdef BN_DEBUG /* BN_CTX_end would complain about the expanded form */ 797 bn_correct_top(c); 798 bn_correct_top(u); 799 bn_correct_top(v); 800 #endif 801 BN_CTX_end(ctx); 802 return ret; 803 } 804 805 /* Invert xx, reduce modulo p, and store the result in r. r could be xx. 806 * 807 * This function calls down to the BN_GF2m_mod_inv implementation; this wrapper 808 * function is only provided for convenience; for best performance, use the 809 * BN_GF2m_mod_inv function. 810 */ 811 int 812 BN_GF2m_mod_inv_arr(BIGNUM *r, const BIGNUM *xx, const int p[], BN_CTX *ctx) 813 { 814 BIGNUM *field; 815 int ret = 0; 816 817 bn_check_top(xx); 818 BN_CTX_start(ctx); 819 if ((field = BN_CTX_get(ctx)) == NULL) 820 goto err; 821 if (!BN_GF2m_arr2poly(p, field)) 822 goto err; 823 824 ret = BN_GF2m_mod_inv(r, xx, field, ctx); 825 bn_check_top(r); 826 827 err: 828 BN_CTX_end(ctx); 829 return ret; 830 } 831 832 833 #ifndef OPENSSL_SUN_GF2M_DIV 834 /* Divide y by x, reduce modulo p, and store the result in r. r could be x 835 * or y, x could equal y. 836 */ 837 int 838 BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x, const BIGNUM *p, 839 BN_CTX *ctx) 840 { 841 BIGNUM *xinv = NULL; 842 int ret = 0; 843 844 bn_check_top(y); 845 bn_check_top(x); 846 bn_check_top(p); 847 848 BN_CTX_start(ctx); 849 if ((xinv = BN_CTX_get(ctx)) == NULL) 850 goto err; 851 852 if (!BN_GF2m_mod_inv(xinv, x, p, ctx)) 853 goto err; 854 if (!BN_GF2m_mod_mul(r, y, xinv, p, ctx)) 855 goto err; 856 bn_check_top(r); 857 ret = 1; 858 859 err: 860 BN_CTX_end(ctx); 861 return ret; 862 } 863 #else 864 /* Divide y by x, reduce modulo p, and store the result in r. r could be x 865 * or y, x could equal y. 866 * Uses algorithm Modular_Division_GF(2^m) from 867 * Chang-Shantz, S. "From Euclid's GCD to Montgomery Multiplication to 868 * the Great Divide". 869 */ 870 int 871 BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x, const BIGNUM *p, 872 BN_CTX *ctx) 873 { 874 BIGNUM *a, *b, *u, *v; 875 int ret = 0; 876 877 bn_check_top(y); 878 bn_check_top(x); 879 bn_check_top(p); 880 881 BN_CTX_start(ctx); 882 883 if ((a = BN_CTX_get(ctx)) == NULL) 884 goto err; 885 if ((b = BN_CTX_get(ctx)) == NULL) 886 goto err; 887 if ((u = BN_CTX_get(ctx)) == NULL) 888 goto err; 889 if ((v = BN_CTX_get(ctx)) == NULL) 890 goto err; 891 892 /* reduce x and y mod p */ 893 if (!BN_GF2m_mod(u, y, p)) 894 goto err; 895 if (!BN_GF2m_mod(a, x, p)) 896 goto err; 897 if (!BN_copy(b, p)) 898 goto err; 899 900 while (!BN_is_odd(a)) { 901 if (!BN_rshift1(a, a)) 902 goto err; 903 if (BN_is_odd(u)) 904 if (!BN_GF2m_add(u, u, p)) 905 goto err; 906 if (!BN_rshift1(u, u)) 907 goto err; 908 } 909 910 do { 911 if (BN_GF2m_cmp(b, a) > 0) { 912 if (!BN_GF2m_add(b, b, a)) 913 goto err; 914 if (!BN_GF2m_add(v, v, u)) 915 goto err; 916 do { 917 if (!BN_rshift1(b, b)) 918 goto err; 919 if (BN_is_odd(v)) 920 if (!BN_GF2m_add(v, v, p)) 921 goto err; 922 if (!BN_rshift1(v, v)) 923 goto err; 924 } while (!BN_is_odd(b)); 925 } else if (BN_abs_is_word(a, 1)) 926 break; 927 else { 928 if (!BN_GF2m_add(a, a, b)) 929 goto err; 930 if (!BN_GF2m_add(u, u, v)) 931 goto err; 932 do { 933 if (!BN_rshift1(a, a)) 934 goto err; 935 if (BN_is_odd(u)) 936 if (!BN_GF2m_add(u, u, p)) 937 goto err; 938 if (!BN_rshift1(u, u)) 939 goto err; 940 } while (!BN_is_odd(a)); 941 } 942 } while (1); 943 944 if (!BN_copy(r, u)) 945 goto err; 946 bn_check_top(r); 947 ret = 1; 948 949 err: 950 BN_CTX_end(ctx); 951 return ret; 952 } 953 #endif 954 955 /* Divide yy by xx, reduce modulo p, and store the result in r. r could be xx 956 * or yy, xx could equal yy. 957 * 958 * This function calls down to the BN_GF2m_mod_div implementation; this wrapper 959 * function is only provided for convenience; for best performance, use the 960 * BN_GF2m_mod_div function. 961 */ 962 int 963 BN_GF2m_mod_div_arr(BIGNUM *r, const BIGNUM *yy, const BIGNUM *xx, 964 const int p[], BN_CTX *ctx) 965 { 966 BIGNUM *field; 967 int ret = 0; 968 969 bn_check_top(yy); 970 bn_check_top(xx); 971 972 BN_CTX_start(ctx); 973 if ((field = BN_CTX_get(ctx)) == NULL) 974 goto err; 975 if (!BN_GF2m_arr2poly(p, field)) 976 goto err; 977 978 ret = BN_GF2m_mod_div(r, yy, xx, field, ctx); 979 bn_check_top(r); 980 981 err: 982 BN_CTX_end(ctx); 983 return ret; 984 } 985 986 987 /* Compute the bth power of a, reduce modulo p, and store 988 * the result in r. r could be a. 989 * Uses simple square-and-multiply algorithm A.5.1 from IEEE P1363. 990 */ 991 int 992 BN_GF2m_mod_exp_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, const int p[], 993 BN_CTX *ctx) 994 { 995 int ret = 0, i, n; 996 BIGNUM *u; 997 998 bn_check_top(a); 999 bn_check_top(b); 1000 1001 if (BN_is_zero(b)) 1002 return (BN_one(r)); 1003 1004 if (BN_abs_is_word(b, 1)) 1005 return (BN_copy(r, a) != NULL); 1006 1007 BN_CTX_start(ctx); 1008 if ((u = BN_CTX_get(ctx)) == NULL) 1009 goto err; 1010 1011 if (!BN_GF2m_mod_arr(u, a, p)) 1012 goto err; 1013 1014 n = BN_num_bits(b) - 1; 1015 for (i = n - 1; i >= 0; i--) { 1016 if (!BN_GF2m_mod_sqr_arr(u, u, p, ctx)) 1017 goto err; 1018 if (BN_is_bit_set(b, i)) { 1019 if (!BN_GF2m_mod_mul_arr(u, u, a, p, ctx)) 1020 goto err; 1021 } 1022 } 1023 if (!BN_copy(r, u)) 1024 goto err; 1025 bn_check_top(r); 1026 ret = 1; 1027 1028 err: 1029 BN_CTX_end(ctx); 1030 return ret; 1031 } 1032 1033 /* Compute the bth power of a, reduce modulo p, and store 1034 * the result in r. r could be a. 1035 * 1036 * This function calls down to the BN_GF2m_mod_exp_arr implementation; this wrapper 1037 * function is only provided for convenience; for best performance, use the 1038 * BN_GF2m_mod_exp_arr function. 1039 */ 1040 int 1041 BN_GF2m_mod_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, const BIGNUM *p, 1042 BN_CTX *ctx) 1043 { 1044 int ret = 0; 1045 const int max = BN_num_bits(p) + 1; 1046 int *arr = NULL; 1047 1048 bn_check_top(a); 1049 bn_check_top(b); 1050 bn_check_top(p); 1051 if ((arr = reallocarray(NULL, max, sizeof(int))) == NULL) 1052 goto err; 1053 ret = BN_GF2m_poly2arr(p, arr, max); 1054 if (!ret || ret > max) { 1055 BNerror(BN_R_INVALID_LENGTH); 1056 goto err; 1057 } 1058 ret = BN_GF2m_mod_exp_arr(r, a, b, arr, ctx); 1059 bn_check_top(r); 1060 1061 err: 1062 free(arr); 1063 return ret; 1064 } 1065 1066 /* Compute the square root of a, reduce modulo p, and store 1067 * the result in r. r could be a. 1068 * Uses exponentiation as in algorithm A.4.1 from IEEE P1363. 1069 */ 1070 int 1071 BN_GF2m_mod_sqrt_arr(BIGNUM *r, const BIGNUM *a, const int p[], BN_CTX *ctx) 1072 { 1073 int ret = 0; 1074 BIGNUM *u; 1075 1076 bn_check_top(a); 1077 1078 if (!p[0]) { 1079 /* reduction mod 1 => return 0 */ 1080 BN_zero(r); 1081 return 1; 1082 } 1083 1084 BN_CTX_start(ctx); 1085 if ((u = BN_CTX_get(ctx)) == NULL) 1086 goto err; 1087 1088 if (!BN_set_bit(u, p[0] - 1)) 1089 goto err; 1090 ret = BN_GF2m_mod_exp_arr(r, a, u, p, ctx); 1091 bn_check_top(r); 1092 1093 err: 1094 BN_CTX_end(ctx); 1095 return ret; 1096 } 1097 1098 /* Compute the square root of a, reduce modulo p, and store 1099 * the result in r. r could be a. 1100 * 1101 * This function calls down to the BN_GF2m_mod_sqrt_arr implementation; this wrapper 1102 * function is only provided for convenience; for best performance, use the 1103 * BN_GF2m_mod_sqrt_arr function. 1104 */ 1105 int 1106 BN_GF2m_mod_sqrt(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 1107 { 1108 int ret = 0; 1109 const int max = BN_num_bits(p) + 1; 1110 int *arr = NULL; 1111 bn_check_top(a); 1112 bn_check_top(p); 1113 if ((arr = reallocarray(NULL, max, sizeof(int))) == NULL) 1114 goto err; 1115 ret = BN_GF2m_poly2arr(p, arr, max); 1116 if (!ret || ret > max) { 1117 BNerror(BN_R_INVALID_LENGTH); 1118 goto err; 1119 } 1120 ret = BN_GF2m_mod_sqrt_arr(r, a, arr, ctx); 1121 bn_check_top(r); 1122 1123 err: 1124 free(arr); 1125 return ret; 1126 } 1127 1128 /* Find r such that r^2 + r = a mod p. r could be a. If no r exists returns 0. 1129 * Uses algorithms A.4.7 and A.4.6 from IEEE P1363. 1130 */ 1131 int 1132 BN_GF2m_mod_solve_quad_arr(BIGNUM *r, const BIGNUM *a_, const int p[], 1133 BN_CTX *ctx) 1134 { 1135 int ret = 0, count = 0, j; 1136 BIGNUM *a, *z, *rho, *w, *w2, *tmp; 1137 1138 bn_check_top(a_); 1139 1140 if (!p[0]) { 1141 /* reduction mod 1 => return 0 */ 1142 BN_zero(r); 1143 return 1; 1144 } 1145 1146 BN_CTX_start(ctx); 1147 if ((a = BN_CTX_get(ctx)) == NULL) 1148 goto err; 1149 if ((z = BN_CTX_get(ctx)) == NULL) 1150 goto err; 1151 if ((w = BN_CTX_get(ctx)) == NULL) 1152 goto err; 1153 1154 if (!BN_GF2m_mod_arr(a, a_, p)) 1155 goto err; 1156 1157 if (BN_is_zero(a)) { 1158 BN_zero(r); 1159 ret = 1; 1160 goto err; 1161 } 1162 1163 if (p[0] & 0x1) /* m is odd */ 1164 { 1165 /* compute half-trace of a */ 1166 if (!BN_copy(z, a)) 1167 goto err; 1168 for (j = 1; j <= (p[0] - 1) / 2; j++) { 1169 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx)) 1170 goto err; 1171 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx)) 1172 goto err; 1173 if (!BN_GF2m_add(z, z, a)) 1174 goto err; 1175 } 1176 1177 } 1178 else /* m is even */ 1179 { 1180 if ((rho = BN_CTX_get(ctx)) == NULL) 1181 goto err; 1182 if ((w2 = BN_CTX_get(ctx)) == NULL) 1183 goto err; 1184 if ((tmp = BN_CTX_get(ctx)) == NULL) 1185 goto err; 1186 do { 1187 if (!BN_rand(rho, p[0], 0, 0)) 1188 goto err; 1189 if (!BN_GF2m_mod_arr(rho, rho, p)) 1190 goto err; 1191 BN_zero(z); 1192 if (!BN_copy(w, rho)) 1193 goto err; 1194 for (j = 1; j <= p[0] - 1; j++) { 1195 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx)) 1196 goto err; 1197 if (!BN_GF2m_mod_sqr_arr(w2, w, p, ctx)) 1198 goto err; 1199 if (!BN_GF2m_mod_mul_arr(tmp, w2, a, p, ctx)) 1200 goto err; 1201 if (!BN_GF2m_add(z, z, tmp)) 1202 goto err; 1203 if (!BN_GF2m_add(w, w2, rho)) 1204 goto err; 1205 } 1206 count++; 1207 } while (BN_is_zero(w) && (count < MAX_ITERATIONS)); 1208 if (BN_is_zero(w)) { 1209 BNerror(BN_R_TOO_MANY_ITERATIONS); 1210 goto err; 1211 } 1212 } 1213 1214 if (!BN_GF2m_mod_sqr_arr(w, z, p, ctx)) 1215 goto err; 1216 if (!BN_GF2m_add(w, z, w)) 1217 goto err; 1218 if (BN_GF2m_cmp(w, a)) { 1219 BNerror(BN_R_NO_SOLUTION); 1220 goto err; 1221 } 1222 1223 if (!BN_copy(r, z)) 1224 goto err; 1225 bn_check_top(r); 1226 1227 ret = 1; 1228 1229 err: 1230 BN_CTX_end(ctx); 1231 return ret; 1232 } 1233 1234 /* Find r such that r^2 + r = a mod p. r could be a. If no r exists returns 0. 1235 * 1236 * This function calls down to the BN_GF2m_mod_solve_quad_arr implementation; this wrapper 1237 * function is only provided for convenience; for best performance, use the 1238 * BN_GF2m_mod_solve_quad_arr function. 1239 */ 1240 int 1241 BN_GF2m_mod_solve_quad(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 1242 { 1243 int ret = 0; 1244 const int max = BN_num_bits(p) + 1; 1245 int *arr = NULL; 1246 1247 bn_check_top(a); 1248 bn_check_top(p); 1249 if ((arr = reallocarray(NULL, max, sizeof(int))) == NULL) 1250 goto err; 1251 ret = BN_GF2m_poly2arr(p, arr, max); 1252 if (!ret || ret > max) { 1253 BNerror(BN_R_INVALID_LENGTH); 1254 goto err; 1255 } 1256 ret = BN_GF2m_mod_solve_quad_arr(r, a, arr, ctx); 1257 bn_check_top(r); 1258 1259 err: 1260 free(arr); 1261 return ret; 1262 } 1263 1264 /* Convert the bit-string representation of a polynomial 1265 * ( \sum_{i=0}^n a_i * x^i) into an array of integers corresponding 1266 * to the bits with non-zero coefficient. Array is terminated with -1. 1267 * Up to max elements of the array will be filled. Return value is total 1268 * number of array elements that would be filled if array was large enough. 1269 */ 1270 int 1271 BN_GF2m_poly2arr(const BIGNUM *a, int p[], int max) 1272 { 1273 int i, j, k = 0; 1274 BN_ULONG mask; 1275 1276 if (BN_is_zero(a)) 1277 return 0; 1278 1279 for (i = a->top - 1; i >= 0; i--) { 1280 if (!a->d[i]) 1281 /* skip word if a->d[i] == 0 */ 1282 continue; 1283 mask = BN_TBIT; 1284 for (j = BN_BITS2 - 1; j >= 0; j--) { 1285 if (a->d[i] & mask) { 1286 if (k < max) 1287 p[k] = BN_BITS2 * i + j; 1288 k++; 1289 } 1290 mask >>= 1; 1291 } 1292 } 1293 1294 if (k < max) { 1295 p[k] = -1; 1296 k++; 1297 } 1298 1299 return k; 1300 } 1301 1302 /* Convert the coefficient array representation of a polynomial to a 1303 * bit-string. The array must be terminated by -1. 1304 */ 1305 int 1306 BN_GF2m_arr2poly(const int p[], BIGNUM *a) 1307 { 1308 int i; 1309 1310 bn_check_top(a); 1311 BN_zero(a); 1312 for (i = 0; p[i] != -1; i++) { 1313 if (BN_set_bit(a, p[i]) == 0) 1314 return 0; 1315 } 1316 bn_check_top(a); 1317 1318 return 1; 1319 } 1320 1321 #endif 1322