1 /* 2 * mpi.c 3 * 4 * Arbitrary precision integer arithmetic library 5 * 6 * ***** BEGIN LICENSE BLOCK ***** 7 * Version: MPL 1.1/GPL 2.0/LGPL 2.1 8 * 9 * The contents of this file are subject to the Mozilla Public License Version 10 * 1.1 (the "License"); you may not use this file except in compliance with 11 * the License. You may obtain a copy of the License at 12 * http://www.mozilla.org/MPL/ 13 * 14 * Software distributed under the License is distributed on an "AS IS" basis, 15 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 16 * for the specific language governing rights and limitations under the 17 * License. 18 * 19 * The Original Code is the MPI Arbitrary Precision Integer Arithmetic library. 20 * 21 * The Initial Developer of the Original Code is 22 * Michael J. Fromberger. 23 * Portions created by the Initial Developer are Copyright (C) 1998 24 * the Initial Developer. All Rights Reserved. 25 * 26 * Contributor(s): 27 * Netscape Communications Corporation 28 * Douglas Stebila <douglas (at) stebila.ca> of Sun Laboratories. 29 * 30 * Alternatively, the contents of this file may be used under the terms of 31 * either the GNU General Public License Version 2 or later (the "GPL"), or 32 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 33 * in which case the provisions of the GPL or the LGPL are applicable instead 34 * of those above. If you wish to allow use of your version of this file only 35 * under the terms of either the GPL or the LGPL, and not to allow others to 36 * use your version of this file under the terms of the MPL, indicate your 37 * decision by deleting the provisions above and replace them with the notice 38 * and other provisions required by the GPL or the LGPL. If you do not delete 39 * the provisions above, a recipient may use your version of this file under 40 * the terms of any one of the MPL, the GPL or the LGPL. 41 * 42 * ***** END LICENSE BLOCK ***** */ 43 /* 44 * Copyright 2007 Sun Microsystems, Inc. All rights reserved. 45 * Use is subject to license terms. 46 * 47 * Sun elects to use this software under the MPL license. 48 */ 49 50 #pragma ident "%Z%%M% %I% %E% SMI" 51 52 /* $Id: mpi.c,v 1.45 2006/09/29 20:12:21 alexei.volkov.bugs%sun.com Exp $ */ 53 54 #include "mpi-priv.h" 55 #if defined(OSF1) 56 #include <c_asm.h> 57 #endif 58 59 #if MP_LOGTAB 60 /* 61 A table of the logs of 2 for various bases (the 0 and 1 entries of 62 this table are meaningless and should not be referenced). 63 64 This table is used to compute output lengths for the mp_toradix() 65 function. Since a number n in radix r takes up about log_r(n) 66 digits, we estimate the output size by taking the least integer 67 greater than log_r(n), where: 68 69 log_r(n) = log_2(n) * log_r(2) 70 71 This table, therefore, is a table of log_r(2) for 2 <= r <= 36, 72 which are the output bases supported. 73 */ 74 #include "logtab.h" 75 #endif 76 77 /* {{{ Constant strings */ 78 79 /* Constant strings returned by mp_strerror() */ 80 static const char *mp_err_string[] = { 81 "unknown result code", /* say what? */ 82 "boolean true", /* MP_OKAY, MP_YES */ 83 "boolean false", /* MP_NO */ 84 "out of memory", /* MP_MEM */ 85 "argument out of range", /* MP_RANGE */ 86 "invalid input parameter", /* MP_BADARG */ 87 "result is undefined" /* MP_UNDEF */ 88 }; 89 90 /* Value to digit maps for radix conversion */ 91 92 /* s_dmap_1 - standard digits and letters */ 93 static const char *s_dmap_1 = 94 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/"; 95 96 /* }}} */ 97 98 unsigned long mp_allocs; 99 unsigned long mp_frees; 100 unsigned long mp_copies; 101 102 /* {{{ Default precision manipulation */ 103 104 /* Default precision for newly created mp_int's */ 105 static mp_size s_mp_defprec = MP_DEFPREC; 106 107 mp_size mp_get_prec(void) 108 { 109 return s_mp_defprec; 110 111 } /* end mp_get_prec() */ 112 113 void mp_set_prec(mp_size prec) 114 { 115 if(prec == 0) 116 s_mp_defprec = MP_DEFPREC; 117 else 118 s_mp_defprec = prec; 119 120 } /* end mp_set_prec() */ 121 122 /* }}} */ 123 124 /*------------------------------------------------------------------------*/ 125 /* {{{ mp_init(mp, kmflag) */ 126 127 /* 128 mp_init(mp, kmflag) 129 130 Initialize a new zero-valued mp_int. Returns MP_OKAY if successful, 131 MP_MEM if memory could not be allocated for the structure. 132 */ 133 134 mp_err mp_init(mp_int *mp, int kmflag) 135 { 136 return mp_init_size(mp, s_mp_defprec, kmflag); 137 138 } /* end mp_init() */ 139 140 /* }}} */ 141 142 /* {{{ mp_init_size(mp, prec, kmflag) */ 143 144 /* 145 mp_init_size(mp, prec, kmflag) 146 147 Initialize a new zero-valued mp_int with at least the given 148 precision; returns MP_OKAY if successful, or MP_MEM if memory could 149 not be allocated for the structure. 150 */ 151 152 mp_err mp_init_size(mp_int *mp, mp_size prec, int kmflag) 153 { 154 ARGCHK(mp != NULL && prec > 0, MP_BADARG); 155 156 prec = MP_ROUNDUP(prec, s_mp_defprec); 157 if((DIGITS(mp) = s_mp_alloc(prec, sizeof(mp_digit), kmflag)) == NULL) 158 return MP_MEM; 159 160 SIGN(mp) = ZPOS; 161 USED(mp) = 1; 162 ALLOC(mp) = prec; 163 164 return MP_OKAY; 165 166 } /* end mp_init_size() */ 167 168 /* }}} */ 169 170 /* {{{ mp_init_copy(mp, from) */ 171 172 /* 173 mp_init_copy(mp, from) 174 175 Initialize mp as an exact copy of from. Returns MP_OKAY if 176 successful, MP_MEM if memory could not be allocated for the new 177 structure. 178 */ 179 180 mp_err mp_init_copy(mp_int *mp, const mp_int *from) 181 { 182 ARGCHK(mp != NULL && from != NULL, MP_BADARG); 183 184 if(mp == from) 185 return MP_OKAY; 186 187 if((DIGITS(mp) = s_mp_alloc(ALLOC(from), sizeof(mp_digit), FLAG(from))) == NULL) 188 return MP_MEM; 189 190 s_mp_copy(DIGITS(from), DIGITS(mp), USED(from)); 191 USED(mp) = USED(from); 192 ALLOC(mp) = ALLOC(from); 193 SIGN(mp) = SIGN(from); 194 FLAG(mp) = FLAG(from); 195 196 return MP_OKAY; 197 198 } /* end mp_init_copy() */ 199 200 /* }}} */ 201 202 /* {{{ mp_copy(from, to) */ 203 204 /* 205 mp_copy(from, to) 206 207 Copies the mp_int 'from' to the mp_int 'to'. It is presumed that 208 'to' has already been initialized (if not, use mp_init_copy() 209 instead). If 'from' and 'to' are identical, nothing happens. 210 */ 211 212 mp_err mp_copy(const mp_int *from, mp_int *to) 213 { 214 ARGCHK(from != NULL && to != NULL, MP_BADARG); 215 216 if(from == to) 217 return MP_OKAY; 218 219 ++mp_copies; 220 { /* copy */ 221 mp_digit *tmp; 222 223 /* 224 If the allocated buffer in 'to' already has enough space to hold 225 all the used digits of 'from', we'll re-use it to avoid hitting 226 the memory allocater more than necessary; otherwise, we'd have 227 to grow anyway, so we just allocate a hunk and make the copy as 228 usual 229 */ 230 if(ALLOC(to) >= USED(from)) { 231 s_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from)); 232 s_mp_copy(DIGITS(from), DIGITS(to), USED(from)); 233 234 } else { 235 if((tmp = s_mp_alloc(ALLOC(from), sizeof(mp_digit), FLAG(from))) == NULL) 236 return MP_MEM; 237 238 s_mp_copy(DIGITS(from), tmp, USED(from)); 239 240 if(DIGITS(to) != NULL) { 241 #if MP_CRYPTO 242 s_mp_setz(DIGITS(to), ALLOC(to)); 243 #endif 244 s_mp_free(DIGITS(to), ALLOC(to)); 245 } 246 247 DIGITS(to) = tmp; 248 ALLOC(to) = ALLOC(from); 249 } 250 251 /* Copy the precision and sign from the original */ 252 USED(to) = USED(from); 253 SIGN(to) = SIGN(from); 254 } /* end copy */ 255 256 return MP_OKAY; 257 258 } /* end mp_copy() */ 259 260 /* }}} */ 261 262 /* {{{ mp_exch(mp1, mp2) */ 263 264 /* 265 mp_exch(mp1, mp2) 266 267 Exchange mp1 and mp2 without allocating any intermediate memory 268 (well, unless you count the stack space needed for this call and the 269 locals it creates...). This cannot fail. 270 */ 271 272 void mp_exch(mp_int *mp1, mp_int *mp2) 273 { 274 #if MP_ARGCHK == 2 275 assert(mp1 != NULL && mp2 != NULL); 276 #else 277 if(mp1 == NULL || mp2 == NULL) 278 return; 279 #endif 280 281 s_mp_exch(mp1, mp2); 282 283 } /* end mp_exch() */ 284 285 /* }}} */ 286 287 /* {{{ mp_clear(mp) */ 288 289 /* 290 mp_clear(mp) 291 292 Release the storage used by an mp_int, and void its fields so that 293 if someone calls mp_clear() again for the same int later, we won't 294 get tollchocked. 295 */ 296 297 void mp_clear(mp_int *mp) 298 { 299 if(mp == NULL) 300 return; 301 302 if(DIGITS(mp) != NULL) { 303 #if MP_CRYPTO 304 s_mp_setz(DIGITS(mp), ALLOC(mp)); 305 #endif 306 s_mp_free(DIGITS(mp), ALLOC(mp)); 307 DIGITS(mp) = NULL; 308 } 309 310 USED(mp) = 0; 311 ALLOC(mp) = 0; 312 313 } /* end mp_clear() */ 314 315 /* }}} */ 316 317 /* {{{ mp_zero(mp) */ 318 319 /* 320 mp_zero(mp) 321 322 Set mp to zero. Does not change the allocated size of the structure, 323 and therefore cannot fail (except on a bad argument, which we ignore) 324 */ 325 void mp_zero(mp_int *mp) 326 { 327 if(mp == NULL) 328 return; 329 330 s_mp_setz(DIGITS(mp), ALLOC(mp)); 331 USED(mp) = 1; 332 SIGN(mp) = ZPOS; 333 334 } /* end mp_zero() */ 335 336 /* }}} */ 337 338 /* {{{ mp_set(mp, d) */ 339 340 void mp_set(mp_int *mp, mp_digit d) 341 { 342 if(mp == NULL) 343 return; 344 345 mp_zero(mp); 346 DIGIT(mp, 0) = d; 347 348 } /* end mp_set() */ 349 350 /* }}} */ 351 352 /* {{{ mp_set_int(mp, z) */ 353 354 mp_err mp_set_int(mp_int *mp, long z) 355 { 356 int ix; 357 unsigned long v = labs(z); 358 mp_err res; 359 360 ARGCHK(mp != NULL, MP_BADARG); 361 362 mp_zero(mp); 363 if(z == 0) 364 return MP_OKAY; /* shortcut for zero */ 365 366 if (sizeof v <= sizeof(mp_digit)) { 367 DIGIT(mp,0) = v; 368 } else { 369 for (ix = sizeof(long) - 1; ix >= 0; ix--) { 370 if ((res = s_mp_mul_d(mp, (UCHAR_MAX + 1))) != MP_OKAY) 371 return res; 372 373 res = s_mp_add_d(mp, (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX)); 374 if (res != MP_OKAY) 375 return res; 376 } 377 } 378 if(z < 0) 379 SIGN(mp) = NEG; 380 381 return MP_OKAY; 382 383 } /* end mp_set_int() */ 384 385 /* }}} */ 386 387 /* {{{ mp_set_ulong(mp, z) */ 388 389 mp_err mp_set_ulong(mp_int *mp, unsigned long z) 390 { 391 int ix; 392 mp_err res; 393 394 ARGCHK(mp != NULL, MP_BADARG); 395 396 mp_zero(mp); 397 if(z == 0) 398 return MP_OKAY; /* shortcut for zero */ 399 400 if (sizeof z <= sizeof(mp_digit)) { 401 DIGIT(mp,0) = z; 402 } else { 403 for (ix = sizeof(long) - 1; ix >= 0; ix--) { 404 if ((res = s_mp_mul_d(mp, (UCHAR_MAX + 1))) != MP_OKAY) 405 return res; 406 407 res = s_mp_add_d(mp, (mp_digit)((z >> (ix * CHAR_BIT)) & UCHAR_MAX)); 408 if (res != MP_OKAY) 409 return res; 410 } 411 } 412 return MP_OKAY; 413 } /* end mp_set_ulong() */ 414 415 /* }}} */ 416 417 /*------------------------------------------------------------------------*/ 418 /* {{{ Digit arithmetic */ 419 420 /* {{{ mp_add_d(a, d, b) */ 421 422 /* 423 mp_add_d(a, d, b) 424 425 Compute the sum b = a + d, for a single digit d. Respects the sign of 426 its primary addend (single digits are unsigned anyway). 427 */ 428 429 mp_err mp_add_d(const mp_int *a, mp_digit d, mp_int *b) 430 { 431 mp_int tmp; 432 mp_err res; 433 434 ARGCHK(a != NULL && b != NULL, MP_BADARG); 435 436 if((res = mp_init_copy(&tmp, a)) != MP_OKAY) 437 return res; 438 439 if(SIGN(&tmp) == ZPOS) { 440 if((res = s_mp_add_d(&tmp, d)) != MP_OKAY) 441 goto CLEANUP; 442 } else if(s_mp_cmp_d(&tmp, d) >= 0) { 443 if((res = s_mp_sub_d(&tmp, d)) != MP_OKAY) 444 goto CLEANUP; 445 } else { 446 mp_neg(&tmp, &tmp); 447 448 DIGIT(&tmp, 0) = d - DIGIT(&tmp, 0); 449 } 450 451 if(s_mp_cmp_d(&tmp, 0) == 0) 452 SIGN(&tmp) = ZPOS; 453 454 s_mp_exch(&tmp, b); 455 456 CLEANUP: 457 mp_clear(&tmp); 458 return res; 459 460 } /* end mp_add_d() */ 461 462 /* }}} */ 463 464 /* {{{ mp_sub_d(a, d, b) */ 465 466 /* 467 mp_sub_d(a, d, b) 468 469 Compute the difference b = a - d, for a single digit d. Respects the 470 sign of its subtrahend (single digits are unsigned anyway). 471 */ 472 473 mp_err mp_sub_d(const mp_int *a, mp_digit d, mp_int *b) 474 { 475 mp_int tmp; 476 mp_err res; 477 478 ARGCHK(a != NULL && b != NULL, MP_BADARG); 479 480 if((res = mp_init_copy(&tmp, a)) != MP_OKAY) 481 return res; 482 483 if(SIGN(&tmp) == NEG) { 484 if((res = s_mp_add_d(&tmp, d)) != MP_OKAY) 485 goto CLEANUP; 486 } else if(s_mp_cmp_d(&tmp, d) >= 0) { 487 if((res = s_mp_sub_d(&tmp, d)) != MP_OKAY) 488 goto CLEANUP; 489 } else { 490 mp_neg(&tmp, &tmp); 491 492 DIGIT(&tmp, 0) = d - DIGIT(&tmp, 0); 493 SIGN(&tmp) = NEG; 494 } 495 496 if(s_mp_cmp_d(&tmp, 0) == 0) 497 SIGN(&tmp) = ZPOS; 498 499 s_mp_exch(&tmp, b); 500 501 CLEANUP: 502 mp_clear(&tmp); 503 return res; 504 505 } /* end mp_sub_d() */ 506 507 /* }}} */ 508 509 /* {{{ mp_mul_d(a, d, b) */ 510 511 /* 512 mp_mul_d(a, d, b) 513 514 Compute the product b = a * d, for a single digit d. Respects the sign 515 of its multiplicand (single digits are unsigned anyway) 516 */ 517 518 mp_err mp_mul_d(const mp_int *a, mp_digit d, mp_int *b) 519 { 520 mp_err res; 521 522 ARGCHK(a != NULL && b != NULL, MP_BADARG); 523 524 if(d == 0) { 525 mp_zero(b); 526 return MP_OKAY; 527 } 528 529 if((res = mp_copy(a, b)) != MP_OKAY) 530 return res; 531 532 res = s_mp_mul_d(b, d); 533 534 return res; 535 536 } /* end mp_mul_d() */ 537 538 /* }}} */ 539 540 /* {{{ mp_mul_2(a, c) */ 541 542 mp_err mp_mul_2(const mp_int *a, mp_int *c) 543 { 544 mp_err res; 545 546 ARGCHK(a != NULL && c != NULL, MP_BADARG); 547 548 if((res = mp_copy(a, c)) != MP_OKAY) 549 return res; 550 551 return s_mp_mul_2(c); 552 553 } /* end mp_mul_2() */ 554 555 /* }}} */ 556 557 /* {{{ mp_div_d(a, d, q, r) */ 558 559 /* 560 mp_div_d(a, d, q, r) 561 562 Compute the quotient q = a / d and remainder r = a mod d, for a 563 single digit d. Respects the sign of its divisor (single digits are 564 unsigned anyway). 565 */ 566 567 mp_err mp_div_d(const mp_int *a, mp_digit d, mp_int *q, mp_digit *r) 568 { 569 mp_err res; 570 mp_int qp; 571 mp_digit rem; 572 int pow; 573 574 ARGCHK(a != NULL, MP_BADARG); 575 576 if(d == 0) 577 return MP_RANGE; 578 579 /* Shortcut for powers of two ... */ 580 if((pow = s_mp_ispow2d(d)) >= 0) { 581 mp_digit mask; 582 583 mask = ((mp_digit)1 << pow) - 1; 584 rem = DIGIT(a, 0) & mask; 585 586 if(q) { 587 mp_copy(a, q); 588 s_mp_div_2d(q, pow); 589 } 590 591 if(r) 592 *r = rem; 593 594 return MP_OKAY; 595 } 596 597 if((res = mp_init_copy(&qp, a)) != MP_OKAY) 598 return res; 599 600 res = s_mp_div_d(&qp, d, &rem); 601 602 if(s_mp_cmp_d(&qp, 0) == 0) 603 SIGN(q) = ZPOS; 604 605 if(r) 606 *r = rem; 607 608 if(q) 609 s_mp_exch(&qp, q); 610 611 mp_clear(&qp); 612 return res; 613 614 } /* end mp_div_d() */ 615 616 /* }}} */ 617 618 /* {{{ mp_div_2(a, c) */ 619 620 /* 621 mp_div_2(a, c) 622 623 Compute c = a / 2, disregarding the remainder. 624 */ 625 626 mp_err mp_div_2(const mp_int *a, mp_int *c) 627 { 628 mp_err res; 629 630 ARGCHK(a != NULL && c != NULL, MP_BADARG); 631 632 if((res = mp_copy(a, c)) != MP_OKAY) 633 return res; 634 635 s_mp_div_2(c); 636 637 return MP_OKAY; 638 639 } /* end mp_div_2() */ 640 641 /* }}} */ 642 643 /* {{{ mp_expt_d(a, d, b) */ 644 645 mp_err mp_expt_d(const mp_int *a, mp_digit d, mp_int *c) 646 { 647 mp_int s, x; 648 mp_err res; 649 650 ARGCHK(a != NULL && c != NULL, MP_BADARG); 651 652 if((res = mp_init(&s, FLAG(a))) != MP_OKAY) 653 return res; 654 if((res = mp_init_copy(&x, a)) != MP_OKAY) 655 goto X; 656 657 DIGIT(&s, 0) = 1; 658 659 while(d != 0) { 660 if(d & 1) { 661 if((res = s_mp_mul(&s, &x)) != MP_OKAY) 662 goto CLEANUP; 663 } 664 665 d /= 2; 666 667 if((res = s_mp_sqr(&x)) != MP_OKAY) 668 goto CLEANUP; 669 } 670 671 s_mp_exch(&s, c); 672 673 CLEANUP: 674 mp_clear(&x); 675 X: 676 mp_clear(&s); 677 678 return res; 679 680 } /* end mp_expt_d() */ 681 682 /* }}} */ 683 684 /* }}} */ 685 686 /*------------------------------------------------------------------------*/ 687 /* {{{ Full arithmetic */ 688 689 /* {{{ mp_abs(a, b) */ 690 691 /* 692 mp_abs(a, b) 693 694 Compute b = |a|. 'a' and 'b' may be identical. 695 */ 696 697 mp_err mp_abs(const mp_int *a, mp_int *b) 698 { 699 mp_err res; 700 701 ARGCHK(a != NULL && b != NULL, MP_BADARG); 702 703 if((res = mp_copy(a, b)) != MP_OKAY) 704 return res; 705 706 SIGN(b) = ZPOS; 707 708 return MP_OKAY; 709 710 } /* end mp_abs() */ 711 712 /* }}} */ 713 714 /* {{{ mp_neg(a, b) */ 715 716 /* 717 mp_neg(a, b) 718 719 Compute b = -a. 'a' and 'b' may be identical. 720 */ 721 722 mp_err mp_neg(const mp_int *a, mp_int *b) 723 { 724 mp_err res; 725 726 ARGCHK(a != NULL && b != NULL, MP_BADARG); 727 728 if((res = mp_copy(a, b)) != MP_OKAY) 729 return res; 730 731 if(s_mp_cmp_d(b, 0) == MP_EQ) 732 SIGN(b) = ZPOS; 733 else 734 SIGN(b) = (SIGN(b) == NEG) ? ZPOS : NEG; 735 736 return MP_OKAY; 737 738 } /* end mp_neg() */ 739 740 /* }}} */ 741 742 /* {{{ mp_add(a, b, c) */ 743 744 /* 745 mp_add(a, b, c) 746 747 Compute c = a + b. All parameters may be identical. 748 */ 749 750 mp_err mp_add(const mp_int *a, const mp_int *b, mp_int *c) 751 { 752 mp_err res; 753 754 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); 755 756 if(SIGN(a) == SIGN(b)) { /* same sign: add values, keep sign */ 757 MP_CHECKOK( s_mp_add_3arg(a, b, c) ); 758 } else if(s_mp_cmp(a, b) >= 0) { /* different sign: |a| >= |b| */ 759 MP_CHECKOK( s_mp_sub_3arg(a, b, c) ); 760 } else { /* different sign: |a| < |b| */ 761 MP_CHECKOK( s_mp_sub_3arg(b, a, c) ); 762 } 763 764 if (s_mp_cmp_d(c, 0) == MP_EQ) 765 SIGN(c) = ZPOS; 766 767 CLEANUP: 768 return res; 769 770 } /* end mp_add() */ 771 772 /* }}} */ 773 774 /* {{{ mp_sub(a, b, c) */ 775 776 /* 777 mp_sub(a, b, c) 778 779 Compute c = a - b. All parameters may be identical. 780 */ 781 782 mp_err mp_sub(const mp_int *a, const mp_int *b, mp_int *c) 783 { 784 mp_err res; 785 int magDiff; 786 787 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); 788 789 if (a == b) { 790 mp_zero(c); 791 return MP_OKAY; 792 } 793 794 if (MP_SIGN(a) != MP_SIGN(b)) { 795 MP_CHECKOK( s_mp_add_3arg(a, b, c) ); 796 } else if (!(magDiff = s_mp_cmp(a, b))) { 797 mp_zero(c); 798 res = MP_OKAY; 799 } else if (magDiff > 0) { 800 MP_CHECKOK( s_mp_sub_3arg(a, b, c) ); 801 } else { 802 MP_CHECKOK( s_mp_sub_3arg(b, a, c) ); 803 MP_SIGN(c) = !MP_SIGN(a); 804 } 805 806 if (s_mp_cmp_d(c, 0) == MP_EQ) 807 MP_SIGN(c) = MP_ZPOS; 808 809 CLEANUP: 810 return res; 811 812 } /* end mp_sub() */ 813 814 /* }}} */ 815 816 /* {{{ mp_mul(a, b, c) */ 817 818 /* 819 mp_mul(a, b, c) 820 821 Compute c = a * b. All parameters may be identical. 822 */ 823 mp_err mp_mul(const mp_int *a, const mp_int *b, mp_int * c) 824 { 825 mp_digit *pb; 826 mp_int tmp; 827 mp_err res; 828 mp_size ib; 829 mp_size useda, usedb; 830 831 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); 832 833 if (a == c) { 834 if ((res = mp_init_copy(&tmp, a)) != MP_OKAY) 835 return res; 836 if (a == b) 837 b = &tmp; 838 a = &tmp; 839 } else if (b == c) { 840 if ((res = mp_init_copy(&tmp, b)) != MP_OKAY) 841 return res; 842 b = &tmp; 843 } else { 844 MP_DIGITS(&tmp) = 0; 845 } 846 847 if (MP_USED(a) < MP_USED(b)) { 848 const mp_int *xch = b; /* switch a and b, to do fewer outer loops */ 849 b = a; 850 a = xch; 851 } 852 853 MP_USED(c) = 1; MP_DIGIT(c, 0) = 0; 854 if((res = s_mp_pad(c, USED(a) + USED(b))) != MP_OKAY) 855 goto CLEANUP; 856 857 #ifdef NSS_USE_COMBA 858 if ((MP_USED(a) == MP_USED(b)) && IS_POWER_OF_2(MP_USED(b))) { 859 if (MP_USED(a) == 4) { 860 s_mp_mul_comba_4(a, b, c); 861 goto CLEANUP; 862 } 863 if (MP_USED(a) == 8) { 864 s_mp_mul_comba_8(a, b, c); 865 goto CLEANUP; 866 } 867 if (MP_USED(a) == 16) { 868 s_mp_mul_comba_16(a, b, c); 869 goto CLEANUP; 870 } 871 if (MP_USED(a) == 32) { 872 s_mp_mul_comba_32(a, b, c); 873 goto CLEANUP; 874 } 875 } 876 #endif 877 878 pb = MP_DIGITS(b); 879 s_mpv_mul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c)); 880 881 /* Outer loop: Digits of b */ 882 useda = MP_USED(a); 883 usedb = MP_USED(b); 884 for (ib = 1; ib < usedb; ib++) { 885 mp_digit b_i = *pb++; 886 887 /* Inner product: Digits of a */ 888 if (b_i) 889 s_mpv_mul_d_add(MP_DIGITS(a), useda, b_i, MP_DIGITS(c) + ib); 890 else 891 MP_DIGIT(c, ib + useda) = b_i; 892 } 893 894 s_mp_clamp(c); 895 896 if(SIGN(a) == SIGN(b) || s_mp_cmp_d(c, 0) == MP_EQ) 897 SIGN(c) = ZPOS; 898 else 899 SIGN(c) = NEG; 900 901 CLEANUP: 902 mp_clear(&tmp); 903 return res; 904 } /* end mp_mul() */ 905 906 /* }}} */ 907 908 /* {{{ mp_sqr(a, sqr) */ 909 910 #if MP_SQUARE 911 /* 912 Computes the square of a. This can be done more 913 efficiently than a general multiplication, because many of the 914 computation steps are redundant when squaring. The inner product 915 step is a bit more complicated, but we save a fair number of 916 iterations of the multiplication loop. 917 */ 918 919 /* sqr = a^2; Caller provides both a and tmp; */ 920 mp_err mp_sqr(const mp_int *a, mp_int *sqr) 921 { 922 mp_digit *pa; 923 mp_digit d; 924 mp_err res; 925 mp_size ix; 926 mp_int tmp; 927 int count; 928 929 ARGCHK(a != NULL && sqr != NULL, MP_BADARG); 930 931 if (a == sqr) { 932 if((res = mp_init_copy(&tmp, a)) != MP_OKAY) 933 return res; 934 a = &tmp; 935 } else { 936 DIGITS(&tmp) = 0; 937 res = MP_OKAY; 938 } 939 940 ix = 2 * MP_USED(a); 941 if (ix > MP_ALLOC(sqr)) { 942 MP_USED(sqr) = 1; 943 MP_CHECKOK( s_mp_grow(sqr, ix) ); 944 } 945 MP_USED(sqr) = ix; 946 MP_DIGIT(sqr, 0) = 0; 947 948 #ifdef NSS_USE_COMBA 949 if (IS_POWER_OF_2(MP_USED(a))) { 950 if (MP_USED(a) == 4) { 951 s_mp_sqr_comba_4(a, sqr); 952 goto CLEANUP; 953 } 954 if (MP_USED(a) == 8) { 955 s_mp_sqr_comba_8(a, sqr); 956 goto CLEANUP; 957 } 958 if (MP_USED(a) == 16) { 959 s_mp_sqr_comba_16(a, sqr); 960 goto CLEANUP; 961 } 962 if (MP_USED(a) == 32) { 963 s_mp_sqr_comba_32(a, sqr); 964 goto CLEANUP; 965 } 966 } 967 #endif 968 969 pa = MP_DIGITS(a); 970 count = MP_USED(a) - 1; 971 if (count > 0) { 972 d = *pa++; 973 s_mpv_mul_d(pa, count, d, MP_DIGITS(sqr) + 1); 974 for (ix = 3; --count > 0; ix += 2) { 975 d = *pa++; 976 s_mpv_mul_d_add(pa, count, d, MP_DIGITS(sqr) + ix); 977 } /* for(ix ...) */ 978 MP_DIGIT(sqr, MP_USED(sqr)-1) = 0; /* above loop stopped short of this. */ 979 980 /* now sqr *= 2 */ 981 s_mp_mul_2(sqr); 982 } else { 983 MP_DIGIT(sqr, 1) = 0; 984 } 985 986 /* now add the squares of the digits of a to sqr. */ 987 s_mpv_sqr_add_prop(MP_DIGITS(a), MP_USED(a), MP_DIGITS(sqr)); 988 989 SIGN(sqr) = ZPOS; 990 s_mp_clamp(sqr); 991 992 CLEANUP: 993 mp_clear(&tmp); 994 return res; 995 996 } /* end mp_sqr() */ 997 #endif 998 999 /* }}} */ 1000 1001 /* {{{ mp_div(a, b, q, r) */ 1002 1003 /* 1004 mp_div(a, b, q, r) 1005 1006 Compute q = a / b and r = a mod b. Input parameters may be re-used 1007 as output parameters. If q or r is NULL, that portion of the 1008 computation will be discarded (although it will still be computed) 1009 */ 1010 mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r) 1011 { 1012 mp_err res; 1013 mp_int *pQ, *pR; 1014 mp_int qtmp, rtmp, btmp; 1015 int cmp; 1016 mp_sign signA; 1017 mp_sign signB; 1018 1019 ARGCHK(a != NULL && b != NULL, MP_BADARG); 1020 1021 signA = MP_SIGN(a); 1022 signB = MP_SIGN(b); 1023 1024 if(mp_cmp_z(b) == MP_EQ) 1025 return MP_RANGE; 1026 1027 DIGITS(&qtmp) = 0; 1028 DIGITS(&rtmp) = 0; 1029 DIGITS(&btmp) = 0; 1030 1031 /* Set up some temporaries... */ 1032 if (!r || r == a || r == b) { 1033 MP_CHECKOK( mp_init_copy(&rtmp, a) ); 1034 pR = &rtmp; 1035 } else { 1036 MP_CHECKOK( mp_copy(a, r) ); 1037 pR = r; 1038 } 1039 1040 if (!q || q == a || q == b) { 1041 MP_CHECKOK( mp_init_size(&qtmp, MP_USED(a), FLAG(a)) ); 1042 pQ = &qtmp; 1043 } else { 1044 MP_CHECKOK( s_mp_pad(q, MP_USED(a)) ); 1045 pQ = q; 1046 mp_zero(pQ); 1047 } 1048 1049 /* 1050 If |a| <= |b|, we can compute the solution without division; 1051 otherwise, we actually do the work required. 1052 */ 1053 if ((cmp = s_mp_cmp(a, b)) <= 0) { 1054 if (cmp) { 1055 /* r was set to a above. */ 1056 mp_zero(pQ); 1057 } else { 1058 mp_set(pQ, 1); 1059 mp_zero(pR); 1060 } 1061 } else { 1062 MP_CHECKOK( mp_init_copy(&btmp, b) ); 1063 MP_CHECKOK( s_mp_div(pR, &btmp, pQ) ); 1064 } 1065 1066 /* Compute the signs for the output */ 1067 MP_SIGN(pR) = signA; /* Sr = Sa */ 1068 /* Sq = ZPOS if Sa == Sb */ /* Sq = NEG if Sa != Sb */ 1069 MP_SIGN(pQ) = (signA == signB) ? ZPOS : NEG; 1070 1071 if(s_mp_cmp_d(pQ, 0) == MP_EQ) 1072 SIGN(pQ) = ZPOS; 1073 if(s_mp_cmp_d(pR, 0) == MP_EQ) 1074 SIGN(pR) = ZPOS; 1075 1076 /* Copy output, if it is needed */ 1077 if(q && q != pQ) 1078 s_mp_exch(pQ, q); 1079 1080 if(r && r != pR) 1081 s_mp_exch(pR, r); 1082 1083 CLEANUP: 1084 mp_clear(&btmp); 1085 mp_clear(&rtmp); 1086 mp_clear(&qtmp); 1087 1088 return res; 1089 1090 } /* end mp_div() */ 1091 1092 /* }}} */ 1093 1094 /* {{{ mp_div_2d(a, d, q, r) */ 1095 1096 mp_err mp_div_2d(const mp_int *a, mp_digit d, mp_int *q, mp_int *r) 1097 { 1098 mp_err res; 1099 1100 ARGCHK(a != NULL, MP_BADARG); 1101 1102 if(q) { 1103 if((res = mp_copy(a, q)) != MP_OKAY) 1104 return res; 1105 } 1106 if(r) { 1107 if((res = mp_copy(a, r)) != MP_OKAY) 1108 return res; 1109 } 1110 if(q) { 1111 s_mp_div_2d(q, d); 1112 } 1113 if(r) { 1114 s_mp_mod_2d(r, d); 1115 } 1116 1117 return MP_OKAY; 1118 1119 } /* end mp_div_2d() */ 1120 1121 /* }}} */ 1122 1123 /* {{{ mp_expt(a, b, c) */ 1124 1125 /* 1126 mp_expt(a, b, c) 1127 1128 Compute c = a ** b, that is, raise a to the b power. Uses a 1129 standard iterative square-and-multiply technique. 1130 */ 1131 1132 mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c) 1133 { 1134 mp_int s, x; 1135 mp_err res; 1136 mp_digit d; 1137 int dig, bit; 1138 1139 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); 1140 1141 if(mp_cmp_z(b) < 0) 1142 return MP_RANGE; 1143 1144 if((res = mp_init(&s, FLAG(a))) != MP_OKAY) 1145 return res; 1146 1147 mp_set(&s, 1); 1148 1149 if((res = mp_init_copy(&x, a)) != MP_OKAY) 1150 goto X; 1151 1152 /* Loop over low-order digits in ascending order */ 1153 for(dig = 0; dig < (USED(b) - 1); dig++) { 1154 d = DIGIT(b, dig); 1155 1156 /* Loop over bits of each non-maximal digit */ 1157 for(bit = 0; bit < DIGIT_BIT; bit++) { 1158 if(d & 1) { 1159 if((res = s_mp_mul(&s, &x)) != MP_OKAY) 1160 goto CLEANUP; 1161 } 1162 1163 d >>= 1; 1164 1165 if((res = s_mp_sqr(&x)) != MP_OKAY) 1166 goto CLEANUP; 1167 } 1168 } 1169 1170 /* Consider now the last digit... */ 1171 d = DIGIT(b, dig); 1172 1173 while(d) { 1174 if(d & 1) { 1175 if((res = s_mp_mul(&s, &x)) != MP_OKAY) 1176 goto CLEANUP; 1177 } 1178 1179 d >>= 1; 1180 1181 if((res = s_mp_sqr(&x)) != MP_OKAY) 1182 goto CLEANUP; 1183 } 1184 1185 if(mp_iseven(b)) 1186 SIGN(&s) = SIGN(a); 1187 1188 res = mp_copy(&s, c); 1189 1190 CLEANUP: 1191 mp_clear(&x); 1192 X: 1193 mp_clear(&s); 1194 1195 return res; 1196 1197 } /* end mp_expt() */ 1198 1199 /* }}} */ 1200 1201 /* {{{ mp_2expt(a, k) */ 1202 1203 /* Compute a = 2^k */ 1204 1205 mp_err mp_2expt(mp_int *a, mp_digit k) 1206 { 1207 ARGCHK(a != NULL, MP_BADARG); 1208 1209 return s_mp_2expt(a, k); 1210 1211 } /* end mp_2expt() */ 1212 1213 /* }}} */ 1214 1215 /* {{{ mp_mod(a, m, c) */ 1216 1217 /* 1218 mp_mod(a, m, c) 1219 1220 Compute c = a (mod m). Result will always be 0 <= c < m. 1221 */ 1222 1223 mp_err mp_mod(const mp_int *a, const mp_int *m, mp_int *c) 1224 { 1225 mp_err res; 1226 int mag; 1227 1228 ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG); 1229 1230 if(SIGN(m) == NEG) 1231 return MP_RANGE; 1232 1233 /* 1234 If |a| > m, we need to divide to get the remainder and take the 1235 absolute value. 1236 1237 If |a| < m, we don't need to do any division, just copy and adjust 1238 the sign (if a is negative). 1239 1240 If |a| == m, we can simply set the result to zero. 1241 1242 This order is intended to minimize the average path length of the 1243 comparison chain on common workloads -- the most frequent cases are 1244 that |a| != m, so we do those first. 1245 */ 1246 if((mag = s_mp_cmp(a, m)) > 0) { 1247 if((res = mp_div(a, m, NULL, c)) != MP_OKAY) 1248 return res; 1249 1250 if(SIGN(c) == NEG) { 1251 if((res = mp_add(c, m, c)) != MP_OKAY) 1252 return res; 1253 } 1254 1255 } else if(mag < 0) { 1256 if((res = mp_copy(a, c)) != MP_OKAY) 1257 return res; 1258 1259 if(mp_cmp_z(a) < 0) { 1260 if((res = mp_add(c, m, c)) != MP_OKAY) 1261 return res; 1262 1263 } 1264 1265 } else { 1266 mp_zero(c); 1267 1268 } 1269 1270 return MP_OKAY; 1271 1272 } /* end mp_mod() */ 1273 1274 /* }}} */ 1275 1276 /* {{{ mp_mod_d(a, d, c) */ 1277 1278 /* 1279 mp_mod_d(a, d, c) 1280 1281 Compute c = a (mod d). Result will always be 0 <= c < d 1282 */ 1283 mp_err mp_mod_d(const mp_int *a, mp_digit d, mp_digit *c) 1284 { 1285 mp_err res; 1286 mp_digit rem; 1287 1288 ARGCHK(a != NULL && c != NULL, MP_BADARG); 1289 1290 if(s_mp_cmp_d(a, d) > 0) { 1291 if((res = mp_div_d(a, d, NULL, &rem)) != MP_OKAY) 1292 return res; 1293 1294 } else { 1295 if(SIGN(a) == NEG) 1296 rem = d - DIGIT(a, 0); 1297 else 1298 rem = DIGIT(a, 0); 1299 } 1300 1301 if(c) 1302 *c = rem; 1303 1304 return MP_OKAY; 1305 1306 } /* end mp_mod_d() */ 1307 1308 /* }}} */ 1309 1310 /* {{{ mp_sqrt(a, b) */ 1311 1312 /* 1313 mp_sqrt(a, b) 1314 1315 Compute the integer square root of a, and store the result in b. 1316 Uses an integer-arithmetic version of Newton's iterative linear 1317 approximation technique to determine this value; the result has the 1318 following two properties: 1319 1320 b^2 <= a 1321 (b+1)^2 >= a 1322 1323 It is a range error to pass a negative value. 1324 */ 1325 mp_err mp_sqrt(const mp_int *a, mp_int *b) 1326 { 1327 mp_int x, t;