mp_gf2m.c (16393B)
1 /* This Source Code Form is subject to the terms of the Mozilla Public 2 * License, v. 2.0. If a copy of the MPL was not distributed with this 3 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ 4 5 #include "mp_gf2m.h" 6 #include "mp_gf2m-priv.h" 7 #include "mplogic.h" 8 #include "mpi-priv.h" 9 10 const mp_digit mp_gf2m_sqr_tb[16] = { 11 0, 1, 4, 5, 16, 17, 20, 21, 12 64, 65, 68, 69, 80, 81, 84, 85 13 }; 14 15 /* Multiply two binary polynomials mp_digits a, b. 16 * Result is a polynomial with degree < 2 * MP_DIGIT_BITS - 1. 17 * Output in two mp_digits rh, rl. 18 */ 19 #if MP_DIGIT_BITS == 32 20 void 21 s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b) 22 { 23 register mp_digit h, l, s; 24 mp_digit tab[8], top2b = a >> 30; 25 register mp_digit a1, a2, a4; 26 27 a1 = a & (0x3FFFFFFF); 28 a2 = a1 << 1; 29 a4 = a2 << 1; 30 31 tab[0] = 0; 32 tab[1] = a1; 33 tab[2] = a2; 34 tab[3] = a1 ^ a2; 35 tab[4] = a4; 36 tab[5] = a1 ^ a4; 37 tab[6] = a2 ^ a4; 38 tab[7] = a1 ^ a2 ^ a4; 39 40 s = tab[b & 0x7]; 41 l = s; 42 s = tab[b >> 3 & 0x7]; 43 l ^= s << 3; 44 h = s >> 29; 45 s = tab[b >> 6 & 0x7]; 46 l ^= s << 6; 47 h ^= s >> 26; 48 s = tab[b >> 9 & 0x7]; 49 l ^= s << 9; 50 h ^= s >> 23; 51 s = tab[b >> 12 & 0x7]; 52 l ^= s << 12; 53 h ^= s >> 20; 54 s = tab[b >> 15 & 0x7]; 55 l ^= s << 15; 56 h ^= s >> 17; 57 s = tab[b >> 18 & 0x7]; 58 l ^= s << 18; 59 h ^= s >> 14; 60 s = tab[b >> 21 & 0x7]; 61 l ^= s << 21; 62 h ^= s >> 11; 63 s = tab[b >> 24 & 0x7]; 64 l ^= s << 24; 65 h ^= s >> 8; 66 s = tab[b >> 27 & 0x7]; 67 l ^= s << 27; 68 h ^= s >> 5; 69 s = tab[b >> 30]; 70 l ^= s << 30; 71 h ^= s >> 2; 72 73 /* compensate for the top two bits of a */ 74 75 if (top2b & 01) { 76 l ^= b << 30; 77 h ^= b >> 2; 78 } 79 if (top2b & 02) { 80 l ^= b << 31; 81 h ^= b >> 1; 82 } 83 84 *rh = h; 85 *rl = l; 86 } 87 #else 88 void 89 s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b) 90 { 91 register mp_digit h, l, s; 92 mp_digit tab[16], top3b = a >> 61; 93 register mp_digit a1, a2, a4, a8; 94 95 a1 = a & (0x1FFFFFFFFFFFFFFFULL); 96 a2 = a1 << 1; 97 a4 = a2 << 1; 98 a8 = a4 << 1; 99 tab[0] = 0; 100 tab[1] = a1; 101 tab[2] = a2; 102 tab[3] = a1 ^ a2; 103 tab[4] = a4; 104 tab[5] = a1 ^ a4; 105 tab[6] = a2 ^ a4; 106 tab[7] = a1 ^ a2 ^ a4; 107 tab[8] = a8; 108 tab[9] = a1 ^ a8; 109 tab[10] = a2 ^ a8; 110 tab[11] = a1 ^ a2 ^ a8; 111 tab[12] = a4 ^ a8; 112 tab[13] = a1 ^ a4 ^ a8; 113 tab[14] = a2 ^ a4 ^ a8; 114 tab[15] = a1 ^ a2 ^ a4 ^ a8; 115 116 s = tab[b & 0xF]; 117 l = s; 118 s = tab[b >> 4 & 0xF]; 119 l ^= s << 4; 120 h = s >> 60; 121 s = tab[b >> 8 & 0xF]; 122 l ^= s << 8; 123 h ^= s >> 56; 124 s = tab[b >> 12 & 0xF]; 125 l ^= s << 12; 126 h ^= s >> 52; 127 s = tab[b >> 16 & 0xF]; 128 l ^= s << 16; 129 h ^= s >> 48; 130 s = tab[b >> 20 & 0xF]; 131 l ^= s << 20; 132 h ^= s >> 44; 133 s = tab[b >> 24 & 0xF]; 134 l ^= s << 24; 135 h ^= s >> 40; 136 s = tab[b >> 28 & 0xF]; 137 l ^= s << 28; 138 h ^= s >> 36; 139 s = tab[b >> 32 & 0xF]; 140 l ^= s << 32; 141 h ^= s >> 32; 142 s = tab[b >> 36 & 0xF]; 143 l ^= s << 36; 144 h ^= s >> 28; 145 s = tab[b >> 40 & 0xF]; 146 l ^= s << 40; 147 h ^= s >> 24; 148 s = tab[b >> 44 & 0xF]; 149 l ^= s << 44; 150 h ^= s >> 20; 151 s = tab[b >> 48 & 0xF]; 152 l ^= s << 48; 153 h ^= s >> 16; 154 s = tab[b >> 52 & 0xF]; 155 l ^= s << 52; 156 h ^= s >> 12; 157 s = tab[b >> 56 & 0xF]; 158 l ^= s << 56; 159 h ^= s >> 8; 160 s = tab[b >> 60]; 161 l ^= s << 60; 162 h ^= s >> 4; 163 164 /* compensate for the top three bits of a */ 165 166 if (top3b & 01) { 167 l ^= b << 61; 168 h ^= b >> 3; 169 } 170 if (top3b & 02) { 171 l ^= b << 62; 172 h ^= b >> 2; 173 } 174 if (top3b & 04) { 175 l ^= b << 63; 176 h ^= b >> 1; 177 } 178 179 *rh = h; 180 *rl = l; 181 } 182 #endif 183 184 /* Compute xor-multiply of two binary polynomials (a1, a0) x (b1, b0) 185 * result is a binary polynomial in 4 mp_digits r[4]. 186 * The caller MUST ensure that r has the right amount of space allocated. 187 */ 188 void 189 s_bmul_2x2(mp_digit *r, const mp_digit a1, const mp_digit a0, const mp_digit b1, 190 const mp_digit b0) 191 { 192 mp_digit m1, m0; 193 /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */ 194 s_bmul_1x1(r + 3, r + 2, a1, b1); 195 s_bmul_1x1(r + 1, r, a0, b0); 196 s_bmul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1); 197 /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */ 198 r[2] ^= m1 ^ r[1] ^ r[3]; /* h0 ^= m1 ^ l1 ^ h1; */ 199 r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */ 200 } 201 202 /* Compute xor-multiply of two binary polynomials (a2, a1, a0) x (b2, b1, b0) 203 * result is a binary polynomial in 6 mp_digits r[6]. 204 * The caller MUST ensure that r has the right amount of space allocated. 205 */ 206 void 207 s_bmul_3x3(mp_digit *r, const mp_digit a2, const mp_digit a1, const mp_digit a0, 208 const mp_digit b2, const mp_digit b1, const mp_digit b0) 209 { 210 mp_digit zm[4]; 211 212 s_bmul_1x1(r + 5, r + 4, a2, b2); /* fill top 2 words */ 213 s_bmul_2x2(zm, a1, a2 ^ a0, b1, b2 ^ b0); /* fill middle 4 words */ 214 s_bmul_2x2(r, a1, a0, b1, b0); /* fill bottom 4 words */ 215 216 zm[3] ^= r[3]; 217 zm[2] ^= r[2]; 218 zm[1] ^= r[1] ^ r[5]; 219 zm[0] ^= r[0] ^ r[4]; 220 221 r[5] ^= zm[3]; 222 r[4] ^= zm[2]; 223 r[3] ^= zm[1]; 224 r[2] ^= zm[0]; 225 } 226 227 /* Compute xor-multiply of two binary polynomials (a3, a2, a1, a0) x (b3, b2, b1, b0) 228 * result is a binary polynomial in 8 mp_digits r[8]. 229 * The caller MUST ensure that r has the right amount of space allocated. 230 */ 231 void 232 s_bmul_4x4(mp_digit *r, const mp_digit a3, const mp_digit a2, const mp_digit a1, 233 const mp_digit a0, const mp_digit b3, const mp_digit b2, const mp_digit b1, 234 const mp_digit b0) 235 { 236 mp_digit zm[4]; 237 238 s_bmul_2x2(r + 4, a3, a2, b3, b2); /* fill top 4 words */ 239 s_bmul_2x2(zm, a3 ^ a1, a2 ^ a0, b3 ^ b1, b2 ^ b0); /* fill middle 4 words */ 240 s_bmul_2x2(r, a1, a0, b1, b0); /* fill bottom 4 words */ 241 242 zm[3] ^= r[3] ^ r[7]; 243 zm[2] ^= r[2] ^ r[6]; 244 zm[1] ^= r[1] ^ r[5]; 245 zm[0] ^= r[0] ^ r[4]; 246 247 r[5] ^= zm[3]; 248 r[4] ^= zm[2]; 249 r[3] ^= zm[1]; 250 r[2] ^= zm[0]; 251 } 252 253 /* Compute addition of two binary polynomials a and b, 254 * store result in c; c could be a or b, a and b could be equal; 255 * c is the bitwise XOR of a and b. 256 */ 257 mp_err 258 mp_badd(const mp_int *a, const mp_int *b, mp_int *c) 259 { 260 mp_digit *pa, *pb, *pc; 261 mp_size ix; 262 mp_size used_pa, used_pb; 263 mp_err res = MP_OKAY; 264 265 /* Add all digits up to the precision of b. If b had more 266 * precision than a initially, swap a, b first 267 */ 268 if (MP_USED(a) >= MP_USED(b)) { 269 pa = MP_DIGITS(a); 270 pb = MP_DIGITS(b); 271 used_pa = MP_USED(a); 272 used_pb = MP_USED(b); 273 } else { 274 pa = MP_DIGITS(b); 275 pb = MP_DIGITS(a); 276 used_pa = MP_USED(b); 277 used_pb = MP_USED(a); 278 } 279 280 /* Make sure c has enough precision for the output value */ 281 MP_CHECKOK(s_mp_pad(c, used_pa)); 282 283 /* Do word-by-word xor */ 284 pc = MP_DIGITS(c); 285 for (ix = 0; ix < used_pb; ix++) { 286 (*pc++) = (*pa++) ^ (*pb++); 287 } 288 289 /* Finish the rest of digits until we're actually done */ 290 for (; ix < used_pa; ++ix) { 291 *pc++ = *pa++; 292 } 293 294 MP_USED(c) = used_pa; 295 MP_SIGN(c) = ZPOS; 296 s_mp_clamp(c); 297 298 CLEANUP: 299 return res; 300 } 301 302 #define s_mp_div2(a) MP_CHECKOK(mpl_rsh((a), (a), 1)); 303 304 /* Compute binary polynomial multiply d = a * b */ 305 static void 306 s_bmul_d(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d) 307 { 308 mp_digit a_i, a0b0, a1b1, carry = 0; 309 while (a_len--) { 310 a_i = *a++; 311 s_bmul_1x1(&a1b1, &a0b0, a_i, b); 312 *d++ = a0b0 ^ carry; 313 carry = a1b1; 314 } 315 *d = carry; 316 } 317 318 /* Compute binary polynomial xor multiply accumulate d ^= a * b */ 319 static void 320 s_bmul_d_add(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d) 321 { 322 mp_digit a_i, a0b0, a1b1, carry = 0; 323 while (a_len--) { 324 a_i = *a++; 325 s_bmul_1x1(&a1b1, &a0b0, a_i, b); 326 *d++ ^= a0b0 ^ carry; 327 carry = a1b1; 328 } 329 *d ^= carry; 330 } 331 332 /* Compute binary polynomial xor multiply c = a * b. 333 * All parameters may be identical. 334 */ 335 mp_err 336 mp_bmul(const mp_int *a, const mp_int *b, mp_int *c) 337 { 338 mp_digit *pb, b_i; 339 mp_int tmp; 340 mp_size ib, a_used, b_used; 341 mp_err res = MP_OKAY; 342 343 MP_DIGITS(&tmp) = 0; 344 345 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); 346 347 if (a == c) { 348 MP_CHECKOK(mp_init_copy(&tmp, a)); 349 if (a == b) 350 b = &tmp; 351 a = &tmp; 352 } else if (b == c) { 353 MP_CHECKOK(mp_init_copy(&tmp, b)); 354 b = &tmp; 355 } 356 357 if (MP_USED(a) < MP_USED(b)) { 358 const mp_int *xch = b; /* switch a and b if b longer */ 359 b = a; 360 a = xch; 361 } 362 363 MP_USED(c) = 1; 364 MP_DIGIT(c, 0) = 0; 365 MP_CHECKOK(s_mp_pad(c, USED(a) + USED(b))); 366 367 pb = MP_DIGITS(b); 368 s_bmul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c)); 369 370 /* Outer loop: Digits of b */ 371 a_used = MP_USED(a); 372 b_used = MP_USED(b); 373 MP_USED(c) = a_used + b_used; 374 for (ib = 1; ib < b_used; ib++) { 375 b_i = *pb++; 376 377 /* Inner product: Digits of a */ 378 if (b_i) 379 s_bmul_d_add(MP_DIGITS(a), a_used, b_i, MP_DIGITS(c) + ib); 380 else 381 MP_DIGIT(c, ib + a_used) = b_i; 382 } 383 384 s_mp_clamp(c); 385 386 SIGN(c) = ZPOS; 387 388 CLEANUP: 389 mp_clear(&tmp); 390 return res; 391 } 392 393 /* Compute modular reduction of a and store result in r. 394 * r could be a. 395 * For modular arithmetic, the irreducible polynomial f(t) is represented 396 * as an array of int[], where f(t) is of the form: 397 * f(t) = t^p[0] + t^p[1] + ... + t^p[k] 398 * where m = p[0] > p[1] > ... > p[k] = 0. 399 */ 400 mp_err 401 mp_bmod(const mp_int *a, const unsigned int p[], mp_int *r) 402 { 403 int j, k; 404 int n, dN, d0, d1; 405 mp_digit zz, *z, tmp; 406 mp_size used; 407 mp_err res = MP_OKAY; 408 409 /* The algorithm does the reduction in place in r, 410 * if a != r, copy a into r first so reduction can be done in r 411 */ 412 if (a != r) { 413 MP_CHECKOK(mp_copy(a, r)); 414 } 415 z = MP_DIGITS(r); 416 417 /* start reduction */ 418 /*dN = p[0] / MP_DIGIT_BITS; */ 419 dN = p[0] >> MP_DIGIT_BITS_LOG_2; 420 used = MP_USED(r); 421 422 for (j = used - 1; j > dN;) { 423 424 zz = z[j]; 425 if (zz == 0) { 426 j--; 427 continue; 428 } 429 z[j] = 0; 430 431 for (k = 1; p[k] > 0; k++) { 432 /* reducing component t^p[k] */ 433 n = p[0] - p[k]; 434 /*d0 = n % MP_DIGIT_BITS; */ 435 d0 = n & MP_DIGIT_BITS_MASK; 436 d1 = MP_DIGIT_BITS - d0; 437 /*n /= MP_DIGIT_BITS; */ 438 n >>= MP_DIGIT_BITS_LOG_2; 439 z[j - n] ^= (zz >> d0); 440 if (d0) 441 z[j - n - 1] ^= (zz << d1); 442 } 443 444 /* reducing component t^0 */ 445 n = dN; 446 /*d0 = p[0] % MP_DIGIT_BITS;*/ 447 d0 = p[0] & MP_DIGIT_BITS_MASK; 448 d1 = MP_DIGIT_BITS - d0; 449 z[j - n] ^= (zz >> d0); 450 if (d0) 451 z[j - n - 1] ^= (zz << d1); 452 } 453 454 /* final round of reduction */ 455 while (j == dN) { 456 457 /* d0 = p[0] % MP_DIGIT_BITS; */ 458 d0 = p[0] & MP_DIGIT_BITS_MASK; 459 zz = z[dN] >> d0; 460 if (zz == 0) 461 break; 462 d1 = MP_DIGIT_BITS - d0; 463 464 /* clear up the top d1 bits */ 465 if (d0) { 466 z[dN] = (z[dN] << d1) >> d1; 467 } else { 468 z[dN] = 0; 469 } 470 *z ^= zz; /* reduction t^0 component */ 471 472 for (k = 1; p[k] > 0; k++) { 473 /* reducing component t^p[k]*/ 474 /* n = p[k] / MP_DIGIT_BITS; */ 475 n = p[k] >> MP_DIGIT_BITS_LOG_2; 476 /* d0 = p[k] % MP_DIGIT_BITS; */ 477 d0 = p[k] & MP_DIGIT_BITS_MASK; 478 d1 = MP_DIGIT_BITS - d0; 479 z[n] ^= (zz << d0); 480 tmp = zz >> d1; 481 if (d0 && tmp) 482 z[n + 1] ^= tmp; 483 } 484 } 485 486 s_mp_clamp(r); 487 CLEANUP: 488 return res; 489 } 490 491 /* Compute the product of two polynomials a and b, reduce modulo p, 492 * Store the result in r. r could be a or b; a could be b. 493 */ 494 mp_err 495 mp_bmulmod(const mp_int *a, const mp_int *b, const unsigned int p[], mp_int *r) 496 { 497 mp_err res; 498 499 if (a == b) 500 return mp_bsqrmod(a, p, r); 501 if ((res = mp_bmul(a, b, r)) != MP_OKAY) 502 return res; 503 return mp_bmod(r, p, r); 504 } 505 506 /* Compute binary polynomial squaring c = a*a mod p . 507 * Parameter r and a can be identical. 508 */ 509 510 mp_err 511 mp_bsqrmod(const mp_int *a, const unsigned int p[], mp_int *r) 512 { 513 mp_digit *pa, *pr, a_i; 514 mp_int tmp; 515 mp_size ia, a_used; 516 mp_err res; 517 518 ARGCHK(a != NULL && r != NULL, MP_BADARG); 519 MP_DIGITS(&tmp) = 0; 520 521 if (a == r) { 522 MP_CHECKOK(mp_init_copy(&tmp, a)); 523 a = &tmp; 524 } 525 526 MP_USED(r) = 1; 527 MP_DIGIT(r, 0) = 0; 528 MP_CHECKOK(s_mp_pad(r, 2 * USED(a))); 529 530 pa = MP_DIGITS(a); 531 pr = MP_DIGITS(r); 532 a_used = MP_USED(a); 533 MP_USED(r) = 2 * a_used; 534 535 for (ia = 0; ia < a_used; ia++) { 536 a_i = *pa++; 537 *pr++ = gf2m_SQR0(a_i); 538 *pr++ = gf2m_SQR1(a_i); 539 } 540 541 MP_CHECKOK(mp_bmod(r, p, r)); 542 s_mp_clamp(r); 543 SIGN(r) = ZPOS; 544 545 CLEANUP: 546 mp_clear(&tmp); 547 return res; 548 } 549 550 /* Compute binary polynomial y/x mod p, y divided by x, reduce modulo p. 551 * Store the result in r. r could be x or y, and x could equal y. 552 * Uses algorithm Modular_Division_GF(2^m) from 553 * Chang-Shantz, S. "From Euclid's GCD to Montgomery Multiplication to 554 * the Great Divide". 555 */ 556 int 557 mp_bdivmod(const mp_int *y, const mp_int *x, const mp_int *pp, 558 const unsigned int p[], mp_int *r) 559 { 560 mp_int aa, bb, uu; 561 mp_int *a, *b, *u, *v; 562 mp_err res = MP_OKAY; 563 564 MP_DIGITS(&aa) = 0; 565 MP_DIGITS(&bb) = 0; 566 MP_DIGITS(&uu) = 0; 567 568 MP_CHECKOK(mp_init_copy(&aa, x)); 569 MP_CHECKOK(mp_init_copy(&uu, y)); 570 MP_CHECKOK(mp_init_copy(&bb, pp)); 571 MP_CHECKOK(s_mp_pad(r, USED(pp))); 572 MP_USED(r) = 1; 573 MP_DIGIT(r, 0) = 0; 574 575 a = &aa; 576 b = &bb; 577 u = &uu; 578 v = r; 579 /* reduce x and y mod p */ 580 MP_CHECKOK(mp_bmod(a, p, a)); 581 MP_CHECKOK(mp_bmod(u, p, u)); 582 583 while (!mp_isodd(a)) { 584 s_mp_div2(a); 585 if (mp_isodd(u)) { 586 MP_CHECKOK(mp_badd(u, pp, u)); 587 } 588 s_mp_div2(u); 589 } 590 591 do { 592 if (mp_cmp_mag(b, a) > 0) { 593 MP_CHECKOK(mp_badd(b, a, b)); 594 MP_CHECKOK(mp_badd(v, u, v)); 595 do { 596 s_mp_div2(b); 597 if (mp_isodd(v)) { 598 MP_CHECKOK(mp_badd(v, pp, v)); 599 } 600 s_mp_div2(v); 601 } while (!mp_isodd(b)); 602 } else if ((MP_DIGIT(a, 0) == 1) && (MP_USED(a) == 1)) 603 break; 604 else { 605 MP_CHECKOK(mp_badd(a, b, a)); 606 MP_CHECKOK(mp_badd(u, v, u)); 607 do { 608 s_mp_div2(a); 609 if (mp_isodd(u)) { 610 MP_CHECKOK(mp_badd(u, pp, u)); 611 } 612 s_mp_div2(u); 613 } while (!mp_isodd(a)); 614 } 615 } while (1); 616 617 MP_CHECKOK(mp_copy(u, r)); 618 619 CLEANUP: 620 mp_clear(&aa); 621 mp_clear(&bb); 622 mp_clear(&uu); 623 return res; 624 } 625 626 /* Convert the bit-string representation of a polynomial a into an array 627 * of integers corresponding to the bits with non-zero coefficient. 628 * Up to max elements of the array will be filled. Return value is total 629 * number of coefficients that would be extracted if array was large enough. 630 */ 631 int 632 mp_bpoly2arr(const mp_int *a, unsigned int p[], int max) 633 { 634 int i, j, k; 635 mp_digit top_bit, mask; 636 637 top_bit = 1; 638 top_bit <<= MP_DIGIT_BIT - 1; 639 640 for (k = 0; k < max; k++) 641 p[k] = 0; 642 k = 0; 643 644 for (i = MP_USED(a) - 1; i >= 0; i--) { 645 mask = top_bit; 646 for (j = MP_DIGIT_BIT - 1; j >= 0; j--) { 647 if (MP_DIGITS(a)[i] & mask) { 648 if (k < max) 649 p[k] = MP_DIGIT_BIT * i + j; 650 k++; 651 } 652 mask >>= 1; 653 } 654 } 655 656 return k; 657 } 658 659 /* Convert the coefficient array representation of a polynomial to a 660 * bit-string. The array must be terminated by 0. 661 */ 662 mp_err 663 mp_barr2poly(const unsigned int p[], mp_int *a) 664 { 665 666 mp_err res = MP_OKAY; 667 int i; 668 669 mp_zero(a); 670 for (i = 0; p[i] > 0; i++) { 671 MP_CHECKOK(mpl_set_bit(a, p[i], 1)); 672 } 673 MP_CHECKOK(mpl_set_bit(a, 0, 1)); 674 675 CLEANUP: 676 return res; 677 }