1 /* 2 * Minimal code for RSA support from LibTomMath 0.3.9 3 * http://math.libtomcrypt.com/ 4 * http://math.libtomcrypt.com/files/ltm-0.39.tar.bz2 5 * This library was released in public domain by Tom St Denis. 6 * 7 * The combination in this file is not using many of the optimized algorithms 8 * (e.g., Montgomery reduction) and is considerable slower than the LibTomMath 9 * with its default of SC_RSA_1 settins. The main purpose of having this 10 * version here is to make it easier to build bignum.c wrapper without having 11 * to install and build an external library. However, it is likely worth the 12 * effort to use the full library with SC_RSA_1 instead of this minimized copy. 13 * Including the optimized algorithms may increase the size requirements by 14 * 15 kB or so (measured with x86 build). 15 * 16 * If CONFIG_INTERNAL_LIBTOMMATH is defined, bignum.c includes this 17 * libtommath.c file instead of using the external LibTomMath library. 18 */ 19 20 #ifndef CHAR_BIT 21 #define CHAR_BIT 8 22 #endif 23 24 #define BN_MP_INVMOD_C 25 #define BN_S_MP_EXPTMOD_C /* Note: #undef in tommath_superclass.h; this would 26 * require BN_MP_EXPTMOD_FAST_C instead */ 27 #define BN_S_MP_MUL_DIGS_C 28 #define BN_MP_INVMOD_SLOW_C 29 #define BN_S_MP_SQR_C 30 #define BN_S_MP_MUL_HIGH_DIGS_C /* Note: #undef in tommath_superclass.h; this 31 * would require other than mp_reduce */ 32 33 34 /* from tommath.h */ 35 36 #ifndef MIN 37 #define MIN(x,y) ((x)<(y)?(x):(y)) 38 #endif 39 40 #ifndef MAX 41 #define MAX(x,y) ((x)>(y)?(x):(y)) 42 #endif 43 44 #define OPT_CAST(x) 45 46 typedef unsigned long mp_digit; 47 typedef u64 mp_word; 48 49 #define DIGIT_BIT 28 50 #define MP_28BIT 51 52 53 #define XMALLOC os_malloc 54 #define XFREE os_free 55 #define XREALLOC os_realloc 56 57 58 #define MP_MASK ((((mp_digit)1)<<((mp_digit)DIGIT_BIT))-((mp_digit)1)) 59 60 #define MP_LT -1 /* less than */ 61 #define MP_EQ 0 /* equal to */ 62 #define MP_GT 1 /* greater than */ 63 64 #define MP_ZPOS 0 /* positive integer */ 65 #define MP_NEG 1 /* negative */ 66 67 #define MP_OKAY 0 /* ok result */ 68 #define MP_MEM -2 /* out of mem */ 69 #define MP_VAL -3 /* invalid input */ 70 71 #define MP_YES 1 /* yes response */ 72 #define MP_NO 0 /* no response */ 73 74 typedef int mp_err; 75 76 /* define this to use lower memory usage routines (exptmods mostly) */ 77 #define MP_LOW_MEM 78 79 /* default precision */ 80 #ifndef MP_PREC 81 #ifndef MP_LOW_MEM 82 #define MP_PREC 32 /* default digits of precision */ 83 #else 84 #define MP_PREC 8 /* default digits of precision */ 85 #endif 86 #endif 87 88 /* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */ 89 #define MP_WARRAY (1 << (sizeof(mp_word) * CHAR_BIT - 2 * DIGIT_BIT + 1)) 90 91 /* the infamous mp_int structure */ 92 typedef struct { 93 int used, alloc, sign; 94 mp_digit *dp; 95 } mp_int; 96 97 98 /* ---> Basic Manipulations <--- */ 99 #define mp_iszero(a) (((a)->used == 0) ? MP_YES : MP_NO) 100 #define mp_iseven(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 0)) ? MP_YES : MP_NO) 101 #define mp_isodd(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? MP_YES : MP_NO) 102 103 104 /* prototypes for copied functions */ 105 #define s_mp_mul(a, b, c) s_mp_mul_digs(a, b, c, (a)->used + (b)->used + 1) 106 static int s_mp_exptmod(mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode); 107 static int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs); 108 static int s_mp_sqr(mp_int * a, mp_int * b); 109 static int s_mp_mul_high_digs(mp_int * a, mp_int * b, mp_int * c, int digs); 110 111 static int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs); 112 113 static int mp_init_multi(mp_int *mp, ...); 114 static void mp_clear_multi(mp_int *mp, ...); 115 static int mp_lshd(mp_int * a, int b); 116 static void mp_set(mp_int * a, mp_digit b); 117 static void mp_clamp(mp_int * a); 118 static void mp_exch(mp_int * a, mp_int * b); 119 static void mp_rshd(mp_int * a, int b); 120 static void mp_zero(mp_int * a); 121 static int mp_mod_2d(mp_int * a, int b, mp_int * c); 122 static int mp_div_2d(mp_int * a, int b, mp_int * c, mp_int * d); 123 static int mp_init_copy(mp_int * a, mp_int * b); 124 static int mp_mul_2d(mp_int * a, int b, mp_int * c); 125 static int mp_div_2(mp_int * a, mp_int * b); 126 static int mp_copy(mp_int * a, mp_int * b); 127 static int mp_count_bits(mp_int * a); 128 static int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d); 129 static int mp_mod(mp_int * a, mp_int * b, mp_int * c); 130 static int mp_grow(mp_int * a, int size); 131 static int mp_cmp_mag(mp_int * a, mp_int * b); 132 static int mp_invmod(mp_int * a, mp_int * b, mp_int * c); 133 static int mp_abs(mp_int * a, mp_int * b); 134 static int mp_invmod_slow(mp_int * a, mp_int * b, mp_int * c); 135 static int mp_sqr(mp_int * a, mp_int * b); 136 static int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d); 137 static int mp_reduce_2k_setup_l(mp_int *a, mp_int *d); 138 static int mp_2expt(mp_int * a, int b); 139 static int mp_reduce_setup(mp_int * a, mp_int * b); 140 static int mp_reduce(mp_int * x, mp_int * m, mp_int * mu); 141 static int mp_init_size(mp_int * a, int size); 142 143 144 145 /* functions from bn_<func name>.c */ 146 147 148 /* reverse an array, used for radix code */ 149 static void bn_reverse (unsigned char *s, int len) 150 { 151 int ix, iy; 152 unsigned char t; 153 154 ix = 0; 155 iy = len - 1; 156 while (ix < iy) { 157 t = s[ix]; 158 s[ix] = s[iy]; 159 s[iy] = t; 160 ++ix; 161 --iy; 162 } 163 } 164 165 166 /* low level addition, based on HAC pp.594, Algorithm 14.7 */ 167 static int s_mp_add (mp_int * a, mp_int * b, mp_int * c) 168 { 169 mp_int *x; 170 int olduse, res, min, max; 171 172 /* find sizes, we let |a| <= |b| which means we have to sort 173 * them. "x" will point to the input with the most digits 174 */ 175 if (a->used > b->used) { 176 min = b->used; 177 max = a->used; 178 x = a; 179 } else { 180 min = a->used; 181 max = b->used; 182 x = b; 183 } 184 185 /* init result */ 186 if (c->alloc < max + 1) { 187 if ((res = mp_grow (c, max + 1)) != MP_OKAY) { 188 return res; 189 } 190 } 191 192 /* get old used digit count and set new one */ 193 olduse = c->used; 194 c->used = max + 1; 195 196 { 197 register mp_digit u, *tmpa, *tmpb, *tmpc; 198 register int i; 199 200 /* alias for digit pointers */ 201 202 /* first input */ 203 tmpa = a->dp; 204 205 /* second input */ 206 tmpb = b->dp; 207 208 /* destination */ 209 tmpc = c->dp; 210 211 /* zero the carry */ 212 u = 0; 213 for (i = 0; i < min; i++) { 214 /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */ 215 *tmpc = *tmpa++ + *tmpb++ + u; 216 217 /* U = carry bit of T[i] */ 218 u = *tmpc >> ((mp_digit)DIGIT_BIT); 219 220 /* take away carry bit from T[i] */ 221 *tmpc++ &= MP_MASK; 222 } 223 224 /* now copy higher words if any, that is in A+B 225 * if A or B has more digits add those in 226 */ 227 if (min != max) { 228 for (; i < max; i++) { 229 /* T[i] = X[i] + U */ 230 *tmpc = x->dp[i] + u; 231 232 /* U = carry bit of T[i] */ 233 u = *tmpc >> ((mp_digit)DIGIT_BIT); 234 235 /* take away carry bit from T[i] */ 236 *tmpc++ &= MP_MASK; 237 } 238 } 239 240 /* add carry */ 241 *tmpc++ = u; 242 243 /* clear digits above oldused */ 244 for (i = c->used; i < olduse; i++) { 245 *tmpc++ = 0; 246 } 247 } 248 249 mp_clamp (c); 250 return MP_OKAY; 251 } 252 253 254 /* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */ 255 static int s_mp_sub (mp_int * a, mp_int * b, mp_int * c) 256 { 257 int olduse, res, min, max; 258 259 /* find sizes */ 260 min = b->used; 261 max = a->used; 262 263 /* init result */ 264 if (c->alloc < max) { 265 if ((res = mp_grow (c, max)) != MP_OKAY) { 266 return res; 267 } 268 } 269 olduse = c->used; 270 c->used = max; 271 272 { 273 register mp_digit u, *tmpa, *tmpb, *tmpc; 274 register int i; 275 276 /* alias for digit pointers */ 277 tmpa = a->dp; 278 tmpb = b->dp; 279 tmpc = c->dp; 280 281 /* set carry to zero */ 282 u = 0; 283 for (i = 0; i < min; i++) { 284 /* T[i] = A[i] - B[i] - U */ 285 *tmpc = *tmpa++ - *tmpb++ - u; 286 287 /* U = carry bit of T[i] 288 * Note this saves performing an AND operation since 289 * if a carry does occur it will propagate all the way to the 290 * MSB. As a result a single shift is enough to get the carry 291 */ 292 u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1)); 293 294 /* Clear carry from T[i] */ 295 *tmpc++ &= MP_MASK; 296 } 297 298 /* now copy higher words if any, e.g. if A has more digits than B */ 299 for (; i < max; i++) { 300 /* T[i] = A[i] - U */ 301 *tmpc = *tmpa++ - u; 302 303 /* U = carry bit of T[i] */ 304 u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1)); 305 306 /* Clear carry from T[i] */ 307 *tmpc++ &= MP_MASK; 308 } 309 310 /* clear digits above used (since we may not have grown result above) */ 311 for (i = c->used; i < olduse; i++) { 312 *tmpc++ = 0; 313 } 314 } 315 316 mp_clamp (c); 317 return MP_OKAY; 318 } 319 320 321 /* init a new mp_int */ 322 static int mp_init (mp_int * a) 323 { 324 int i; 325 326 /* allocate memory required and clear it */ 327 a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * MP_PREC); 328 if (a->dp == NULL) { 329 return MP_MEM; 330 } 331 332 /* set the digits to zero */ 333 for (i = 0; i < MP_PREC; i++) { 334 a->dp[i] = 0; 335 } 336 337 /* set the used to zero, allocated digits to the default precision 338 * and sign to positive */ 339 a->used = 0; 340 a->alloc = MP_PREC; 341 a->sign = MP_ZPOS; 342 343 return MP_OKAY; 344 } 345 346 347 /* clear one (frees) */ 348 static void mp_clear (mp_int * a) 349 { 350 int i; 351 352 /* only do anything if a hasn't been freed previously */ 353 if (a->dp != NULL) { 354 /* first zero the digits */ 355 for (i = 0; i < a->used; i++) { 356 a->dp[i] = 0; 357 } 358 359 /* free ram */ 360 XFREE(a->dp); 361 362 /* reset members to make debugging easier */ 363 a->dp = NULL; 364 a->alloc = a->used = 0; 365 a->sign = MP_ZPOS; 366 } 367 } 368 369 370 /* high level addition (handles signs) */ 371 static int mp_add (mp_int * a, mp_int * b, mp_int * c) 372 { 373 int sa, sb, res; 374 375 /* get sign of both inputs */ 376 sa = a->sign; 377 sb = b->sign; 378 379 /* handle two cases, not four */ 380 if (sa == sb) { 381 /* both positive or both negative */ 382 /* add their magnitudes, copy the sign */ 383 c->sign = sa; 384 res = s_mp_add (a, b, c); 385 } else { 386 /* one positive, the other negative */ 387 /* subtract the one with the greater magnitude from */ 388 /* the one of the lesser magnitude. The result gets */ 389 /* the sign of the one with the greater magnitude. */ 390 if (mp_cmp_mag (a, b) == MP_LT) { 391 c->sign = sb; 392 res = s_mp_sub (b, a, c); 393 } else { 394 c->sign = sa; 395 res = s_mp_sub (a, b, c); 396 } 397 } 398 return res; 399 } 400 401 402 /* high level subtraction (handles signs) */ 403 static int mp_sub (mp_int * a, mp_int * b, mp_int * c) 404 { 405 int sa, sb, res; 406 407 sa = a->sign; 408 sb = b->sign; 409 410 if (sa != sb) { 411 /* subtract a negative from a positive, OR */ 412 /* subtract a positive from a negative. */ 413 /* In either case, ADD their magnitudes, */ 414 /* and use the sign of the first number. */ 415 c->sign = sa; 416 res = s_mp_add (a, b, c); 417 } else { 418 /* subtract a positive from a positive, OR */ 419 /* subtract a negative from a negative. */ 420 /* First, take the difference between their */ 421 /* magnitudes, then... */ 422 if (mp_cmp_mag (a, b) != MP_LT) { 423 /* Copy the sign from the first */ 424 c->sign = sa; 425 /* The first has a larger or equal magnitude */ 426 res = s_mp_sub (a, b, c); 427 } else { 428 /* The result has the *opposite* sign from */ 429 /* the first number. */ 430 c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS; 431 /* The second has a larger magnitude */ 432 res = s_mp_sub (b, a, c); 433 } 434 } 435 return res; 436 } 437 438 439 /* high level multiplication (handles sign) */ 440 static int mp_mul (mp_int * a, mp_int * b, mp_int * c) 441 { 442 int res, neg; 443 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 444 445 /* use Toom-Cook? */ 446 #ifdef BN_MP_TOOM_MUL_C 447 if (MIN (a->used, b->used) >= TOOM_MUL_CUTOFF) { 448 res = mp_toom_mul(a, b, c); 449 } else 450 #endif 451 #ifdef BN_MP_KARATSUBA_MUL_C 452 /* use Karatsuba? */ 453 if (MIN (a->used, b->used) >= KARATSUBA_MUL_CUTOFF) { 454 res = mp_karatsuba_mul (a, b, c); 455 } else 456 #endif 457 { 458 /* can we use the fast multiplier? 459 * 460 * The fast multiplier can be used if the output will 461 * have less than MP_WARRAY digits and the number of 462 * digits won't affect carry propagation 463 */ 464 #ifdef BN_FAST_S_MP_MUL_DIGS_C 465 int digs = a->used + b->used + 1; 466 467 if ((digs < MP_WARRAY) && 468 MIN(a->used, b->used) <= 469 (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { 470 res = fast_s_mp_mul_digs (a, b, c, digs); 471 } else 472 #endif 473 #ifdef BN_S_MP_MUL_DIGS_C 474 res = s_mp_mul (a, b, c); /* uses s_mp_mul_digs */ 475 #else 476 #error mp_mul could fail 477 res = MP_VAL; 478 #endif 479 480 } 481 c->sign = (c->used > 0) ? neg : MP_ZPOS; 482 return res; 483 } 484 485 486 /* d = a * b (mod c) */ 487 static int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d) 488 { 489 int res; 490 mp_int t; 491 492 if ((res = mp_init (&t)) != MP_OKAY) { 493 return res; 494 } 495 496 if ((res = mp_mul (a, b, &t)) != MP_OKAY) { 497 mp_clear (&t); 498 return res; 499 } 500 res = mp_mod (&t, c, d); 501 mp_clear (&t); 502 return res; 503 } 504 505 506 /* c = a mod b, 0 <= c < b */ 507 static int mp_mod (mp_int * a, mp_int * b, mp_int * c) 508 { 509 mp_int t; 510 int res; 511 512 if ((res = mp_init (&t)) != MP_OKAY) { 513 return res; 514 } 515 516 if ((res = mp_div (a, b, NULL, &t)) != MP_OKAY) { 517 mp_clear (&t); 518 return res; 519 } 520 521 if (t.sign != b->sign) { 522 res = mp_add (b, &t, c); 523 } else { 524 res = MP_OKAY; 525 mp_exch (&t, c); 526 } 527 528 mp_clear (&t); 529 return res; 530 } 531 532 533 /* this is a shell function that calls either the normal or Montgomery 534 * exptmod functions. Originally the call to the montgomery code was 535 * embedded in the normal function but that wasted alot of stack space 536 * for nothing (since 99% of the time the Montgomery code would be called) 537 */ 538 static int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y) 539 { 540 int dr; 541 542 /* modulus P must be positive */ 543 if (P->sign == MP_NEG) { 544 return MP_VAL; 545 } 546 547 /* if exponent X is negative we have to recurse */ 548 if (X->sign == MP_NEG) { 549 #ifdef BN_MP_INVMOD_C 550 mp_int tmpG, tmpX; 551 int err; 552 553 /* first compute 1/G mod P */ 554 if ((err = mp_init(&tmpG)) != MP_OKAY) { 555 return err; 556 } 557 if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) { 558 mp_clear(&tmpG); 559 return err; 560 } 561 562 /* now get |X| */ 563 if ((err = mp_init(&tmpX)) != MP_OKAY) { 564 mp_clear(&tmpG); 565 return err; 566 } 567 if ((err = mp_abs(X, &tmpX)) != MP_OKAY) { 568 mp_clear_multi(&tmpG, &tmpX, NULL); 569 return err; 570 } 571 572 /* and now compute (1/G)**|X| instead of G**X [X < 0] */ 573 err = mp_exptmod(&tmpG, &tmpX, P, Y); 574 mp_clear_multi(&tmpG, &tmpX, NULL); 575 return err; 576 #else 577 #error mp_exptmod would always fail 578 /* no invmod */ 579 return MP_VAL; 580 #endif 581 } 582 583 /* modified diminished radix reduction */ 584 #if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C) && defined(BN_S_MP_EXPTMOD_C) 585 if (mp_reduce_is_2k_l(P) == MP_YES) { 586 return s_mp_exptmod(G, X, P, Y, 1); 587 } 588 #endif 589 590 #ifdef BN_MP_DR_IS_MODULUS_C 591 /* is it a DR modulus? */ 592 dr = mp_dr_is_modulus(P); 593 #else 594 /* default to no */ 595 dr = 0; 596 #endif 597 598 #ifdef BN_MP_REDUCE_IS_2K_C 599 /* if not, is it a unrestricted DR modulus? */ 600 if (dr == 0) { 601 dr = mp_reduce_is_2k(P) << 1; 602 } 603 #endif 604 605 /* if the modulus is odd or dr != 0 use the montgomery method */ 606 #ifdef BN_MP_EXPTMOD_FAST_C 607 if (mp_isodd (P) == 1 || dr != 0) { 608 return mp_exptmod_fast (G, X, P, Y, dr); 609 } else { 610 #endif 611 #ifdef BN_S_MP_EXPTMOD_C 612 /* otherwise use the generic Barrett reduction technique */ 613 return s_mp_exptmod (G, X, P, Y, 0); 614 #else 615 #error mp_exptmod could fail 616 /* no exptmod for evens */ 617 return MP_VAL; 618 #endif 619 #ifdef BN_MP_EXPTMOD_FAST_C 620 } 621 #endif 622 } 623 624 625 /* compare two ints (signed)*/ 626 static int mp_cmp (mp_int * a, mp_int * b) 627 { 628 /* compare based on sign */ 629 if (a->sign != b->sign) { 630 if (a->sign == MP_NEG) { 631 return MP_LT; 632 } else { 633 return MP_GT; 634 } 635 } 636 637 /* compare digits */ 638 if (a->sign == MP_NEG) { 639 /* if negative compare opposite direction */ 640 return mp_cmp_mag(b, a); 641 } else { 642 return mp_cmp_mag(a, b); 643 } 644 } 645 646 647 /* compare a digit */ 648 static int mp_cmp_d(mp_int * a, mp_digit b) 649 { 650 /* compare based on sign */ 651 if (a->sign == MP_NEG) { 652 return MP_LT; 653 } 654 655 /* compare based on magnitude */ 656 if (a->used > 1) { 657 return MP_GT; 658 } 659 660 /* compare the only digit of a to b */ 661 if (a->dp[0] > b) { 662 return MP_GT; 663 } else if (a->dp[0] < b) { 664 return MP_LT; 665 } else { 666 return MP_EQ; 667 } 668 } 669 670 671 /* hac 14.61, pp608 */ 672 static int mp_invmod (mp_int * a, mp_int * b, mp_int * c) 673 { 674 /* b cannot be negative */ 675 if (b->sign == MP_NEG || mp_iszero(b) == 1) { 676 return MP_VAL; 677 } 678 679 #ifdef BN_FAST_MP_INVMOD_C 680 /* if the modulus is odd we can use a faster routine instead */ 681 if (mp_isodd (b) == 1) { 682 return fast_mp_invmod (a, b, c); 683 } 684 #endif 685 686 #ifdef BN_MP_INVMOD_SLOW_C 687 return mp_invmod_slow(a, b, c); 688 #endif 689 690 #ifndef BN_FAST_MP_INVMOD_C 691 #ifndef BN_MP_INVMOD_SLOW_C 692 #error mp_invmod would always fail 693 #endif 694 #endif 695 return MP_VAL; 696 } 697 698 699 /* get the size for an unsigned equivalent */ 700 static int mp_unsigned_bin_size (mp_int * a) 701 { 702 int size = mp_count_bits (a); 703 return (size / 8 + ((size & 7) != 0 ? 1 : 0)); 704 } 705 706 707 /* hac 14.61, pp608 */ 708 static int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c) 709 { 710 mp_int x, y, u, v, A, B, C, D; 711 int res; 712 713 /* b cannot be negative */ 714 if (b->sign == MP_NEG || mp_iszero(b) == 1) { 715 return MP_VAL; 716 } 717 718 /* init temps */ 719 if ((res = mp_init_multi(&x, &y, &u, &v, 720 &A, &B, &C, &D, NULL)) != MP_OKAY) { 721 return res; 722 } 723 724 /* x = a, y = b */ 725 if ((res = mp_mod(a, b, &x)) != MP_OKAY) { 726 goto LBL_ERR; 727 } 728 if ((res = mp_copy (b, &y)) != MP_OKAY) { 729 goto LBL_ERR; 730 } 731 732 /* 2. [modified] if x,y are both even then return an error! */ 733 if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) { 734 res = MP_VAL; 735 goto LBL_ERR; 736 } 737 738 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */ 739 if ((res = mp_copy (&x, &u)) != MP_OKAY) { 740 goto LBL_ERR; 741 } 742 if ((res = mp_copy (&y, &v)) != MP_OKAY) { 743 goto LBL_ERR; 744 } 745 mp_set (&A, 1); 746 mp_set (&D, 1); 747 748 top: 749 /* 4. while u is even do */ 750 while (mp_iseven (&u) == 1) { 751 /* 4.1 u = u/2 */ 752 if ((res = mp_div_2 (&u, &u)) != MP_OKAY) { 753 goto LBL_ERR; 754 } 755 /* 4.2 if A or B is odd then */ 756 if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) { 757 /* A = (A+y)/2, B = (B-x)/2 */ 758 if ((res = mp_add (&A, &y, &A)) != MP_OKAY) { 759 goto LBL_ERR; 760 } 761 if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) { 762 goto LBL_ERR; 763 } 764 } 765 /* A = A/2, B = B/2 */ 766 if ((res = mp_div_2 (&A, &A)) != MP_OKAY) { 767 goto LBL_ERR; 768 } 769 if ((res = mp_div_2 (&B, &B)) != MP_OKAY) { 770 goto LBL_ERR; 771 } 772 } 773 774 /* 5. while v is even do */ 775 while (mp_iseven (&v) == 1) { 776 /* 5.1 v = v/2 */ 777 if ((res = mp_div_2 (&v, &v)) != MP_OKAY) { 778 goto LBL_ERR; 779 } 780 /* 5.2 if C or D is odd then */ 781 if (mp_isodd (&C) == 1 || mp_isodd (&D) == 1) { 782 /* C = (C+y)/2, D = (D-x)/2 */ 783 if ((res = mp_add (&C, &y, &C)) != MP_OKAY) { 784 goto LBL_ERR; 785 } 786 if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) { 787 goto LBL_ERR; 788 } 789 } 790 /* C = C/2, D = D/2 */ 791 if ((res = mp_div_2 (&C, &C)) != MP_OKAY) { 792 goto LBL_ERR; 793 } 794 if ((res = mp_div_2 (&D, &D)) != MP_OKAY) { 795 goto LBL_ERR; 796 } 797 } 798 799 /* 6. if u >= v then */ 800 if (mp_cmp (&u, &v) != MP_LT) { 801 /* u = u - v, A = A - C, B = B - D */ 802 if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) { 803 goto LBL_ERR; 804 } 805 806 if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) { 807 goto LBL_ERR; 808 } 809 810 if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) { 811 goto LBL_ERR; 812 } 813 } else { 814 /* v - v - u, C = C - A, D = D - B */ 815 if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) { 816 goto LBL_ERR; 817 } 818 819 if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) { 820 goto LBL_ERR; 821 } 822 823 if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) { 824 goto LBL_ERR; 825 } 826 } 827 828 /* if not zero goto step 4 */ 829 if (mp_iszero (&u) == 0) 830 goto top; 831 832 /* now a = C, b = D, gcd == g*v */ 833 834 /* if v != 1 then there is no inverse */ 835 if (mp_cmp_d (&v, 1) != MP_EQ) { 836 res = MP_VAL; 837 goto LBL_ERR; 838 } 839 840 /* if its too low */ 841 while (mp_cmp_d(&C, 0) == MP_LT) { 842 if ((res = mp_add(&C, b, &C)) != MP_OKAY) { 843 goto LBL_ERR; 844 } 845 } 846 847 /* too big */ 848 while (mp_cmp_mag(&C, b) != MP_LT) { 849 if ((res = mp_sub(&C, b, &C)) != MP_OKAY) { 850 goto LBL_ERR; 851 } 852 } 853 854 /* C is now the inverse */ 855 mp_exch (&C, c); 856 res = MP_OKAY; 857 LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL); 858 return res; 859 } 860 861 862 /* compare maginitude of two ints (unsigned) */ 863 static int mp_cmp_mag (mp_int * a, mp_int * b) 864 { 865 int n; 866 mp_digit *tmpa, *tmpb; 867 868 /* compare based on # of non-zero digits */ 869 if (a->used > b->used) { 870 return MP_GT; 871 } 872 873 if (a->used < b->used) { 874 return MP_LT; 875 } 876 877 /* alias for a */ 878 tmpa = a->dp + (a->used - 1); 879 880 /* alias for b */ 881 tmpb = b->dp + (a->used - 1); 882 883 /* compare based on digits */ 884 for (n = 0; n < a->used; ++n, --tmpa, --tmpb) { 885 if (*tmpa > *tmpb) { 886 return MP_GT; 887 } 888 889 if (*tmpa < *tmpb) { 890 return MP_LT; 891 } 892 } 893 return MP_EQ; 894 } 895 896 897 /* reads a unsigned char array, assumes the msb is stored first [big endian] */ 898 static int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c) 899 { 900 int res; 901 902 /* make sure there are at least two digits */ 903 if (a->alloc < 2) { 904 if ((res = mp_grow(a, 2)) != MP_OKAY) { 905 return res; 906 } 907 } 908 909 /* zero the int */ 910 mp_zero (a); 911 912 /* read the bytes in */ 913 while (c-- > 0) { 914 if ((res = mp_mul_2d (a, 8, a)) != MP_OKAY) { 915 return res; 916 } 917 918 #ifndef MP_8BIT 919 a->dp[0] |= *b++; 920 a->used += 1; 921 #else 922 a->dp[0] = (*b & MP_MASK); 923 a->dp[1] |= ((*b++ >> 7U) & 1); 924 a->used += 2; 925 #endif 926 } 927 mp_clamp (a); 928 return MP_OKAY; 929 } 930 931 932 /* store in unsigned [big endian] format */ 933 static int mp_to_unsigned_bin (mp_int * a, unsigned char *b) 934 { 935 int x, res; 936 mp_int t; 937 938 if ((res = mp_init_copy (&t, a)) != MP_OKAY) { 939 return res; 940 } 941 942 x = 0; 943 while (mp_iszero (&t) == 0) { 944 #ifndef MP_8BIT 945 b[x++] = (unsigned char) (t.dp[0] & 255); 946 #else 947 b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7)); 948 #endif 949 if ((res = mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) { 950 mp_clear (&t); 951 return res; 952 } 953 } 954 bn_reverse (b, x); 955 mp_clear (&t); 956 return MP_OKAY; 957 } 958 959 960 /* shift right by a certain bit count (store quotient in c, optional remainder in d) */ 961 static int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d) 962 { 963 mp_digit D, r, rr; 964 int x, res; 965 mp_int t; 966 967 968 /* if the shift count is <= 0 then we do no work */ 969 if (b <= 0) { 970 res = mp_copy (a, c); 971 if (d != NULL) { 972 mp_zero (d); 973 } 974 return res; 975 } 976 977 if ((res = mp_init (&t)) != MP_OKAY) { 978 return res; 979 } 980 981 /* get the remainder */ 982 if (d != NULL) { 983 if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) { 984 mp_clear (&t); 985 return res; 986 } 987 } 988 989 /* copy */ 990 if ((res = mp_copy (a, c)) != MP_OKAY) { 991 mp_clear (&t); 992 return res; 993 } 994 995 /* shift by as many digits in the bit count */ 996 if (b >= (int)DIGIT_BIT) { 997 mp_rshd (c, b / DIGIT_BIT); 998 } 999 1000 /* shift any bit count < DIGIT_BIT */ 1001 D = (mp_digit) (b % DIGIT_BIT); 1002 if (D != 0) { 1003 register mp_digit *tmpc, mask, shift; 1004 1005 /* mask */ 1006 mask = (((mp_digit)1) << D) - 1; 1007 1008 /* shift for lsb */ 1009 shift = DIGIT_BIT - D; 1010 1011 /* alias */ 1012 tmpc = c->dp + (c->used - 1); 1013 1014 /* carry */ 1015 r = 0; 1016 for (x = c->used - 1; x >= 0; x--) { 1017 /* get the lower bits of this word in a temp */ 1018 rr = *tmpc & mask; 1019 1020 /* shift the current word and mix in the carry bits from the previous word */ 1021 *tmpc = (*tmpc >> D) | (r << shift); 1022 --tmpc; 1023 1024 /* set the carry to the carry bits of the current word found above */ 1025 r = rr; 1026 } 1027 } 1028 mp_clamp (c); 1029 if (d != NULL) { 1030 mp_exch (&t, d); 1031 } 1032 mp_clear (&t); 1033 return MP_OKAY; 1034 } 1035 1036 1037 static int mp_init_copy (mp_int * a, mp_int * b) 1038 { 1039 int res; 1040 1041 if ((res = mp_init (a)) != MP_OKAY) { 1042 return res; 1043 } 1044 return mp_copy (b, a); 1045 } 1046 1047 1048 /* set to zero */ 1049 static void mp_zero (mp_int * a) 1050 { 1051 int n; 1052 mp_digit *tmp; 1053 1054 a->sign = MP_ZPOS; 1055 a->used = 0; 1056 1057 tmp = a->dp; 1058 for (n = 0; n < a->alloc; n++) { 1059 *tmp++ = 0; 1060 } 1061 } 1062 1063 1064 /* copy, b = a */ 1065 static int mp_copy (mp_int * a, mp_int * b) 1066 { 1067 int res, n; 1068 1069 /* if dst == src do nothing */ 1070 if (a == b) { 1071 return MP_OKAY; 1072 } 1073 1074 /* grow dest */ 1075 if (b->alloc < a->used) { 1076 if ((res = mp_grow (b, a->used)) != MP_OKAY) { 1077 return res; 1078 } 1079 } 1080 1081 /* zero b and copy the parameters over */ 1082 { 1083 register mp_digit *tmpa, *tmpb; 1084 1085 /* pointer aliases */ 1086 1087 /* source */ 1088 tmpa = a->dp; 1089 1090 /* destination */ 1091 tmpb = b->dp; 1092 1093 /* copy all the digits */ 1094 for (n = 0; n < a->used; n++) { 1095 *tmpb++ = *tmpa++; 1096 } 1097 1098 /* clear high digits */ 1099 for (; n < b->used; n++) { 1100 *tmpb++ = 0; 1101 } 1102 } 1103 1104 /* copy used count and sign */ 1105 b->used = a->used; 1106 b->sign = a->sign; 1107 return MP_OKAY; 1108 } 1109 1110 1111 /* shift right a certain amount of digits */ 1112 static void mp_rshd (mp_int * a, int b) 1113 { 1114 int x; 1115 1116 /* if b <= 0 then ignore it */ 1117 if (b <= 0) { 1118 return; 1119 } 1120 1121 /* if b > used then simply zero it and return */ 1122 if (a->used <= b) { 1123 mp_zero (a); 1124 return; 1125 } 1126 1127 { 1128 register mp_digit *bottom, *top; 1129 1130 /* shift the digits down */ 1131 1132 /* bottom */ 1133 bottom = a->dp; 1134 1135 /* top [offset into digits] */ 1136 top = a->dp + b; 1137 1138 /* this is implemented as a sliding window where 1139 * the window is b-digits long and digits from 1140 * the top of the window are copied to the bottom 1141 * 1142 * e.g. 1143 1144 b-2 | b-1 | b0 | b1 | b2 | ... | bb | ----> 1145 /\ | ----> 1146 \-------------------/ ----> 1147 */ 1148 for (x = 0; x < (a->used - b); x++) { 1149 *bottom++ = *top++; 1150 } 1151 1152 /* zero the top digits */ 1153 for (; x < a->used; x++) { 1154 *bottom++ = 0; 1155 } 1156 } 1157 1158 /* remove excess digits */ 1159 a->used -= b; 1160 } 1161 1162 1163 /* swap the elements of two integers, for cases where you can't simply swap the 1164 * mp_int pointers around 1165 */ 1166 static void mp_exch (mp_int * a, mp_int * b) 1167 { 1168 mp_int t; 1169 1170 t = *a; 1171 *a = *b; 1172 *b = t; 1173 } 1174 1175 1176 /* trim unused digits 1177 * 1178 * This is used to ensure that leading zero digits are 1179 * trimed and the leading "used" digit will be non-zero 1180 * Typically very fast. Also fixes the sign if there 1181 * are no more leading digits 1182 */ 1183 static void mp_clamp (mp_int * a) 1184 { 1185 /* decrease used while the most significant digit is 1186 * zero. 1187 */ 1188 while (a->used > 0 && a->dp[a->used - 1] == 0) { 1189 --(a->used); 1190 } 1191 1192 /* reset the sign flag if used == 0 */ 1193 if (a->used == 0) { 1194 a->sign = MP_ZPOS; 1195 } 1196 } 1197 1198 1199 /* grow as required */ 1200 static int mp_grow (mp_int * a, int size) 1201 { 1202 int i; 1203 mp_digit *tmp; 1204 1205 /* if the alloc size is smaller alloc more ram */ 1206 if (a->alloc < size) { 1207 /* ensure there are always at least MP_PREC digits extra on top */ 1208 size += (MP_PREC * 2) - (size % MP_PREC); 1209 1210 /* reallocate the array a->dp 1211 * 1212 * We store the return in a temporary variable 1213 * in case the operation failed we don't want 1214 * to overwrite the dp member of a. 1215 */ 1216 tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * size); 1217 if (tmp == NULL) { 1218 /* reallocation failed but "a" is still valid [can be freed] */ 1219 return MP_MEM; 1220 } 1221 1222 /* reallocation succeeded so set a->dp */ 1223 a->dp = tmp; 1224 1225 /* zero excess digits */ 1226 i = a->alloc; 1227 a->alloc = size; 1228 for (; i < a->alloc; i++) { 1229 a->dp[i] = 0; 1230 } 1231 } 1232 return MP_OKAY; 1233 } 1234 1235 1236 /* b = |a| 1237 * 1238 * Simple function copies the input and fixes the sign to positive 1239 */ 1240 static int mp_abs (mp_int * a, mp_int * b) 1241 { 1242 int res; 1243 1244 /* copy a to b */ 1245 if (a != b) { 1246 if ((res = mp_copy (a, b)) != MP_OKAY) { 1247 return res; 1248 } 1249 } 1250 1251 /* force the sign of b to positive */ 1252 b->sign = MP_ZPOS; 1253 1254 return MP_OKAY; 1255 } 1256 1257 1258 /* set to a digit */ 1259 static void mp_set (mp_int * a, mp_digit b) 1260 { 1261 mp_zero (a); 1262 a->dp[0] = b & MP_MASK; 1263 a->used = (a->dp[0] != 0) ? 1 : 0; 1264 } 1265 1266 1267 /* b = a/2 */ 1268 static int mp_div_2(mp_int * a, mp_int * b) 1269 { 1270 int x, res, oldused; 1271 1272 /* copy */ 1273 if (b->alloc < a->used) { 1274 if ((res = mp_grow (b, a->used)) != MP_OKAY) { 1275 return res; 1276 } 1277 } 1278 1279 oldused = b->used; 1280 b->used = a->used; 1281 { 1282 register mp_digit r, rr, *tmpa, *tmpb; 1283 1284 /* source alias */ 1285 tmpa = a->dp + b->used - 1; 1286 1287 /* dest alias */ 1288 tmpb = b->dp + b->used - 1; 1289 1290 /* carry */ 1291 r = 0; 1292 for (x = b->used - 1; x >= 0; x--) { 1293 /* get the carry for the next iteration */ 1294 rr = *tmpa & 1; 1295 1296 /* shift the current digit, add in carry and store */ 1297 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1)); 1298 1299 /* forward carry to next iteration */ 1300 r = rr; 1301 } 1302 1303 /* zero excess digits */ 1304 tmpb = b->dp + b->used; 1305 for (x = b->used; x < oldused; x++) { 1306 *tmpb++ = 0; 1307 } 1308 } 1309 b->sign = a->sign; 1310 mp_clamp (b); 1311 return MP_OKAY; 1312 } 1313 1314 1315 /* shift left by a certain bit count */ 1316 static int mp_mul_2d (mp_int * a, int b, mp_int * c) 1317 { 1318 mp_digit d; 1319 int res; 1320 1321 /* copy */ 1322 if (a != c) { 1323 if ((res = mp_copy (a, c)) != MP_OKAY) { 1324 return res; 1325 } 1326 } 1327 1328 if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) { 1329 if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) { 1330 return res; 1331 } 1332 } 1333 1334 /* shift by as many digits in the bit count */ 1335 if (b >= (int)DIGIT_BIT) { 1336 if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) { 1337 return res; 1338 } 1339 } 1340 1341 /* shift any bit count < DIGIT_BIT */ 1342 d = (mp_digit) (b % DIGIT_BIT); 1343 if (d != 0) { 1344 register mp_digit *tmpc, shift, mask, r, rr; 1345 register int x; 1346 1347 /* bitmask for carries */ 1348 mask = (((mp_digit)1) << d) - 1; 1349 1350 /* shift for msbs */ 1351 shift = DIGIT_BIT - d; 1352 1353 /* alias */ 1354 tmpc = c->dp; 1355 1356 /* carry */ 1357 r = 0; 1358 for (x = 0; x < c->used; x++) { 1359 /* get the higher bits of the current word */ 1360 rr = (*tmpc >> shift) & mask; 1361 1362 /* shift the current word and OR in the carry */ 1363 *tmpc = ((*tmpc << d) | r) & MP_MASK; 1364 ++tmpc; 1365 1366 /* set the carry to the carry bits of the current word */ 1367 r = rr; 1368 } 1369 1370 /* set final carry */ 1371 if (r != 0) { 1372 c->dp[(c->used)++] = r; 1373 } 1374 } 1375 mp_clamp (c); 1376 return MP_OKAY; 1377 } 1378 1379 1380 static int mp_init_multi(mp_int *mp, ...) 1381 { 1382 mp_err res = MP_OKAY; /* Assume ok until proven otherwise */ 1383 int n = 0; /* Number of ok inits */ 1384 mp_int* cur_arg = mp; 1385 va_list args; 1386 1387 va_start(args, mp); /* init args to next argument from caller */ 1388 while (cur_arg != NULL) { 1389 if (mp_init(cur_arg) != MP_OKAY) { 1390 /* Oops - error! Back-track and mp_clear what we already 1391 succeeded in init-ing, then return error. 1392 */ 1393 va_list clean_args; 1394 1395 /* end the current list */ 1396 va_end(args); 1397 1398 /* now start cleaning up */ 1399 cur_arg = mp; 1400 va_start(clean_args, mp); 1401 while (n--) { 1402 mp_clear(cur_arg); 1403 cur_arg = va_arg(clean_args, mp_int*); 1404 } 1405 va_end(clean_args); 1406 res = MP_MEM; 1407 break; 1408 } 1409 n++; 1410 cur_arg = va_arg(args, mp_int*); 1411 } 1412 va_end(args); 1413 return res; /* Assumed ok, if error flagged above. */ 1414 } 1415 1416 1417 static void mp_clear_multi(mp_int *mp, ...) 1418 { 1419 mp_int* next_mp = mp; 1420 va_list args; 1421 va_start(args, mp); 1422 while (next_mp != NULL) { 1423 mp_clear(next_mp); 1424 next_mp = va_arg(args, mp_int*); 1425 } 1426 va_end(args); 1427 } 1428 1429 1430 /* shift left a certain amount of digits */ 1431 static int mp_lshd (mp_int * a, int b) 1432 { 1433 int x, res; 1434 1435 /* if its less than zero return */ 1436 if (b <= 0) { 1437 return MP_OKAY; 1438 } 1439 1440 /* grow to fit the new digits */ 1441 if (a->alloc < a->used + b) { 1442 if ((res = mp_grow (a, a->used + b)) != MP_OKAY) { 1443 return res; 1444 } 1445 } 1446 1447 { 1448 register mp_digit *top, *bottom; 1449 1450 /* increment the used by the shift amount then copy upwards */ 1451 a->used += b; 1452 1453 /* top */ 1454 top = a->dp + a->used - 1; 1455 1456 /* base */ 1457 bottom = a->dp + a->used - 1 - b; 1458 1459 /* much like mp_rshd this is implemented using a sliding window 1460 * except the window goes the otherway around. Copying from 1461 * the bottom to the top. see bn_mp_rshd.c for more info. 1462 */ 1463 for (x = a->used - 1; x >= b; x--) { 1464 *top-- = *bottom--; 1465 } 1466 1467 /* zero the lower digits */ 1468 top = a->dp; 1469 for (x = 0; x < b; x++) { 1470 *top++ = 0; 1471 } 1472 } 1473 return MP_OKAY; 1474 } 1475 1476 1477 /* returns the number of bits in an int */ 1478 static int mp_count_bits (mp_int * a) 1479 { 1480 int r; 1481 mp_digit q; 1482 1483 /* shortcut */ 1484 if (a->used == 0) { 1485 return 0; 1486 } 1487 1488 /* get number of digits and add that */ 1489 r = (a->used - 1) * DIGIT_BIT; 1490 1491 /* take the last digit and count the bits in it */ 1492 q = a->dp[a->used - 1]; 1493 while (q > ((mp_digit) 0)) { 1494 ++r; 1495 q >>= ((mp_digit) 1); 1496 } 1497 return r; 1498 } 1499 1500 1501 /* calc a value mod 2**b */ 1502 static int mp_mod_2d (mp_int * a, int b, mp_int * c) 1503 { 1504 int x, res; 1505 1506 /* if b is <= 0 then zero the int */ 1507 if (b <= 0) { 1508 mp_zero (c); 1509 return MP_OKAY; 1510 } 1511 1512 /* if the modulus is larger than the value than return */ 1513 if (b >= (int) (a->used * DIGIT_BIT)) { 1514 res = mp_copy (a, c); 1515 return res; 1516 } 1517 1518 /* copy */ 1519 if ((res = mp_copy (a, c)) != MP_OKAY) { 1520 return res; 1521 } 1522 1523 /* zero digits above the last digit of the modulus */ 1524 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) { 1525 c->dp[x] = 0; 1526 } 1527 /* clear the digit that is not completely outside/inside the modulus */ 1528 c->dp[b / DIGIT_BIT] &= 1529 (mp_digit) ((((mp_digit) 1) << (((mp_digit) b) % DIGIT_BIT)) - ((mp_digit) 1)); 1530 mp_clamp (c); 1531 return MP_OKAY; 1532 } 1533 1534 1535 /* slower bit-bang division... also smaller */ 1536 static int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d) 1537 { 1538 mp_int ta, tb, tq, q; 1539 int res, n, n2; 1540 1541 /* is divisor zero ? */ 1542 if (mp_iszero (b) == 1) { 1543 return MP_VAL; 1544 } 1545 1546 /* if a < b then q=0, r = a */ 1547 if (mp_cmp_mag (a, b) == MP_LT) { 1548 if (d != NULL) { 1549 res = mp_copy (a, d); 1550 } else { 1551 res = MP_OKAY; 1552 } 1553 if (c != NULL) { 1554 mp_zero (c); 1555 } 1556 return res; 1557 } 1558 1559 /* init our temps */ 1560 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) { 1561 return res; 1562 } 1563 1564 1565 mp_set(&tq, 1); 1566 n = mp_count_bits(a) - mp_count_bits(b); 1567 if (((res = mp_abs(a, &ta)) != MP_OKAY) || 1568 ((res = mp_abs(b, &tb)) != MP_OKAY) || 1569 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || 1570 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { 1571 goto LBL_ERR; 1572 } 1573 1574 while (n-- >= 0) { 1575 if (mp_cmp(&tb, &ta) != MP_GT) { 1576 if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || 1577 ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { 1578 goto LBL_ERR; 1579 } 1580 } 1581 if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) || 1582 ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) { 1583 goto LBL_ERR; 1584 } 1585 } 1586 1587 /* now q == quotient and ta == remainder */ 1588 n = a->sign; 1589 n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG); 1590 if (c != NULL) { 1591 mp_exch(c, &q); 1592 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; 1593 } 1594 if (d != NULL) { 1595 mp_exch(d, &ta); 1596 d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; 1597 } 1598 LBL_ERR: 1599 mp_clear_multi(&ta, &tb, &tq, &q, NULL); 1600 return res; 1601 } 1602 1603 1604 #ifdef MP_LOW_MEM 1605 #define TAB_SIZE 32 1606 #else 1607 #define TAB_SIZE 256 1608 #endif 1609 1610 static int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode) 1611 { 1612 mp_int M[TAB_SIZE], res, mu; 1613 mp_digit buf; 1614 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; 1615 int (*redux)(mp_int*,mp_int*,mp_int*); 1616 1617 /* find window size */ 1618 x = mp_count_bits (X); 1619 if (x <= 7) { 1620 winsize = 2; 1621 } else if (x <= 36) { 1622 winsize = 3; 1623 } else if (x <= 140) { 1624 winsize = 4; 1625 } else if (x <= 450) { 1626 winsize = 5; 1627 } else if (x <= 1303) { 1628 winsize = 6; 1629 } else if (x <= 3529) { 1630 winsize = 7; 1631 } else { 1632 winsize = 8; 1633 } 1634 1635 #ifdef MP_LOW_MEM 1636 if (winsize > 5) { 1637 winsize = 5; 1638 } 1639 #endif 1640 1641 /* init M array */ 1642 /* init first cell */ 1643 if ((err = mp_init(&M[1])) != MP_OKAY) { 1644 return err; 1645 } 1646 1647 /* now init the second half of the array */ 1648 for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 1649 if ((err = mp_init(&M[x])) != MP_OKAY) { 1650 for (y = 1<<(winsize-1); y < x; y++) { 1651 mp_clear (&M[y]); 1652 } 1653 mp_clear(&M[1]); 1654 return err; 1655 } 1656 } 1657 1658 /* create mu, used for Barrett reduction */ 1659 if ((err = mp_init (&mu)) != MP_OKAY) { 1660 goto LBL_M; 1661 } 1662 1663 if (redmode == 0) { 1664 if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) { 1665 goto LBL_MU; 1666 } 1667 redux = mp_reduce; 1668 } else { 1669 if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) { 1670 goto LBL_MU; 1671 } 1672 redux = mp_reduce_2k_l; 1673 } 1674 1675 /* create M table 1676 * 1677 * The M table contains powers of the base, 1678 * e.g. M[x] = G**x mod P 1679 * 1680 * The first half of the table is not 1681 * computed though accept for M[0] and M[1] 1682 */ 1683 if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) { 1684 goto LBL_MU; 1685 } 1686 1687 /* compute the value at M[1<<(winsize-1)] by squaring 1688 * M[1] (winsize-1) times 1689 */ 1690 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) { 1691 goto LBL_MU; 1692 } 1693 1694 for (x = 0; x < (winsize - 1); x++) { 1695 /* square it */ 1696 if ((err = mp_sqr (&M[1 << (winsize - 1)], 1697 &M[1 << (winsize - 1)])) != MP_OKAY) { 1698 goto LBL_MU; 1699 } 1700 1701 /* reduce modulo P */ 1702 if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) { 1703 goto LBL_MU; 1704 } 1705 } 1706 1707 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P) 1708 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1) 1709 */ 1710 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) { 1711 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) { 1712 goto LBL_MU; 1713 } 1714 if ((err = redux (&M[x], P, &mu)) != MP_OKAY) { 1715 goto LBL_MU; 1716 } 1717 } 1718 1719 /* setup result */ 1720 if ((err = mp_init (&res)) != MP_OKAY) { 1721 goto LBL_MU; 1722 } 1723 mp_set (&res, 1); 1724 1725 /* set initial mode and bit cnt */ 1726 mode = 0; 1727 bitcnt = 1; 1728 buf = 0; 1729 digidx = X->used - 1; 1730 bitcpy = 0; 1731 bitbuf = 0; 1732 1733 for (;;) { 1734 /* grab next digit as required */ 1735 if (--bitcnt == 0) { 1736 /* if digidx == -1 we are out of digits */ 1737 if (digidx == -1) { 1738 break; 1739 } 1740 /* read next digit and reset the bitcnt */ 1741 buf = X->dp[digidx--]; 1742 bitcnt = (int) DIGIT_BIT; 1743 } 1744 1745 /* grab the next msb from the exponent */ 1746 y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1; 1747 buf <<= (mp_digit)1; 1748 1749 /* if the bit is zero and mode == 0 then we ignore it 1750 * These represent the leading zero bits before the first 1 bit 1751 * in the exponent. Technically this opt is not required but it 1752 * does lower the # of trivial squaring/reductions used 1753 */ 1754 if (mode == 0 && y == 0) { 1755 continue; 1756 } 1757 1758 /* if the bit is zero and mode == 1 then we square */ 1759 if (mode == 1 && y == 0) { 1760 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 1761 goto LBL_RES; 1762 } 1763 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 1764 goto LBL_RES; 1765 } 1766 continue; 1767 } 1768 1769 /* else we add it to the window */ 1770 bitbuf |= (y << (winsize - ++bitcpy)); 1771 mode = 2; 1772 1773 if (bitcpy == winsize) { 1774 /* ok window is filled so square as required and multiply */ 1775 /* square first */ 1776 for (x = 0; x < winsize; x++) { 1777 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 1778 goto LBL_RES; 1779 } 1780 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 1781 goto LBL_RES; 1782 } 1783 } 1784 1785 /* then multiply */ 1786 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) { 1787 goto LBL_RES; 1788 } 1789 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 1790 goto LBL_RES; 1791 } 1792 1793 /* empty window and reset */ 1794 bitcpy = 0; 1795 bitbuf = 0; 1796 mode = 1; 1797 } 1798 } 1799 1800 /* if bits remain then square/multiply */ 1801 if (mode == 2 && bitcpy > 0) { 1802 /* square then multiply if the bit is set */ 1803 for (x = 0; x < bitcpy; x++) { 1804 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 1805 goto LBL_RES; 1806 } 1807 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 1808 goto LBL_RES; 1809 } 1810 1811 bitbuf <<= 1; 1812 if ((bitbuf & (1 << winsize)) != 0) { 1813 /* then multiply */ 1814 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) { 1815 goto LBL_RES; 1816 } 1817 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 1818 goto LBL_RES; 1819 } 1820 } 1821 } 1822 } 1823 1824 mp_exch (&res, Y); 1825 err = MP_OKAY; 1826 LBL_RES:mp_clear (&res); 1827 LBL_MU:mp_clear (&mu); 1828 LBL_M: 1829 mp_clear(&M[1]); 1830 for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 1831 mp_clear (&M[x]); 1832 } 1833 return err; 1834 } 1835 1836 1837 /* computes b = a*a */ 1838 static int mp_sqr (mp_int * a, mp_int * b) 1839 { 1840 int res; 1841 1842 #ifdef BN_MP_TOOM_SQR_C 1843 /* use Toom-Cook? */ 1844 if (a->used >= TOOM_SQR_CUTOFF) { 1845 res = mp_toom_sqr(a, b); 1846 /* Karatsuba? */ 1847 } else 1848 #endif 1849 #ifdef BN_MP_KARATSUBA_SQR_C 1850 if (a->used >= KARATSUBA_SQR_CUTOFF) { 1851 res = mp_karatsuba_sqr (a, b); 1852 } else 1853 #endif 1854 { 1855 #ifdef BN_FAST_S_MP_SQR_C 1856 /* can we use the fast comba multiplier? */ 1857 if ((a->used * 2 + 1) < MP_WARRAY && 1858 a->used < 1859 (1 << (sizeof(mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) { 1860 res = fast_s_mp_sqr (a, b); 1861 } else 1862 #endif 1863 #ifdef BN_S_MP_SQR_C 1864 res = s_mp_sqr (a, b); 1865 #else 1866 #error mp_sqr could fail 1867 res = MP_VAL; 1868 #endif 1869 } 1870 b->sign = MP_ZPOS; 1871 return res; 1872 } 1873 1874 1875 /* reduces a modulo n where n is of the form 2**p - d 1876 This differs from reduce_2k since "d" can be larger 1877 than a single digit. 1878 */ 1879 static int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d) 1880 { 1881 mp_int q; 1882 int p, res; 1883 1884 if ((res = mp_init(&q)) != MP_OKAY) { 1885 return res; 1886 } 1887 1888 p = mp_count_bits(n); 1889 top: 1890 /* q = a/2**p, a = a mod 2**p */ 1891 if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) { 1892 goto ERR; 1893 } 1894 1895 /* q = q * d */ 1896 if ((res = mp_mul(&q, d, &q)) != MP_OKAY) { 1897 goto ERR; 1898 } 1899 1900 /* a = a + q */ 1901 if ((res = s_mp_add(a, &q, a)) != MP_OKAY) { 1902 goto ERR; 1903 } 1904 1905 if (mp_cmp_mag(a, n) != MP_LT) { 1906 s_mp_sub(a, n, a); 1907 goto top; 1908 } 1909 1910 ERR: 1911 mp_clear(&q); 1912 return res; 1913 } 1914 1915 1916 /* determines the setup value */ 1917 static int mp_reduce_2k_setup_l(mp_int *a, mp_int *d) 1918 { 1919 int res; 1920 mp_int tmp; 1921 1922 if ((res = mp_init(&tmp)) != MP_OKAY) { 1923 return res; 1924 } 1925 1926 if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) { 1927 goto ERR; 1928 } 1929 1930 if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) { 1931 goto ERR; 1932 } 1933 1934 ERR: 1935 mp_clear(&tmp); 1936 return res; 1937 } 1938 1939 1940 /* computes a = 2**b 1941 * 1942 * Simple algorithm which zeroes the int, grows it then just sets one bit 1943 * as required. 1944 */ 1945 static int mp_2expt (mp_int * a, int b) 1946 { 1947 int res; 1948 1949 /* zero a as per default */ 1950 mp_zero (a); 1951 1952 /* grow a to accomodate the single bit */ 1953 if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) { 1954 return res; 1955 } 1956 1957 /* set the used count of where the bit will go */ 1958 a->used = b / DIGIT_BIT + 1; 1959 1960 /* put the single bit in its place */ 1961 a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT); 1962 1963 return MP_OKAY; 1964 } 1965 1966 1967 /* pre-calculate the value required for Barrett reduction 1968 * For a given modulus "b" it calulates the value required in "a" 1969 */ 1970 static int mp_reduce_setup (mp_int * a, mp_int * b) 1971 { 1972 int res; 1973 1974 if ((res = mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) { 1975 return res; 1976 } 1977 return mp_div (a, b, a, NULL); 1978 } 1979 1980 1981 /* reduces x mod m, assumes 0 < x < m**2, mu is 1982 * precomputed via mp_reduce_setup. 1983 * From HAC pp.604 Algorithm 14.42 1984 */ 1985 static int mp_reduce (mp_int * x, mp_int * m, mp_int * mu) 1986 { 1987 mp_int q; 1988 int res, um = m->used; 1989 1990 /* q = x */ 1991 if ((res = mp_init_copy (&q, x)) != MP_OKAY) { 1992 return res; 1993 } 1994 1995 /* q1 = x / b**(k-1) */ 1996 mp_rshd (&q, um - 1); 1997 1998 /* according to HAC this optimization is ok */ 1999 if (((unsigned long) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) { 2000 if ((res = mp_mul (&q, mu, &q)) != MP_OKAY) { 2001 goto CLEANUP; 2002 } 2003 } else { 2004 #ifdef BN_S_MP_MUL_HIGH_DIGS_C 2005 if ((res = s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) { 2006 goto CLEANUP; 2007 } 2008 #elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C) 2009 if ((res = fast_s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) { 2010 goto CLEANUP; 2011 } 2012 #else 2013 { 2014 #error mp_reduce would always fail 2015 res = MP_VAL; 2016 goto CLEANUP; 2017 } 2018 #endif 2019 } 2020 2021 /* q3 = q2 / b**(k+1) */ 2022 mp_rshd (&q, um + 1); 2023 2024 /* x = x mod b**(k+1), quick (no division) */ 2025 if ((res = mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) { 2026 goto CLEANUP; 2027 } 2028 2029 /* q = q * m mod b**(k+1), quick (no division) */ 2030 if ((res = s_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) { 2031 goto CLEANUP; 2032 } 2033 2034 /* x = x - q */ 2035 if ((res = mp_sub (x, &q, x)) != MP_OKAY) { 2036 goto CLEANUP; 2037 } 2038 2039 /* If x < 0, add b**(k+1) to it */ 2040 if (mp_cmp_d (x, 0) == MP_LT) { 2041 mp_set (&q, 1); 2042 if ((res = mp_lshd (&q, um + 1)) != MP_OKAY) { 2043 goto CLEANUP; 2044 } 2045 if ((res = mp_add (x, &q, x)) != MP_OKAY) { 2046 goto CLEANUP; 2047 } 2048 } 2049 2050 /* Back off if it's too big */ 2051 while (mp_cmp (x, m) != MP_LT) { 2052 if ((res = s_mp_sub (x, m, x)) != MP_OKAY) { 2053 goto CLEANUP; 2054 } 2055 } 2056 2057 CLEANUP: 2058 mp_clear (&q); 2059 2060 return res; 2061 } 2062 2063 2064 /* multiplies |a| * |b| and only computes upto digs digits of result 2065 * HAC pp. 595, Algorithm 14.12 Modified so you can control how 2066 * many digits of output are created. 2067 */ 2068 static int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs) 2069 { 2070 mp_int t; 2071 int res, pa, pb, ix, iy; 2072 mp_digit u; 2073 mp_word r; 2074 mp_digit tmpx, *tmpt, *tmpy; 2075 2076 /* can we use the fast multiplier? */ 2077 if (((digs) < MP_WARRAY) && 2078 MIN (a->used, b->used) < 2079 (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { 2080 return fast_s_mp_mul_digs (a, b, c, digs); 2081 } 2082 2083 if ((res = mp_init_size (&t, digs)) != MP_OKAY) { 2084 return res; 2085 } 2086 t.used = digs; 2087 2088 /* compute the digits of the product directly */ 2089 pa = a->used; 2090 for (ix = 0; ix < pa; ix++) { 2091 /* set the carry to zero */ 2092 u = 0; 2093 2094 /* limit ourselves to making digs digits of output */ 2095 pb = MIN (b->used, digs - ix); 2096 2097 /* setup some aliases */ 2098 /* copy of the digit from a used within the nested loop */ 2099 tmpx = a->dp[ix]; 2100 2101 /* an alias for the destination shifted ix places */ 2102 tmpt = t.dp + ix; 2103 2104 /* an alias for the digits of b */ 2105 tmpy = b->dp; 2106 2107 /* compute the columns of the output and propagate the carry */ 2108 for (iy = 0; iy < pb; iy++) { 2109 /* compute the column as a mp_word */ 2110 r = ((mp_word)*tmpt) + 2111 ((mp_word)tmpx) * ((mp_word)*tmpy++) + 2112 ((mp_word) u); 2113 2114 /* the new column is the lower part of the result */ 2115 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2116 2117 /* get the carry word from the result */ 2118 u = (mp_digit) (r >> ((mp_word) DIGIT_BIT)); 2119 } 2120 /* set carry if it is placed below digs */ 2121 if (ix + iy < digs) { 2122 *tmpt = u; 2123 } 2124 } 2125 2126 mp_clamp (&t); 2127 mp_exch (&t, c); 2128 2129 mp_clear (&t); 2130 return MP_OKAY; 2131 } 2132 2133 2134 /* Fast (comba) multiplier 2135 * 2136 * This is the fast column-array [comba] multiplier. It is 2137 * designed to compute the columns of the product first 2138 * then handle the carries afterwards. This has the effect 2139 * of making the nested loops that compute the columns very 2140 * simple and schedulable on super-scalar processors. 2141 * 2142 * This has been modified to produce a variable number of 2143 * digits of output so if say only a half-product is required 2144 * you don't have to compute the upper half (a feature 2145 * required for fast Barrett reduction). 2146 * 2147 * Based on Algorithm 14.12 on pp.595 of HAC. 2148 * 2149 */ 2150 static int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs) 2151 { 2152 int olduse, res, pa, ix, iz; 2153 mp_digit W[MP_WARRAY]; 2154 register mp_word _W; 2155 2156 /* grow the destination as required */ 2157 if (c->alloc < digs) { 2158 if ((res = mp_grow (c, digs)) != MP_OKAY) { 2159 return res; 2160 } 2161 } 2162 2163 /* number of output digits to produce */ 2164 pa = MIN(digs, a->used + b->used); 2165 2166 /* clear the carry */ 2167 _W = 0; 2168 for (ix = 0; ix < pa; ix++) { 2169 int tx, ty; 2170 int iy; 2171 mp_digit *tmpx, *tmpy; 2172 2173 /* get offsets into the two bignums */ 2174 ty = MIN(b->used-1, ix); 2175 tx = ix - ty; 2176 2177 /* setup temp aliases */ 2178 tmpx = a->dp + tx; 2179 tmpy = b->dp + ty; 2180 2181 /* this is the number of times the loop will iterrate, essentially 2182 while (tx++ < a->used && ty-- >= 0) { ... } 2183 */ 2184 iy = MIN(a->used-tx, ty+1); 2185 2186 /* execute loop */ 2187 for (iz = 0; iz < iy; ++iz) { 2188 _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--); 2189 2190 } 2191 2192 /* store term */ 2193 W[ix] = ((mp_digit)_W) & MP_MASK; 2194 2195 /* make next carry */ 2196 _W = _W >> ((mp_word)DIGIT_BIT); 2197 } 2198 2199 /* setup dest */ 2200 olduse = c->used; 2201 c->used = pa; 2202 2203 { 2204 register mp_digit *tmpc; 2205 tmpc = c->dp; 2206 for (ix = 0; ix < pa+1; ix++) { 2207 /* now extract the previous digit [below the carry] */ 2208 *tmpc++ = W[ix]; 2209 } 2210 2211 /* clear unused digits [that existed in the old copy of c] */ 2212 for (; ix < olduse; ix++) { 2213 *tmpc++ = 0; 2214 } 2215 } 2216 mp_clamp (c); 2217 return MP_OKAY; 2218 } 2219 2220 2221 /* init an mp_init for a given size */ 2222 static int mp_init_size (mp_int * a, int size) 2223 { 2224 int x; 2225 2226 /* pad size so there are always extra digits */ 2227 size += (MP_PREC * 2) - (size % MP_PREC); 2228 2229 /* alloc mem */ 2230 a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * size); 2231 if (a->dp == NULL) { 2232 return MP_MEM; 2233 } 2234 2235 /* set the members */ 2236 a->used = 0; 2237 a->alloc = size; 2238 a->sign = MP_ZPOS; 2239 2240 /* zero the digits */ 2241 for (x = 0; x < size; x++) { 2242 a->dp[x] = 0; 2243 } 2244 2245 return MP_OKAY; 2246 } 2247 2248 2249 /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */ 2250 static int s_mp_sqr (mp_int * a, mp_int * b) 2251 { 2252 mp_int t; 2253 int res, ix, iy, pa; 2254 mp_word r; 2255 mp_digit u, tmpx, *tmpt; 2256 2257 pa = a->used; 2258 if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) { 2259 return res; 2260 } 2261 2262 /* default used is maximum possible size */ 2263 t.used = 2*pa + 1; 2264 2265 for (ix = 0; ix < pa; ix++) { 2266 /* first calculate the digit at 2*ix */ 2267 /* calculate double precision result */ 2268 r = ((mp_word) t.dp[2*ix]) + 2269 ((mp_word)a->dp[ix])*((mp_word)a->dp[ix]); 2270 2271 /* store lower part in result */ 2272 t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK)); 2273 2274 /* get the carry */ 2275 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT)); 2276 2277 /* left hand side of A[ix] * A[iy] */ 2278 tmpx = a->dp[ix]; 2279 2280 /* alias for where to store the results */ 2281 tmpt = t.dp + (2*ix + 1); 2282 2283 for (iy = ix + 1; iy < pa; iy++) { 2284 /* first calculate the product */ 2285 r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]); 2286 2287 /* now calculate the double precision result, note we use 2288 * addition instead of *2 since it's easier to optimize 2289 */ 2290 r = ((mp_word) *tmpt) + r + r + ((mp_word) u); 2291 2292 /* store lower part */ 2293 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2294 2295 /* get carry */ 2296 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT)); 2297 } 2298 /* propagate upwards */ 2299 while (u != ((mp_digit) 0)) { 2300 r = ((mp_word) *tmpt) + ((mp_word) u); 2301 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2302 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT)); 2303 } 2304 } 2305 2306 mp_clamp (&t); 2307 mp_exch (&t, b); 2308 mp_clear (&t); 2309 return MP_OKAY; 2310 } 2311 2312 2313 /* multiplies |a| * |b| and does not compute the lower digs digits 2314 * [meant to get the higher part of the product] 2315 */ 2316 static int s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs) 2317 { 2318 mp_int t; 2319 int res, pa, pb, ix, iy; 2320 mp_digit u; 2321 mp_word r; 2322 mp_digit tmpx, *tmpt, *tmpy; 2323 2324 /* can we use the fast multiplier? */ 2325 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C 2326 if (((a->used + b->used + 1) < MP_WARRAY) 2327 && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { 2328 return fast_s_mp_mul_high_digs (a, b, c, digs); 2329 } 2330 #endif 2331 2332 if ((res = mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) { 2333 return res; 2334 } 2335 t.used = a->used + b->used + 1; 2336 2337 pa = a->used; 2338 pb = b->used; 2339 for (ix = 0; ix < pa; ix++) { 2340 /* clear the carry */ 2341 u = 0; 2342 2343 /* left hand side of A[ix] * B[iy] */ 2344 tmpx = a->dp[ix]; 2345 2346 /* alias to the address of where the digits will be stored */ 2347 tmpt = &(t.dp[digs]); 2348 2349 /* alias for where to read the right hand side from */ 2350 tmpy = b->dp + (digs - ix); 2351 2352 for (iy = digs - ix; iy < pb; iy++) { 2353 /* calculate the double precision result */ 2354 r = ((mp_word)*tmpt) + 2355 ((mp_word)tmpx) * ((mp_word)*tmpy++) + 2356 ((mp_word) u); 2357 2358 /* get the lower part */ 2359 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2360 2361 /* carry the carry */ 2362 u = (mp_digit) (r >> ((mp_word) DIGIT_BIT)); 2363 } 2364 *tmpt = u; 2365 } 2366 mp_clamp (&t); 2367 mp_exch (&t, c); 2368 mp_clear (&t); 2369 return MP_OKAY; 2370 } 2371
