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