git.sur5r.net Git - freertos/blob - FreeRTOS-Plus/Source/CyaSSL/ctaocrypt/src/integer.c

   1 /* integer.c
   2  *
   3  * Copyright (C) 2006-2014 wolfSSL Inc.
   4  *
   5  * This file is part of CyaSSL.
   6  *
   7  * CyaSSL is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 2 of the License, or
  10  * (at your option) any later version.
  11  *
  12  * CyaSSL is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  20  */
  21
  22
  23 /*
  24  * Based on public domain LibTomMath 0.38 by Tom St Denis, tomstdenis@iahu.ca,
  25  * http://math.libtomcrypt.com
  26  */
  27
  28
  29 #ifdef HAVE_CONFIG_H
  30     #include <config.h>
  31 #endif
  32
  33 /* in case user set USE_FAST_MATH there */
  34 #include <cyassl/ctaocrypt/settings.h>
  35
  36 #ifndef NO_BIG_INT
  37
  38 #ifndef USE_FAST_MATH
  39
  40 #include <cyassl/ctaocrypt/integer.h>
  41
  42 #ifndef NO_CYASSL_SMALL_STACK
  43     #ifndef CYASSL_SMALL_STACK
  44         #define CYASSL_SMALL_STACK
  45     #endif
  46 #endif
  47
  48 static void bn_reverse (unsigned char *s, int len);
  49
  50 /* math settings check */
  51 word32 CheckRunTimeSettings(void)
  52 {
  53     return CTC_SETTINGS;
  54 }
  55
  56
  57 /* handle up to 6 inits */
  58 int mp_init_multi(mp_int* a, mp_int* b, mp_int* c, mp_int* d, mp_int* e,
  59                   mp_int* f)
  60 {
  61     int res = MP_OKAY;
  62
  63     if (a && ((res = mp_init(a)) != MP_OKAY))
  64         return res;
  65
  66     if (b && ((res = mp_init(b)) != MP_OKAY)) {
  67         mp_clear(a);
  68         return res;
  69     }
  70
  71     if (c && ((res = mp_init(c)) != MP_OKAY)) {
  72         mp_clear(a); mp_clear(b);
  73         return res;
  74     }
  75
  76     if (d && ((res = mp_init(d)) != MP_OKAY)) {
  77         mp_clear(a); mp_clear(b); mp_clear(c);
  78         return res;
  79     }
  80
  81     if (e && ((res = mp_init(e)) != MP_OKAY)) {
  82         mp_clear(a); mp_clear(b); mp_clear(c); mp_clear(d);
  83         return res;
  84     }
  85
  86     if (f && ((res = mp_init(f)) != MP_OKAY)) {
  87         mp_clear(a); mp_clear(b); mp_clear(c); mp_clear(d); mp_clear(e);
  88         return res;
  89     }
  90
  91     return res;
  92 }
  93
  94
  95 /* init a new mp_int */
  96 int mp_init (mp_int * a)
  97 {
  98   int i;
  99
 100   /* allocate memory required and clear it */
 101   a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * MP_PREC, 0,
 102                                       DYNAMIC_TYPE_BIGINT);
 103   if (a->dp == NULL) {
 104     return MP_MEM;
 105   }
 106
 107   /* set the digits to zero */
 108   for (i = 0; i < MP_PREC; i++) {
 109       a->dp[i] = 0;
 110   }
 111
 112   /* set the used to zero, allocated digits to the default precision
 113    * and sign to positive */
 114   a->used  = 0;
 115   a->alloc = MP_PREC;
 116   a->sign  = MP_ZPOS;
 117
 118   return MP_OKAY;
 119 }
 120
 121
 122 /* clear one (frees)  */
 123 void
 124 mp_clear (mp_int * a)
 125 {
 126   int i;
 127
 128   if (a == NULL)
 129       return;
 130
 131   /* only do anything if a hasn't been freed previously */
 132   if (a->dp != NULL) {
 133     /* first zero the digits */
 134     for (i = 0; i < a->used; i++) {
 135         a->dp[i] = 0;
 136     }
 137
 138     /* free ram */
 139     XFREE(a->dp, 0, DYNAMIC_TYPE_BIGINT);
 140
 141     /* reset members to make debugging easier */
 142     a->dp    = NULL;
 143     a->alloc = a->used = 0;
 144     a->sign  = MP_ZPOS;
 145   }
 146 }
 147
 148
 149 /* get the size for an unsigned equivalent */
 150 int mp_unsigned_bin_size (mp_int * a)
 151 {
 152   int     size = mp_count_bits (a);
 153   return (size / 8 + ((size & 7) != 0 ? 1 : 0));
 154 }
 155
 156
 157 /* returns the number of bits in an int */
 158 int
 159 mp_count_bits (mp_int * a)
 160 {
 161   int     r;
 162   mp_digit q;
 163
 164   /* shortcut */
 165   if (a->used == 0) {
 166     return 0;
 167   }
 168
 169   /* get number of digits and add that */
 170   r = (a->used - 1) * DIGIT_BIT;
 171
 172   /* take the last digit and count the bits in it */
 173   q = a->dp[a->used - 1];
 174   while (q > ((mp_digit) 0)) {
 175     ++r;
 176     q >>= ((mp_digit) 1);
 177   }
 178   return r;
 179 }
 180
 181
 182 int mp_leading_bit (mp_int * a)
 183 {
 184     int bit = 0;
 185     mp_int t;
 186
 187     if (mp_init_copy(&t, a) != MP_OKAY)
 188         return 0;
 189
 190     while (mp_iszero(&t) == 0) {
 191 #ifndef MP_8BIT
 192         bit = (t.dp[0] & 0x80) != 0;
 193 #else
 194         bit = (t.dp[0] | ((t.dp[1] & 0x01) << 7)) & 0x80 != 0;
 195 #endif
 196         if (mp_div_2d (&t, 8, &t, NULL) != MP_OKAY)
 197             break;
 198     }
 199     mp_clear(&t);
 200     return bit;
 201 }
 202
 203
 204 /* store in unsigned [big endian] format */
 205 int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
 206 {
 207   int     x, res;
 208   mp_int  t;
 209
 210   if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
 211     return res;
 212   }
 213
 214   x = 0;
 215   while (mp_iszero (&t) == 0) {
 216 #ifndef MP_8BIT
 217       b[x++] = (unsigned char) (t.dp[0] & 255);
 218 #else
 219       b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7));
 220 #endif
 221     if ((res = mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) {
 222       mp_clear (&t);
 223       return res;
 224     }
 225   }
 226   bn_reverse (b, x);
 227   mp_clear (&t);
 228   return MP_OKAY;
 229 }
 230
 231
 232 /* creates "a" then copies b into it */
 233 int mp_init_copy (mp_int * a, mp_int * b)
 234 {
 235   int     res;
 236
 237   if ((res = mp_init (a)) != MP_OKAY) {
 238     return res;
 239   }
 240   return mp_copy (b, a);
 241 }
 242
 243
 244 /* copy, b = a */
 245 int
 246 mp_copy (mp_int * a, mp_int * b)
 247 {
 248   int     res, n;
 249
 250   /* if dst == src do nothing */
 251   if (a == b) {
 252     return MP_OKAY;
 253   }
 254
 255   /* grow dest */
 256   if (b->alloc < a->used) {
 257      if ((res = mp_grow (b, a->used)) != MP_OKAY) {
 258         return res;
 259      }
 260   }
 261
 262   /* zero b and copy the parameters over */
 263   {
 264     register mp_digit *tmpa, *tmpb;
 265
 266     /* pointer aliases */
 267
 268     /* source */
 269     tmpa = a->dp;
 270
 271     /* destination */
 272     tmpb = b->dp;
 273
 274     /* copy all the digits */
 275     for (n = 0; n < a->used; n++) {
 276       *tmpb++ = *tmpa++;
 277     }
 278
 279     /* clear high digits */
 280     for (; n < b->used; n++) {
 281       *tmpb++ = 0;
 282     }
 283   }
 284
 285   /* copy used count and sign */
 286   b->used = a->used;
 287   b->sign = a->sign;
 288   return MP_OKAY;
 289 }
 290
 291
 292 /* grow as required */
 293 int mp_grow (mp_int * a, int size)
 294 {
 295   int     i;
 296   mp_digit *tmp;
 297
 298   /* if the alloc size is smaller alloc more ram */
 299   if (a->alloc < size) {
 300     /* ensure there are always at least MP_PREC digits extra on top */
 301     size += (MP_PREC * 2) - (size % MP_PREC);
 302
 303     /* reallocate the array a->dp
 304      *
 305      * We store the return in a temporary variable
 306      * in case the operation failed we don't want
 307      * to overwrite the dp member of a.
 308      */
 309     tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * size, 0,
 310                                        DYNAMIC_TYPE_BIGINT);
 311     if (tmp == NULL) {
 312       /* reallocation failed but "a" is still valid [can be freed] */
 313       return MP_MEM;
 314     }
 315
 316     /* reallocation succeeded so set a->dp */
 317     a->dp = tmp;
 318
 319     /* zero excess digits */
 320     i        = a->alloc;
 321     a->alloc = size;
 322     for (; i < a->alloc; i++) {
 323       a->dp[i] = 0;
 324     }
 325   }
 326   return MP_OKAY;
 327 }
 328
 329
 330 /* reverse an array, used for radix code */
 331 void
 332 bn_reverse (unsigned char *s, int len)
 333 {
 334   int     ix, iy;
 335   unsigned char t;
 336
 337   ix = 0;
 338   iy = len - 1;
 339   while (ix < iy) {
 340     t     = s[ix];
 341     s[ix] = s[iy];
 342     s[iy] = t;
 343     ++ix;
 344     --iy;
 345   }
 346 }
 347
 348
 349 /* shift right by a certain bit count (store quotient in c, optional
 350    remainder in d) */
 351 int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
 352 {
 353   int     D, res;
 354   mp_int  t;
 355
 356
 357   /* if the shift count is <= 0 then we do no work */
 358   if (b <= 0) {
 359     res = mp_copy (a, c);
 360     if (d != NULL) {
 361       mp_zero (d);
 362     }
 363     return res;
 364   }
 365
 366   if ((res = mp_init (&t)) != MP_OKAY) {
 367     return res;
 368   }
 369
 370   /* get the remainder */
 371   if (d != NULL) {
 372     if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) {
 373       mp_clear (&t);
 374       return res;
 375     }
 376   }
 377
 378   /* copy */
 379   if ((res = mp_copy (a, c)) != MP_OKAY) {
 380     mp_clear (&t);
 381     return res;
 382   }
 383
 384   /* shift by as many digits in the bit count */
 385   if (b >= (int)DIGIT_BIT) {
 386     mp_rshd (c, b / DIGIT_BIT);
 387   }
 388
 389   /* shift any bit count < DIGIT_BIT */
 390   D = (b % DIGIT_BIT);
 391   if (D != 0) {
 392     mp_rshb(c, D);
 393   }
 394   mp_clamp (c);
 395   if (d != NULL) {
 396     mp_exch (&t, d);
 397   }
 398   mp_clear (&t);
 399   return MP_OKAY;
 400 }
 401
 402
 403 /* set to zero */
 404 void mp_zero (mp_int * a)
 405 {
 406   int       n;
 407   mp_digit *tmp;
 408
 409   a->sign = MP_ZPOS;
 410   a->used = 0;
 411
 412   tmp = a->dp;
 413   for (n = 0; n < a->alloc; n++) {
 414      *tmp++ = 0;
 415   }
 416 }
 417
 418
 419 /* trim unused digits
 420  *
 421  * This is used to ensure that leading zero digits are
 422  * trimed and the leading "used" digit will be non-zero
 423  * Typically very fast.  Also fixes the sign if there
 424  * are no more leading digits
 425  */
 426 void
 427 mp_clamp (mp_int * a)
 428 {
 429   /* decrease used while the most significant digit is
 430    * zero.
 431    */
 432   while (a->used > 0 && a->dp[a->used - 1] == 0) {
 433     --(a->used);
 434   }
 435
 436   /* reset the sign flag if used == 0 */
 437   if (a->used == 0) {
 438     a->sign = MP_ZPOS;
 439   }
 440 }
 441
 442
 443 /* swap the elements of two integers, for cases where you can't simply swap the
 444  * mp_int pointers around
 445  */
 446 void
 447 mp_exch (mp_int * a, mp_int * b)
 448 {
 449   mp_int  t;
 450
 451   t  = *a;
 452   *a = *b;
 453   *b = t;
 454 }
 455
 456
 457 /* shift right a certain number of bits */
 458 void mp_rshb (mp_int *c, int x)
 459 {
 460     register mp_digit *tmpc, mask, shift;
 461     mp_digit r, rr;
 462     mp_digit D = x;
 463
 464     /* mask */
 465     mask = (((mp_digit)1) << D) - 1;
 466
 467     /* shift for lsb */
 468     shift = DIGIT_BIT - D;
 469
 470     /* alias */
 471     tmpc = c->dp + (c->used - 1);
 472
 473     /* carry */
 474     r = 0;
 475     for (x = c->used - 1; x >= 0; x--) {
 476       /* get the lower  bits of this word in a temp */
 477       rr = *tmpc & mask;
 478
 479       /* shift the current word and mix in the carry bits from previous word */
 480       *tmpc = (*tmpc >> D) | (r << shift);
 481       --tmpc;
 482
 483       /* set the carry to the carry bits of the current word found above */
 484       r = rr;
 485     }
 486 }
 487
 488
 489 /* shift right a certain amount of digits */
 490 void mp_rshd (mp_int * a, int b)
 491 {
 492   int     x;
 493
 494   /* if b <= 0 then ignore it */
 495   if (b <= 0) {
 496     return;
 497   }
 498
 499   /* if b > used then simply zero it and return */
 500   if (a->used <= b) {
 501     mp_zero (a);
 502     return;
 503   }
 504
 505   {
 506     register mp_digit *bottom, *top;
 507
 508     /* shift the digits down */
 509
 510     /* bottom */
 511     bottom = a->dp;
 512
 513     /* top [offset into digits] */
 514     top = a->dp + b;
 515
 516     /* this is implemented as a sliding window where
 517      * the window is b-digits long and digits from
 518      * the top of the window are copied to the bottom
 519      *
 520      * e.g.
 521
 522      b-2 | b-1 | b0 | b1 | b2 | ... | bb |   ---->
 523                  /\                   |      ---->
 524                   \-------------------/      ---->
 525      */
 526     for (x = 0; x < (a->used - b); x++) {
 527       *bottom++ = *top++;
 528     }
 529
 530     /* zero the top digits */
 531     for (; x < a->used; x++) {
 532       *bottom++ = 0;
 533     }
 534   }
 535
 536   /* remove excess digits */
 537   a->used -= b;
 538 }
 539
 540
 541 /* calc a value mod 2**b */
 542 int
 543 mp_mod_2d (mp_int * a, int b, mp_int * c)
 544 {
 545   int     x, res;
 546
 547   /* if b is <= 0 then zero the int */
 548   if (b <= 0) {
 549     mp_zero (c);
 550     return MP_OKAY;
 551   }
 552
 553   /* if the modulus is larger than the value than return */
 554   if (b >= (int) (a->used * DIGIT_BIT)) {
 555     res = mp_copy (a, c);
 556     return res;
 557   }
 558
 559   /* copy */
 560   if ((res = mp_copy (a, c)) != MP_OKAY) {
 561     return res;
 562   }
 563
 564   /* zero digits above the last digit of the modulus */
 565   for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
 566     c->dp[x] = 0;
 567   }
 568   /* clear the digit that is not completely outside/inside the modulus */
 569   c->dp[b / DIGIT_BIT] &= (mp_digit) ((((mp_digit) 1) <<
 570               (((mp_digit) b) % DIGIT_BIT)) - ((mp_digit) 1));
 571   mp_clamp (c);
 572   return MP_OKAY;
 573 }
 574
 575
 576 /* reads a unsigned char array, assumes the msb is stored first [big endian] */
 577 int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c)
 578 {
 579   int     res;
 580
 581   /* make sure there are at least two digits */
 582   if (a->alloc < 2) {
 583      if ((res = mp_grow(a, 2)) != MP_OKAY) {
 584         return res;
 585      }
 586   }
 587
 588   /* zero the int */
 589   mp_zero (a);
 590
 591   /* read the bytes in */
 592   while (c-- > 0) {
 593     if ((res = mp_mul_2d (a, 8, a)) != MP_OKAY) {
 594       return res;
 595     }
 596
 597 #ifndef MP_8BIT
 598       a->dp[0] |= *b++;
 599       a->used += 1;
 600 #else
 601       a->dp[0] = (*b & MP_MASK);
 602       a->dp[1] |= ((*b++ >> 7U) & 1);
 603       a->used += 2;
 604 #endif
 605   }
 606   mp_clamp (a);
 607   return MP_OKAY;
 608 }
 609
 610
 611 /* shift left by a certain bit count */
 612 int mp_mul_2d (mp_int * a, int b, mp_int * c)
 613 {
 614   mp_digit d;
 615   int      res;
 616
 617   /* copy */
 618   if (a != c) {
 619      if ((res = mp_copy (a, c)) != MP_OKAY) {
 620        return res;
 621      }
 622   }
 623
 624   if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) {
 625      if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
 626        return res;
 627      }
 628   }
 629
 630   /* shift by as many digits in the bit count */
 631   if (b >= (int)DIGIT_BIT) {
 632     if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
 633       return res;
 634     }
 635   }
 636
 637   /* shift any bit count < DIGIT_BIT */
 638   d = (mp_digit) (b % DIGIT_BIT);
 639   if (d != 0) {
 640     register mp_digit *tmpc, shift, mask, r, rr;
 641     register int x;
 642
 643     /* bitmask for carries */
 644     mask = (((mp_digit)1) << d) - 1;
 645
 646     /* shift for msbs */
 647     shift = DIGIT_BIT - d;
 648
 649     /* alias */
 650     tmpc = c->dp;
 651
 652     /* carry */
 653     r    = 0;
 654     for (x = 0; x < c->used; x++) {
 655       /* get the higher bits of the current word */
 656       rr = (*tmpc >> shift) & mask;
 657
 658       /* shift the current word and OR in the carry */
 659       *tmpc = ((*tmpc << d) | r) & MP_MASK;
 660       ++tmpc;
 661
 662       /* set the carry to the carry bits of the current word */
 663       r = rr;
 664     }
 665
 666     /* set final carry */
 667     if (r != 0) {
 668        c->dp[(c->used)++] = r;
 669     }
 670   }
 671   mp_clamp (c);
 672   return MP_OKAY;
 673 }
 674
 675
 676 /* shift left a certain amount of digits */
 677 int mp_lshd (mp_int * a, int b)
 678 {
 679   int     x, res;
 680
 681   /* if its less than zero return */
 682   if (b <= 0) {
 683     return MP_OKAY;
 684   }
 685
 686   /* grow to fit the new digits */
 687   if (a->alloc < a->used + b) {
 688      if ((res = mp_grow (a, a->used + b)) != MP_OKAY) {
 689        return res;
 690      }
 691   }
 692
 693   {
 694     register mp_digit *top, *bottom;
 695
 696     /* increment the used by the shift amount then copy upwards */
 697     a->used += b;
 698
 699     /* top */
 700     top = a->dp + a->used - 1;
 701
 702     /* base */
 703     bottom = a->dp + a->used - 1 - b;
 704
 705     /* much like mp_rshd this is implemented using a sliding window
 706      * except the window goes the otherway around.  Copying from
 707      * the bottom to the top.  see bn_mp_rshd.c for more info.
 708      */
 709     for (x = a->used - 1; x >= b; x--) {
 710       *top-- = *bottom--;
 711     }
 712
 713     /* zero the lower digits */
 714     top = a->dp;
 715     for (x = 0; x < b; x++) {
 716       *top++ = 0;
 717     }
 718   }
 719   return MP_OKAY;
 720 }
 721
 722
 723 /* this is a shell function that calls either the normal or Montgomery
 724  * exptmod functions.  Originally the call to the montgomery code was
 725  * embedded in the normal function but that wasted alot of stack space
 726  * for nothing (since 99% of the time the Montgomery code would be called)
 727  */
 728 int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 729 {
 730   int dr;
 731
 732   /* modulus P must be positive */
 733   if (P->sign == MP_NEG) {
 734      return MP_VAL;
 735   }
 736
 737   /* if exponent X is negative we have to recurse */
 738   if (X->sign == MP_NEG) {
 739 #ifdef BN_MP_INVMOD_C
 740      mp_int tmpG, tmpX;
 741      int err;
 742
 743      /* first compute 1/G mod P */
 744      if ((err = mp_init(&tmpG)) != MP_OKAY) {
 745         return err;
 746      }
 747      if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {
 748         mp_clear(&tmpG);
 749         return err;
 750      }
 751
 752      /* now get |X| */
 753      if ((err = mp_init(&tmpX)) != MP_OKAY) {
 754         mp_clear(&tmpG);
 755         return err;
 756      }
 757      if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
 758         mp_clear(&tmpG);
 759         mp_clear(&tmpX);
 760         return err;
 761      }
 762
 763      /* and now compute (1/G)**|X| instead of G**X [X < 0] */
 764      err = mp_exptmod(&tmpG, &tmpX, P, Y);
 765      mp_clear(&tmpG);
 766      mp_clear(&tmpX);
 767      return err;
 768 #else
 769      /* no invmod */
 770      return MP_VAL;
 771 #endif
 772   }
 773
 774 /* modified diminished radix reduction */
 775 #if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C) && \
 776   defined(BN_S_MP_EXPTMOD_C)
 777   if (mp_reduce_is_2k_l(P) == MP_YES) {
 778      return s_mp_exptmod(G, X, P, Y, 1);
 779   }
 780 #endif
 781
 782 #ifdef BN_MP_DR_IS_MODULUS_C
 783   /* is it a DR modulus? */
 784   dr = mp_dr_is_modulus(P);
 785 #else
 786   /* default to no */
 787   dr = 0;
 788 #endif
 789
 790 #ifdef BN_MP_REDUCE_IS_2K_C
 791   /* if not, is it a unrestricted DR modulus? */
 792   if (dr == 0) {
 793      dr = mp_reduce_is_2k(P) << 1;
 794   }
 795 #endif
 796
 797   /* if the modulus is odd or dr != 0 use the montgomery method */
 798 #ifdef BN_MP_EXPTMOD_FAST_C
 799   if (mp_isodd (P) == 1 || dr !=  0) {
 800     return mp_exptmod_fast (G, X, P, Y, dr);
 801   } else {
 802 #endif
 803 #ifdef BN_S_MP_EXPTMOD_C
 804     /* otherwise use the generic Barrett reduction technique */
 805     return s_mp_exptmod (G, X, P, Y, 0);
 806 #else
 807     /* no exptmod for evens */
 808     return MP_VAL;
 809 #endif
 810 #ifdef BN_MP_EXPTMOD_FAST_C
 811   }
 812 #endif
 813 }
 814
 815
 816 /* b = |a|
 817  *
 818  * Simple function copies the input and fixes the sign to positive
 819  */
 820 int
 821 mp_abs (mp_int * a, mp_int * b)
 822 {
 823   int     res;
 824
 825   /* copy a to b */
 826   if (a != b) {
 827      if ((res = mp_copy (a, b)) != MP_OKAY) {
 828        return res;
 829      }
 830   }
 831
 832   /* force the sign of b to positive */
 833   b->sign = MP_ZPOS;
 834
 835   return MP_OKAY;
 836 }
 837
 838
 839 /* hac 14.61, pp608 */
 840 int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
 841 {
 842   /* b cannot be negative */
 843   if (b->sign == MP_NEG || mp_iszero(b) == 1) {
 844     return MP_VAL;
 845   }
 846
 847 #ifdef BN_FAST_MP_INVMOD_C
 848   /* if the modulus is odd we can use a faster routine instead */
 849   if (mp_isodd (b) == 1) {
 850     return fast_mp_invmod (a, b, c);
 851   }
 852 #endif
 853
 854 #ifdef BN_MP_INVMOD_SLOW_C
 855   return mp_invmod_slow(a, b, c);
 856 #endif
 857 }
 858
 859
 860 /* computes the modular inverse via binary extended euclidean algorithm,
 861  * that is c = 1/a mod b
 862  *
 863  * Based on slow invmod except this is optimized for the case where b is
 864  * odd as per HAC Note 14.64 on pp. 610
 865  */
 866 int fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
 867 {
 868   mp_int  x, y, u, v, B, D;
 869   int     res, neg;
 870
 871   /* 2. [modified] b must be odd   */
 872   if (mp_iseven (b) == 1) {
 873     return MP_VAL;
 874   }
 875
 876   /* init all our temps */
 877   if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D)) != MP_OKAY) {
 878      return res;
 879   }
 880
 881   /* x == modulus, y == value to invert */
 882   if ((res = mp_copy (b, &x)) != MP_OKAY) {
 883     goto LBL_ERR;
 884   }
 885
 886   /* we need y = |a| */
 887   if ((res = mp_mod (a, b, &y)) != MP_OKAY) {
 888     goto LBL_ERR;
 889   }
 890
 891   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
 892   if ((res = mp_copy (&x, &u)) != MP_OKAY) {
 893     goto LBL_ERR;
 894   }
 895   if ((res = mp_copy (&y, &v)) != MP_OKAY) {
 896     goto LBL_ERR;
 897   }
 898   mp_set (&D, 1);
 899
 900 top:
 901   /* 4.  while u is even do */
 902   while (mp_iseven (&u) == 1) {
 903     /* 4.1 u = u/2 */
 904     if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
 905       goto LBL_ERR;
 906     }
 907     /* 4.2 if B is odd then */
 908     if (mp_isodd (&B) == 1) {
 909       if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
 910         goto LBL_ERR;
 911       }
 912     }
 913     /* B = B/2 */
 914     if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
 915       goto LBL_ERR;
 916     }
 917   }
 918
 919   /* 5.  while v is even do */
 920   while (mp_iseven (&v) == 1) {
 921     /* 5.1 v = v/2 */
 922     if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
 923       goto LBL_ERR;
 924     }
 925     /* 5.2 if D is odd then */
 926     if (mp_isodd (&D) == 1) {
 927       /* D = (D-x)/2 */
 928       if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
 929         goto LBL_ERR;
 930       }
 931     }
 932     /* D = D/2 */
 933     if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
 934       goto LBL_ERR;
 935     }
 936   }
 937
 938   /* 6.  if u >= v then */
 939   if (mp_cmp (&u, &v) != MP_LT) {
 940     /* u = u - v, B = B - D */
 941     if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
 942       goto LBL_ERR;
 943     }
 944
 945     if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
 946       goto LBL_ERR;
 947     }
 948   } else {
 949     /* v - v - u, D = D - B */
 950     if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
 951       goto LBL_ERR;
 952     }
 953
 954     if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
 955       goto LBL_ERR;
 956     }
 957   }
 958
 959   /* if not zero goto step 4 */
 960   if (mp_iszero (&u) == 0) {
 961     goto top;
 962   }
 963
 964   /* now a = C, b = D, gcd == g*v */
 965
 966   /* if v != 1 then there is no inverse */
 967   if (mp_cmp_d (&v, 1) != MP_EQ) {
 968     res = MP_VAL;
 969     goto LBL_ERR;
 970   }
 971
 972   /* b is now the inverse */
 973   neg = a->sign;
 974   while (D.sign == MP_NEG) {
 975     if ((res = mp_add (&D, b, &D)) != MP_OKAY) {
 976       goto LBL_ERR;
 977     }
 978   }
 979   mp_exch (&D, c);
 980   c->sign = neg;
 981   res = MP_OKAY;
 982
 983 LBL_ERR:mp_clear(&x);
 984         mp_clear(&y);
 985         mp_clear(&u);
 986         mp_clear(&v);
 987         mp_clear(&B);
 988         mp_clear(&D);
 989   return res;
 990 }
 991
 992
 993 /* hac 14.61, pp608 */
 994 int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c)
 995 {
 996   mp_int  x, y, u, v, A, B, C, D;
 997   int     res;
 998
 999   /* b cannot be negative */
1000   if (b->sign == MP_NEG || mp_iszero(b) == 1) {
1001     return MP_VAL;
1002   }
1003
1004   /* init temps */
1005   if ((res = mp_init_multi(&x, &y, &u, &v,
1006                            &A, &B)) != MP_OKAY) {
1007      return res;
1008   }
1009
1010   /* init rest of tmps temps */
1011   if ((res = mp_init_multi(&C, &D, 0, 0, 0, 0)) != MP_OKAY) {
1012      return res;
1013   }
1014
1015   /* x = a, y = b */
1016   if ((res = mp_mod(a, b, &x)) != MP_OKAY) {
1017       goto LBL_ERR;
1018   }
1019   if ((res = mp_copy (b, &y)) != MP_OKAY) {
1020     goto LBL_ERR;
1021   }
1022
1023   /* 2. [modified] if x,y are both even then return an error! */
1024   if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) {
1025     res = MP_VAL;
1026     goto LBL_ERR;
1027   }
1028
1029   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
1030   if ((res = mp_copy (&x, &u)) != MP_OKAY) {
1031     goto LBL_ERR;
1032   }
1033   if ((res = mp_copy (&y, &v)) != MP_OKAY) {
1034     goto LBL_ERR;
1035   }
1036   mp_set (&A, 1);
1037   mp_set (&D, 1);
1038
1039 top:
1040   /* 4.  while u is even do */
1041   while (mp_iseven (&u) == 1) {
1042     /* 4.1 u = u/2 */
1043     if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
1044       goto LBL_ERR;
1045     }
1046     /* 4.2 if A or B is odd then */
1047     if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) {
1048       /* A = (A+y)/2, B = (B-x)/2 */
1049       if ((res = mp_add (&A, &y, &A)) != MP_OKAY) {
1050          goto LBL_ERR;
1051       }
1052       if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
1053          goto LBL_ERR;
1054       }
1055     }
1056     /* A = A/2, B = B/2 */
1057     if ((res = mp_div_2 (&A, &A)) != MP_OKAY) {
1058       goto LBL_ERR;
1059     }
1060     if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
1061       goto LBL_ERR;
1062     }
1063   }
1064
1065   /* 5.  while v is even do */
1066   while (mp_iseven (&v) == 1) {
1067     /* 5.1 v = v/2 */
1068     if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
1069       goto LBL_ERR;
1070     }
1071     /* 5.2 if C or D is odd then */
1072     if (mp_isodd (&C) == 1 || mp_isodd (&D) == 1) {
1073       /* C = (C+y)/2, D = (D-x)/2 */
1074       if ((res = mp_add (&C, &y, &C)) != MP_OKAY) {
1075          goto LBL_ERR;
1076       }
1077       if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
1078          goto LBL_ERR;
1079       }
1080     }
1081     /* C = C/2, D = D/2 */
1082     if ((res = mp_div_2 (&C, &C)) != MP_OKAY) {
1083       goto LBL_ERR;
1084     }
1085     if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
1086       goto LBL_ERR;
1087     }
1088   }
1089
1090   /* 6.  if u >= v then */
1091   if (mp_cmp (&u, &v) != MP_LT) {
1092     /* u = u - v, A = A - C, B = B - D */
1093     if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
1094       goto LBL_ERR;
1095     }
1096
1097     if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) {
1098       goto LBL_ERR;
1099     }
1100
1101     if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
1102       goto LBL_ERR;
1103     }
1104   } else {
1105     /* v - v - u, C = C - A, D = D - B */
1106     if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
1107       goto LBL_ERR;
1108     }
1109
1110     if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) {
1111       goto LBL_ERR;
1112     }
1113
1114     if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
1115       goto LBL_ERR;
1116     }
1117   }
1118
1119   /* if not zero goto step 4 */
1120   if (mp_iszero (&u) == 0)
1121     goto top;
1122
1123   /* now a = C, b = D, gcd == g*v */
1124
1125   /* if v != 1 then there is no inverse */
1126   if (mp_cmp_d (&v, 1) != MP_EQ) {
1127     res = MP_VAL;
1128     goto LBL_ERR;
1129   }
1130
1131   /* if its too low */
1132   while (mp_cmp_d(&C, 0) == MP_LT) {
1133       if ((res = mp_add(&C, b, &C)) != MP_OKAY) {
1134          goto LBL_ERR;
1135       }
1136   }
1137
1138   /* too big */
1139   while (mp_cmp_mag(&C, b) != MP_LT) {
1140       if ((res = mp_sub(&C, b, &C)) != MP_OKAY) {
1141          goto LBL_ERR;
1142       }
1143   }
1144
1145   /* C is now the inverse */
1146   mp_exch (&C, c);
1147   res = MP_OKAY;
1148 LBL_ERR:mp_clear(&x);
1149         mp_clear(&y);
1150         mp_clear(&u);
1151         mp_clear(&v);
1152         mp_clear(&A);
1153         mp_clear(&B);
1154         mp_clear(&C);
1155         mp_clear(&D);
1156   return res;
1157 }
1158
1159
1160 /* compare maginitude of two ints (unsigned) */
1161 int mp_cmp_mag (mp_int * a, mp_int * b)
1162 {
1163   int     n;
1164   mp_digit *tmpa, *tmpb;
1165
1166   /* compare based on # of non-zero digits */
1167   if (a->used > b->used) {
1168     return MP_GT;
1169   }
1170
1171   if (a->used < b->used) {
1172     return MP_LT;
1173   }
1174
1175   /* alias for a */
1176   tmpa = a->dp + (a->used - 1);
1177
1178   /* alias for b */
1179   tmpb = b->dp + (a->used - 1);
1180
1181   /* compare based on digits  */
1182   for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
1183     if (*tmpa > *tmpb) {
1184       return MP_GT;
1185     }
1186
1187     if (*tmpa < *tmpb) {
1188       return MP_LT;
1189     }
1190   }
1191   return MP_EQ;
1192 }
1193
1194
1195 /* compare two ints (signed)*/
1196 int
1197 mp_cmp (mp_int * a, mp_int * b)
1198 {
1199   /* compare based on sign */
1200   if (a->sign != b->sign) {
1201      if (a->sign == MP_NEG) {
1202         return MP_LT;
1203      } else {
1204         return MP_GT;
1205      }
1206   }
1207
1208   /* compare digits */
1209   if (a->sign == MP_NEG) {
1210      /* if negative compare opposite direction */
1211      return mp_cmp_mag(b, a);
1212   } else {
1213      return mp_cmp_mag(a, b);
1214   }
1215 }
1216
1217
1218 /* compare a digit */
1219 int mp_cmp_d(mp_int * a, mp_digit b)
1220 {
1221   /* compare based on sign */
1222   if (a->sign == MP_NEG) {
1223     return MP_LT;
1224   }
1225
1226   /* compare based on magnitude */
1227   if (a->used > 1) {
1228     return MP_GT;
1229   }
1230
1231   /* compare the only digit of a to b */
1232   if (a->dp[0] > b) {
1233     return MP_GT;
1234   } else if (a->dp[0] < b) {
1235     return MP_LT;
1236   } else {
1237     return MP_EQ;
1238   }
1239 }
1240
1241
1242 /* set to a digit */
1243 void mp_set (mp_int * a, mp_digit b)
1244 {
1245   mp_zero (a);
1246   a->dp[0] = b & MP_MASK;
1247   a->used  = (a->dp[0] != 0) ? 1 : 0;
1248 }
1249
1250
1251 /* c = a mod b, 0 <= c < b */
1252 int
1253 mp_mod (mp_int * a, mp_int * b, mp_int * c)
1254 {
1255   mp_int  t;
1256   int     res;
1257
1258   if ((res = mp_init (&t)) != MP_OKAY) {
1259     return res;
1260   }
1261
1262   if ((res = mp_div (a, b, NULL, &t)) != MP_OKAY) {
1263     mp_clear (&t);
1264     return res;
1265   }
1266
1267   if (t.sign != b->sign) {
1268     res = mp_add (b, &t, c);
1269   } else {
1270     res = MP_OKAY;
1271     mp_exch (&t, c);
1272   }
1273
1274   mp_clear (&t);
1275   return res;
1276 }
1277
1278
1279 /* slower bit-bang division... also smaller */
1280 int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d)
1281 {
1282    mp_int ta, tb, tq, q;
1283    int    res, n, n2;
1284
1285   /* is divisor zero ? */
1286   if (mp_iszero (b) == 1) {
1287     return MP_VAL;
1288   }
1289
1290   /* if a < b then q=0, r = a */
1291   if (mp_cmp_mag (a, b) == MP_LT) {
1292     if (d != NULL) {
1293       res = mp_copy (a, d);
1294     } else {
1295       res = MP_OKAY;
1296     }
1297     if (c != NULL) {
1298       mp_zero (c);
1299     }
1300     return res;
1301   }
1302
1303   /* init our temps */
1304   if ((res = mp_init_multi(&ta, &tb, &tq, &q, 0, 0)) != MP_OKAY) {
1305      return res;
1306   }
1307
1308
1309   mp_set(&tq, 1);
1310   n = mp_count_bits(a) - mp_count_bits(b);
1311   if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
1312       ((res = mp_abs(b, &tb)) != MP_OKAY) ||
1313       ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
1314       ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
1315       goto LBL_ERR;
1316   }
1317
1318   while (n-- >= 0) {
1319      if (mp_cmp(&tb, &ta) != MP_GT) {
1320         if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
1321             ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
1322            goto LBL_ERR;
1323         }
1324      }
1325      if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
1326          ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
1327            goto LBL_ERR;
1328      }
1329   }
1330
1331   /* now q == quotient and ta == remainder */
1332   n  = a->sign;
1333   n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
1334   if (c != NULL) {
1335      mp_exch(c, &q);
1336      c->sign  = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
1337   }
1338   if (d != NULL) {
1339      mp_exch(d, &ta);
1340      d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
1341   }
1342 LBL_ERR:
1343    mp_clear(&ta);
1344    mp_clear(&tb);
1345    mp_clear(&tq);
1346    mp_clear(&q);
1347    return res;
1348 }
1349
1350
1351 /* b = a/2 */
1352 int mp_div_2(mp_int * a, mp_int * b)
1353 {
1354   int     x, res, oldused;
1355
1356   /* copy */
1357   if (b->alloc < a->used) {
1358     if ((res = mp_grow (b, a->used)) != MP_OKAY) {
1359       return res;
1360     }
1361   }
1362
1363   oldused = b->used;
1364   b->used = a->used;
1365   {
1366     register mp_digit r, rr, *tmpa, *tmpb;
1367
1368     /* source alias */
1369     tmpa = a->dp + b->used - 1;
1370
1371     /* dest alias */
1372     tmpb = b->dp + b->used - 1;
1373
1374     /* carry */
1375     r = 0;
1376     for (x = b->used - 1; x >= 0; x--) {
1377       /* get the carry for the next iteration */
1378       rr = *tmpa & 1;
1379
1380       /* shift the current digit, add in carry and store */
1381       *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1382
1383       /* forward carry to next iteration */
1384       r = rr;
1385     }
1386
1387     /* zero excess digits */
1388     tmpb = b->dp + b->used;
1389     for (x = b->used; x < oldused; x++) {
1390       *tmpb++ = 0;
1391     }
1392   }
1393   b->sign = a->sign;
1394   mp_clamp (b);
1395   return MP_OKAY;
1396 }
1397
1398
1399 /* high level addition (handles signs) */
1400 int mp_add (mp_int * a, mp_int * b, mp_int * c)
1401 {
1402   int     sa, sb, res;
1403
1404   /* get sign of both inputs */
1405   sa = a->sign;
1406   sb = b->sign;
1407
1408   /* handle two cases, not four */
1409   if (sa == sb) {
1410     /* both positive or both negative */
1411     /* add their magnitudes, copy the sign */
1412     c->sign = sa;
1413     res = s_mp_add (a, b, c);
1414   } else {
1415     /* one positive, the other negative */
1416     /* subtract the one with the greater magnitude from */
1417     /* the one of the lesser magnitude.  The result gets */
1418     /* the sign of the one with the greater magnitude. */
1419     if (mp_cmp_mag (a, b) == MP_LT) {
1420       c->sign = sb;
1421       res = s_mp_sub (b, a, c);
1422     } else {
1423       c->sign = sa;
1424       res = s_mp_sub (a, b, c);
1425     }
1426   }
1427   return res;
1428 }
1429
1430
1431 /* low level addition, based on HAC pp.594, Algorithm 14.7 */
1432 int
1433 s_mp_add (mp_int * a, mp_int * b, mp_int * c)
1434 {
1435   mp_int *x;
1436   int     olduse, res, min, max;
1437
1438   /* find sizes, we let |a| <= |b| which means we have to sort
1439    * them.  "x" will point to the input with the most digits
1440    */
1441   if (a->used > b->used) {
1442     min = b->used;
1443     max = a->used;
1444     x = a;
1445   } else {
1446     min = a->used;
1447     max = b->used;
1448     x = b;
1449   }
1450
1451   /* init result */
1452   if (c->alloc < max + 1) {
1453     if ((res = mp_grow (c, max + 1)) != MP_OKAY) {
1454       return res;
1455     }
1456   }
1457
1458   /* get old used digit count and set new one */
1459   olduse = c->used;
1460   c->used = max + 1;
1461
1462   {
1463     register mp_digit u, *tmpa, *tmpb, *tmpc;
1464     register int i;
1465
1466     /* alias for digit pointers */
1467
1468     /* first input */
1469     tmpa = a->dp;
1470
1471     /* second input */
1472     tmpb = b->dp;
1473
1474     /* destination */
1475     tmpc = c->dp;
1476
1477     /* zero the carry */
1478     u = 0;
1479     for (i = 0; i < min; i++) {
1480       /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
1481       *tmpc = *tmpa++ + *tmpb++ + u;
1482
1483       /* U = carry bit of T[i] */
1484       u = *tmpc >> ((mp_digit)DIGIT_BIT);
1485
1486       /* take away carry bit from T[i] */
1487       *tmpc++ &= MP_MASK;
1488     }
1489
1490     /* now copy higher words if any, that is in A+B
1491      * if A or B has more digits add those in
1492      */
1493     if (min != max) {
1494       for (; i < max; i++) {
1495         /* T[i] = X[i] + U */
1496         *tmpc = x->dp[i] + u;
1497
1498         /* U = carry bit of T[i] */
1499         u = *tmpc >> ((mp_digit)DIGIT_BIT);
1500
1501         /* take away carry bit from T[i] */
1502         *tmpc++ &= MP_MASK;
1503       }
1504     }
1505
1506     /* add carry */
1507     *tmpc++ = u;
1508
1509     /* clear digits above oldused */
1510     for (i = c->used; i < olduse; i++) {
1511       *tmpc++ = 0;
1512     }
1513   }
1514
1515   mp_clamp (c);
1516   return MP_OKAY;
1517 }
1518
1519
1520 /* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
1521 int
1522 s_mp_sub (mp_int * a, mp_int * b, mp_int * c)
1523 {
1524   int     olduse, res, min, max;
1525
1526   /* find sizes */
1527   min = b->used;
1528   max = a->used;
1529
1530   /* init result */
1531   if (c->alloc < max) {
1532     if ((res = mp_grow (c, max)) != MP_OKAY) {
1533       return res;
1534     }
1535   }
1536   olduse = c->used;
1537   c->used = max;
1538
1539   {
1540     register mp_digit u, *tmpa, *tmpb, *tmpc;
1541     register int i;
1542
1543     /* alias for digit pointers */
1544     tmpa = a->dp;
1545     tmpb = b->dp;
1546     tmpc = c->dp;
1547
1548     /* set carry to zero */
1549     u = 0;
1550     for (i = 0; i < min; i++) {
1551       /* T[i] = A[i] - B[i] - U */
1552       *tmpc = *tmpa++ - *tmpb++ - u;
1553
1554       /* U = carry bit of T[i]
1555        * Note this saves performing an AND operation since
1556        * if a carry does occur it will propagate all the way to the
1557        * MSB.  As a result a single shift is enough to get the carry
1558        */
1559       u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
1560
1561       /* Clear carry from T[i] */
1562       *tmpc++ &= MP_MASK;
1563     }
1564
1565     /* now copy higher words if any, e.g. if A has more digits than B  */
1566     for (; i < max; i++) {
1567       /* T[i] = A[i] - U */
1568       *tmpc = *tmpa++ - u;
1569
1570       /* U = carry bit of T[i] */
1571       u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
1572
1573       /* Clear carry from T[i] */
1574       *tmpc++ &= MP_MASK;
1575     }
1576
1577     /* clear digits above used (since we may not have grown result above) */
1578     for (i = c->used; i < olduse; i++) {
1579       *tmpc++ = 0;
1580     }
1581   }
1582
1583   mp_clamp (c);
1584   return MP_OKAY;
1585 }
1586
1587
1588 /* high level subtraction (handles signs) */
1589 int
1590 mp_sub (mp_int * a, mp_int * b, mp_int * c)
1591 {
1592   int     sa, sb, res;
1593
1594   sa = a->sign;
1595   sb = b->sign;
1596
1597   if (sa != sb) {
1598     /* subtract a negative from a positive, OR */
1599     /* subtract a positive from a negative. */
1600     /* In either case, ADD their magnitudes, */
1601     /* and use the sign of the first number. */
1602     c->sign = sa;
1603     res = s_mp_add (a, b, c);
1604   } else {
1605     /* subtract a positive from a positive, OR */
1606     /* subtract a negative from a negative. */
1607     /* First, take the difference between their */
1608     /* magnitudes, then... */
1609     if (mp_cmp_mag (a, b) != MP_LT) {
1610       /* Copy the sign from the first */
1611       c->sign = sa;
1612       /* The first has a larger or equal magnitude */
1613       res = s_mp_sub (a, b, c);
1614     } else {
1615       /* The result has the *opposite* sign from */
1616       /* the first number. */
1617       c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
1618       /* The second has a larger magnitude */
1619       res = s_mp_sub (b, a, c);
1620     }
1621   }
1622   return res;
1623 }
1624
1625
1626 /* determines if reduce_2k_l can be used */
1627 int mp_reduce_is_2k_l(mp_int *a)
1628 {
1629    int ix, iy;
1630
1631    if (a->used == 0) {
1632       return MP_NO;
1633    } else if (a->used == 1) {
1634       return MP_YES;
1635    } else if (a->used > 1) {
1636       /* if more than half of the digits are -1 we're sold */
1637       for (iy = ix = 0; ix < a->used; ix++) {
1638           if (a->dp[ix] == MP_MASK) {
1639               ++iy;
1640           }
1641       }
1642       return (iy >= (a->used/2)) ? MP_YES : MP_NO;
1643
1644    }
1645    return MP_NO;
1646 }
1647
1648
1649 /* determines if mp_reduce_2k can be used */
1650 int mp_reduce_is_2k(mp_int *a)
1651 {
1652    int ix, iy, iw;
1653    mp_digit iz;
1654
1655    if (a->used == 0) {
1656       return MP_NO;
1657    } else if (a->used == 1) {
1658       return MP_YES;
1659    } else if (a->used > 1) {
1660       iy = mp_count_bits(a);
1661       iz = 1;
1662       iw = 1;
1663
1664       /* Test every bit from the second digit up, must be 1 */
1665       for (ix = DIGIT_BIT; ix < iy; ix++) {
1666           if ((a->dp[iw] & iz) == 0) {
1667              return MP_NO;
1668           }
1669           iz <<= 1;
1670           if (iz > (mp_digit)MP_MASK) {
1671              ++iw;
1672              iz = 1;
1673           }
1674       }
1675    }
1676    return MP_YES;
1677 }
1678
1679
1680 /* determines if a number is a valid DR modulus */
1681 int mp_dr_is_modulus(mp_int *a)
1682 {
1683    int ix;
1684
1685    /* must be at least two digits */
1686    if (a->used < 2) {
1687       return 0;
1688    }
1689
1690    /* must be of the form b**k - a [a <= b] so all
1691     * but the first digit must be equal to -1 (mod b).
1692     */
1693    for (ix = 1; ix < a->used; ix++) {
1694        if (a->dp[ix] != MP_MASK) {
1695           return 0;
1696        }
1697    }
1698    return 1;
1699 }
1700
1701
1702 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
1703  *
1704  * Uses a left-to-right k-ary sliding window to compute the modular
1705  * exponentiation.
1706  * The value of k changes based on the size of the exponent.
1707  *
1708  * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
1709  */
1710
1711 #ifdef MP_LOW_MEM
1712    #define TAB_SIZE 32
1713 #else
1714    #define TAB_SIZE 256
1715 #endif
1716
1717 int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y,
1718                      int redmode)
1719 {
1720   mp_int  M[TAB_SIZE], res;
1721   mp_digit buf, mp;
1722   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
1723
1724   /* use a pointer to the reduction algorithm.  This allows us to use
1725    * one of many reduction algorithms without modding the guts of
1726    * the code with if statements everywhere.
1727    */
1728   int     (*redux)(mp_int*,mp_int*,mp_digit);
1729
1730   /* find window size */
1731   x = mp_count_bits (X);
1732   if (x <= 7) {
1733     winsize = 2;
1734   } else if (x <= 36) {
1735     winsize = 3;
1736   } else if (x <= 140) {
1737     winsize = 4;
1738   } else if (x <= 450) {
1739     winsize = 5;
1740   } else if (x <= 1303) {
1741     winsize = 6;
1742   } else if (x <= 3529) {
1743     winsize = 7;
1744   } else {
1745     winsize = 8;
1746   }
1747
1748 #ifdef MP_LOW_MEM
1749   if (winsize > 5) {
1750      winsize = 5;
1751   }
1752 #endif
1753
1754   /* init M array */
1755   /* init first cell */
1756   if ((err = mp_init(&M[1])) != MP_OKAY) {
1757      return err;
1758   }
1759
1760   /* now init the second half of the array */
1761   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
1762     if ((err = mp_init(&M[x])) != MP_OKAY) {
1763       for (y = 1<<(winsize-1); y < x; y++) {
1764         mp_clear (&M[y]);
1765       }
1766       mp_clear(&M[1]);
1767       return err;
1768     }
1769   }
1770
1771   /* determine and setup reduction code */
1772   if (redmode == 0) {
1773 #ifdef BN_MP_MONTGOMERY_SETUP_C
1774      /* now setup montgomery  */
1775      if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
1776         goto LBL_M;
1777      }
1778 #else
1779      err = MP_VAL;
1780      goto LBL_M;
1781 #endif
1782
1783      /* automatically pick the comba one if available (saves quite a few
1784         calls/ifs) */
1785 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
1786      if (((P->used * 2 + 1) < MP_WARRAY) &&
1787           P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
1788         redux = fast_mp_montgomery_reduce;
1789      } else
1790 #endif
1791      {
1792 #ifdef BN_MP_MONTGOMERY_REDUCE_C
1793         /* use slower baseline Montgomery method */
1794         redux = mp_montgomery_reduce;
1795 #else
1796         err = MP_VAL;
1797         goto LBL_M;
1798 #endif
1799      }
1800   } else if (redmode == 1) {
1801 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
1802      /* setup DR reduction for moduli of the form B**k - b */
1803      mp_dr_setup(P, &mp);
1804      redux = mp_dr_reduce;
1805 #else
1806      err = MP_VAL;
1807      goto LBL_M;
1808 #endif
1809   } else {
1810 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
1811      /* setup DR reduction for moduli of the form 2**k - b */
1812      if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
1813         goto LBL_M;
1814      }
1815      redux = mp_reduce_2k;
1816 #else
1817      err = MP_VAL;
1818      goto LBL_M;
1819 #endif
1820   }
1821
1822   /* setup result */
1823   if ((err = mp_init (&res)) != MP_OKAY) {
1824     goto LBL_M;
1825   }
1826
1827   /* create M table
1828    *
1829
1830    *
1831    * The first half of the table is not computed though accept for M[0] and M[1]
1832    */
1833
1834   if (redmode == 0) {
1835 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
1836      /* now we need R mod m */
1837      if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
1838        goto LBL_RES;
1839      }
1840 #else
1841      err = MP_VAL;
1842      goto LBL_RES;
1843 #endif
1844
1845      /* now set M[1] to G * R mod m */
1846      if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
1847        goto LBL_RES;
1848      }
1849   } else {
1850      mp_set(&res, 1);
1851      if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
1852         goto LBL_RES;
1853      }
1854   }
1855
1856   /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times*/
1857   if ((err = mp_copy (&M[1], &M[(mp_digit)(1 << (winsize - 1))])) != MP_OKAY) {
1858     goto LBL_RES;
1859   }
1860
1861   for (x = 0; x < (winsize - 1); x++) {
1862     if ((err = mp_sqr (&M[(mp_digit)(1 << (winsize - 1))], &M[(mp_digit)(1 << (winsize - 1))])) != MP_OKAY) {
1863       goto LBL_RES;
1864     }
1865     if ((err = redux (&M[(mp_digit)(1 << (winsize - 1))], P, mp)) != MP_OKAY) {
1866       goto LBL_RES;
1867     }
1868   }
1869
1870   /* create upper table */
1871   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
1872     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
1873       goto LBL_RES;
1874     }
1875     if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
1876       goto LBL_RES;
1877     }
1878   }
1879
1880   /* set initial mode and bit cnt */
1881   mode   = 0;
1882   bitcnt = 1;
1883   buf    = 0;
1884   digidx = X->used - 1;
1885   bitcpy = 0;
1886   bitbuf = 0;
1887
1888   for (;;) {
1889     /* grab next digit as required */
1890     if (--bitcnt == 0) {
1891       /* if digidx == -1 we are out of digits so break */
1892       if (digidx == -1) {
1893         break;
1894       }
1895       /* read next digit and reset bitcnt */
1896       buf    = X->dp[digidx--];
1897       bitcnt = (int)DIGIT_BIT;
1898     }
1899
1900     /* grab the next msb from the exponent */
1901     y     = (int)(buf >> (DIGIT_BIT - 1)) & 1;
1902     buf <<= (mp_digit)1;
1903
1904     /* if the bit is zero and mode == 0 then we ignore it
1905      * These represent the leading zero bits before the first 1 bit
1906      * in the exponent.  Technically this opt is not required but it
1907      * does lower the # of trivial squaring/reductions used
1908      */
1909     if (mode == 0 && y == 0) {
1910       continue;
1911     }
1912
1913     /* if the bit is zero and mode == 1 then we square */
1914     if (mode == 1 && y == 0) {
1915       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
1916         goto LBL_RES;
1917       }
1918       if ((err = redux (&res, P, mp)) != MP_OKAY) {
1919         goto LBL_RES;
1920       }
1921       continue;
1922     }
1923
1924     /* else we add it to the window */
1925     bitbuf |= (y << (winsize - ++bitcpy));
1926     mode    = 2;
1927
1928     if (bitcpy == winsize) {
1929       /* ok window is filled so square as required and multiply  */
1930       /* square first */
1931       for (x = 0; x < winsize; x++) {
1932         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
1933           goto LBL_RES;
1934         }
1935         if ((err = redux (&res, P, mp)) != MP_OKAY) {
1936           goto LBL_RES;
1937         }
1938       }
1939
1940       /* then multiply */
1941       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
1942         goto LBL_RES;
1943       }
1944       if ((err = redux (&res, P, mp)) != MP_OKAY) {
1945         goto LBL_RES;
1946       }
1947
1948       /* empty window and reset */
1949       bitcpy = 0;
1950       bitbuf = 0;
1951       mode   = 1;
1952     }
1953   }
1954
1955   /* if bits remain then square/multiply */
1956   if (mode == 2 && bitcpy > 0) {
1957     /* square then multiply if the bit is set */
1958     for (x = 0; x < bitcpy; x++) {
1959       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
1960         goto LBL_RES;
1961       }
1962       if ((err = redux (&res, P, mp)) != MP_OKAY) {
1963         goto LBL_RES;
1964       }
1965
1966       /* get next bit of the window */
1967       bitbuf <<= 1;
1968       if ((bitbuf & (1 << winsize)) != 0) {
1969         /* then multiply */
1970         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
1971           goto LBL_RES;
1972         }
1973         if ((err = redux (&res, P, mp)) != MP_OKAY) {
1974           goto LBL_RES;
1975         }
1976       }
1977     }
1978   }
1979
1980   if (redmode == 0) {
1981      /* fixup result if Montgomery reduction is used
1982       * recall that any value in a Montgomery system is
1983       * actually multiplied by R mod n.  So we have
1984       * to reduce one more time to cancel out the factor
1985       * of R.
1986       */
1987      if ((err = redux(&res, P, mp)) != MP_OKAY) {
1988        goto LBL_RES;
1989      }
1990   }
1991
1992   /* swap res with Y */
1993   mp_exch (&res, Y);
1994   err = MP_OKAY;
1995 LBL_RES:mp_clear (&res);
1996 LBL_M:
1997   mp_clear(&M[1]);
1998   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
1999     mp_clear (&M[x]);
2000   }
2001   return err;
2002 }
2003
2004
2005 /* setups the montgomery reduction stuff */
2006 int
2007 mp_montgomery_setup (mp_int * n, mp_digit * rho)
2008 {
2009   mp_digit x, b;
2010
2011 /* fast inversion mod 2**k
2012  *
2013  * Based on the fact that
2014  *
2015  * XA = 1 (mod 2**n)  =>  (X(2-XA)) A = 1 (mod 2**2n)
2016  *                    =>  2*X*A - X*X*A*A = 1
2017  *                    =>  2*(1) - (1)     = 1
2018  */
2019   b = n->dp[0];
2020
2021   if ((b & 1) == 0) {
2022     return MP_VAL;
2023   }
2024
2025   x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
2026   x *= 2 - b * x;               /* here x*a==1 mod 2**8 */
2027 #if !defined(MP_8BIT)
2028   x *= 2 - b * x;               /* here x*a==1 mod 2**16 */
2029 #endif
2030 #if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT))
2031   x *= 2 - b * x;               /* here x*a==1 mod 2**32 */
2032 #endif
2033 #ifdef MP_64BIT
2034   x *= 2 - b * x;               /* here x*a==1 mod 2**64 */
2035 #endif
2036
2037   /* rho = -1/m mod b */
2038   /* TAO, switched mp_word casts to mp_digit to shut up compiler */
2039   *rho = (((mp_digit)1 << ((mp_digit) DIGIT_BIT)) - x) & MP_MASK;
2040
2041   return MP_OKAY;
2042 }
2043
2044
2045 /* computes xR**-1 == x (mod N) via Montgomery Reduction
2046  *
2047  * This is an optimized implementation of montgomery_reduce
2048  * which uses the comba method to quickly calculate the columns of the
2049  * reduction.
2050  *
2051  * Based on Algorithm 14.32 on pp.601 of HAC.
2052 */
2053 int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
2054 {
2055   int     ix, res, olduse;
2056 #ifdef CYASSL_SMALL_STACK
2057   mp_word* W;    /* uses dynamic memory and slower */
2058 #else
2059   mp_word W[MP_WARRAY];
2060 #endif
2061
2062   /* get old used count */
2063   olduse = x->used;
2064
2065   /* grow a as required */
2066   if (x->alloc < n->used + 1) {
2067     if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) {
2068       return res;
2069     }
2070   }
2071
2072 #ifdef CYASSL_SMALL_STACK
2073   W = (mp_word*)XMALLOC(sizeof(mp_word) * MP_WARRAY, 0, DYNAMIC_TYPE_BIGINT);
2074   if (W == NULL)
2075     return MP_MEM;
2076 #endif
2077
2078   /* first we have to get the digits of the input into
2079    * an array of double precision words W[...]
2080    */
2081   {
2082     register mp_word *_W;
2083     register mp_digit *tmpx;
2084
2085     /* alias for the W[] array */
2086     _W   = W;
2087
2088     /* alias for the digits of  x*/
2089     tmpx = x->dp;
2090
2091     /* copy the digits of a into W[0..a->used-1] */
2092     for (ix = 0; ix < x->used; ix++) {
2093       *_W++ = *tmpx++;
2094     }
2095
2096     /* zero the high words of W[a->used..m->used*2] */
2097     for (; ix < n->used * 2 + 1; ix++) {
2098       *_W++ = 0;
2099     }
2100   }
2101
2102   /* now we proceed to zero successive digits
2103    * from the least significant upwards
2104    */
2105   for (ix = 0; ix < n->used; ix++) {
2106     /* mu = ai * m' mod b
2107      *
2108      * We avoid a double precision multiplication (which isn't required)
2109      * by casting the value down to a mp_digit.  Note this requires
2110      * that W[ix-1] have  the carry cleared (see after the inner loop)
2111      */
2112     register mp_digit mu;
2113     mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
2114
2115     /* a = a + mu * m * b**i
2116      *
2117      * This is computed in place and on the fly.  The multiplication
2118      * by b**i is handled by offseting which columns the results
2119      * are added to.
2120      *
2121      * Note the comba method normally doesn't handle carries in the
2122      * inner loop In this case we fix the carry from the previous
2123      * column since the Montgomery reduction requires digits of the
2124      * result (so far) [see above] to work.  This is
2125      * handled by fixing up one carry after the inner loop.  The
2126      * carry fixups are done in order so after these loops the
2127      * first m->used words of W[] have the carries fixed
2128      */
2129     {
2130       register int iy;
2131       register mp_digit *tmpn;
2132       register mp_word *_W;
2133
2134       /* alias for the digits of the modulus */
2135       tmpn = n->dp;
2136
2137       /* Alias for the columns set by an offset of ix */
2138       _W = W + ix;
2139
2140       /* inner loop */
2141       for (iy = 0; iy < n->used; iy++) {
2142           *_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
2143       }
2144     }
2145
2146     /* now fix carry for next digit, W[ix+1] */
2147     W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
2148   }
2149
2150   /* now we have to propagate the carries and
2151    * shift the words downward [all those least
2152    * significant digits we zeroed].
2153    */
2154   {
2155     register mp_digit *tmpx;
2156     register mp_word *_W, *_W1;
2157
2158     /* nox fix rest of carries */
2159
2160     /* alias for current word */
2161     _W1 = W + ix;
2162
2163     /* alias for next word, where the carry goes */
2164     _W = W + ++ix;
2165
2166     for (; ix <= n->used * 2 + 1; ix++) {
2167       *_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
2168     }
2169
2170     /* copy out, A = A/b**n
2171      *
2172      * The result is A/b**n but instead of converting from an
2173      * array of mp_word to mp_digit than calling mp_rshd
2174      * we just copy them in the right order
2175      */
2176
2177     /* alias for destination word */
2178     tmpx = x->dp;
2179
2180     /* alias for shifted double precision result */
2181     _W = W + n->used;
2182
2183     for (ix = 0; ix < n->used + 1; ix++) {
2184       *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
2185     }
2186
2187     /* zero oldused digits, if the input a was larger than
2188      * m->used+1 we'll have to clear the digits
2189      */
2190     for (; ix < olduse; ix++) {
2191       *tmpx++ = 0;
2192     }
2193   }
2194
2195   /* set the max used and clamp */
2196   x->used = n->used + 1;
2197   mp_clamp (x);
2198
2199 #ifdef CYASSL_SMALL_STACK
2200   XFREE(W, 0, DYNAMIC_TYPE_BIGINT);
2201 #endif
2202
2203   /* if A >= m then A = A - m */
2204   if (mp_cmp_mag (x, n) != MP_LT) {
2205     return s_mp_sub (x, n, x);
2206   }
2207   return MP_OKAY;
2208 }
2209
2210
2211 /* computes xR**-1 == x (mod N) via Montgomery Reduction */
2212 int
2213 mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
2214 {
2215   int     ix, res, digs;
2216   mp_digit mu;
2217
2218   /* can the fast reduction [comba] method be used?
2219    *
2220    * Note that unlike in mul you're safely allowed *less*
2221    * than the available columns [255 per default] since carries
2222    * are fixed up in the inner loop.
2223    */
2224   digs = n->used * 2 + 1;
2225   if ((digs < MP_WARRAY) &&
2226       n->used <
2227       (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
2228     return fast_mp_montgomery_reduce (x, n, rho);
2229   }
2230
2231   /* grow the input as required */
2232   if (x->alloc < digs) {
2233     if ((res = mp_grow (x, digs)) != MP_OKAY) {
2234       return res;
2235     }
2236   }
2237   x->used = digs;
2238
2239   for (ix = 0; ix < n->used; ix++) {
2240     /* mu = ai * rho mod b
2241      *
2242      * The value of rho must be precalculated via
2243      * montgomery_setup() such that
2244      * it equals -1/n0 mod b this allows the
2245      * following inner loop to reduce the
2246      * input one digit at a time
2247      */
2248     mu = (mp_digit) (((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK);
2249
2250     /* a = a + mu * m * b**i */
2251     {
2252       register int iy;
2253       register mp_digit *tmpn, *tmpx, u;
2254       register mp_word r;
2255
2256       /* alias for digits of the modulus */
2257       tmpn = n->dp;
2258
2259       /* alias for the digits of x [the input] */
2260       tmpx = x->dp + ix;
2261
2262       /* set the carry to zero */
2263       u = 0;
2264
2265       /* Multiply and add in place */
2266       for (iy = 0; iy < n->used; iy++) {
2267         /* compute product and sum */
2268         r       = ((mp_word)mu) * ((mp_word)*tmpn++) +
2269                   ((mp_word) u) + ((mp_word) * tmpx);
2270
2271         /* get carry */
2272         u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
2273
2274         /* fix digit */
2275         *tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK));
2276       }
2277       /* At this point the ix'th digit of x should be zero */
2278
2279
2280       /* propagate carries upwards as required*/
2281       while (u) {
2282         *tmpx   += u;
2283         u        = *tmpx >> DIGIT_BIT;
2284         *tmpx++ &= MP_MASK;
2285       }
2286     }
2287   }
2288
2289   /* at this point the n.used'th least
2290    * significant digits of x are all zero
2291    * which means we can shift x to the
2292    * right by n.used digits and the
2293    * residue is unchanged.
2294    */
2295
2296   /* x = x/b**n.used */
2297   mp_clamp(x);
2298   mp_rshd (x, n->used);
2299
2300   /* if x >= n then x = x - n */
2301   if (mp_cmp_mag (x, n) != MP_LT) {
2302     return s_mp_sub (x, n, x);
2303   }
2304
2305   return MP_OKAY;
2306 }
2307
2308
2309 /* determines the setup value */
2310 void mp_dr_setup(mp_int *a, mp_digit *d)
2311 {
2312    /* the casts are required if DIGIT_BIT is one less than
2313     * the number of bits in a mp_digit [e.g. DIGIT_BIT==31]
2314     */
2315    *d = (mp_digit)((((mp_word)1) << ((mp_word)DIGIT_BIT)) -
2316         ((mp_word)a->dp[0]));
2317 }
2318
2319
2320 /* reduce "x" in place modulo "n" using the Diminished Radix algorithm.
2321  *
2322  * Based on algorithm from the paper
2323  *
2324  * "Generating Efficient Primes for Discrete Log Cryptosystems"
2325  *                 Chae Hoon Lim, Pil Joong Lee,
2326  *          POSTECH Information Research Laboratories
2327  *
2328  * The modulus must be of a special format [see manual]
2329  *
2330  * Has been modified to use algorithm 7.10 from the LTM book instead
2331  *
2332  * Input x must be in the range 0 <= x <= (n-1)**2
2333  */
2334 int
2335 mp_dr_reduce (mp_int * x, mp_int * n, mp_digit k)
2336 {
2337   int      err, i, m;
2338   mp_word  r;
2339   mp_digit mu, *tmpx1, *tmpx2;
2340
2341   /* m = digits in modulus */
2342   m = n->used;
2343
2344   /* ensure that "x" has at least 2m digits */
2345   if (x->alloc < m + m) {
2346     if ((err = mp_grow (x, m + m)) != MP_OKAY) {
2347       return err;
2348     }
2349   }
2350
2351 /* top of loop, this is where the code resumes if
2352  * another reduction pass is required.
2353  */
2354 top:
2355   /* aliases for digits */
2356   /* alias for lower half of x */
2357   tmpx1 = x->dp;
2358
2359   /* alias for upper half of x, or x/B**m */
2360   tmpx2 = x->dp + m;
2361
2362   /* set carry to zero */
2363   mu = 0;
2364
2365   /* compute (x mod B**m) + k * [x/B**m] inline and inplace */
2366   for (i = 0; i < m; i++) {
2367       r         = ((mp_word)*tmpx2++) * ((mp_word)k) + *tmpx1 + mu;
2368       *tmpx1++  = (mp_digit)(r & MP_MASK);
2369       mu        = (mp_digit)(r >> ((mp_word)DIGIT_BIT));
2370   }
2371
2372   /* set final carry */
2373   *tmpx1++ = mu;
2374
2375   /* zero words above m */
2376   for (i = m + 1; i < x->used; i++) {
2377       *tmpx1++ = 0;
2378   }
2379
2380   /* clamp, sub and return */
2381   mp_clamp (x);
2382
2383   /* if x >= n then subtract and reduce again
2384    * Each successive "recursion" makes the input smaller and smaller.
2385    */
2386   if (mp_cmp_mag (x, n) != MP_LT) {
2387     s_mp_sub(x, n, x);
2388     goto top;
2389   }
2390   return MP_OKAY;
2391 }
2392
2393
2394 /* reduces a modulo n where n is of the form 2**p - d */
2395 int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d)
2396 {
2397    mp_int q;
2398    int    p, res;
2399
2400    if ((res = mp_init(&q)) != MP_OKAY) {
2401       return res;
2402    }
2403
2404    p = mp_count_bits(n);
2405 top:
2406    /* q = a/2**p, a = a mod 2**p */
2407    if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
2408       goto ERR;
2409    }
2410
2411    if (d != 1) {
2412       /* q = q * d */
2413       if ((res = mp_mul_d(&q, d, &q)) != MP_OKAY) {
2414          goto ERR;
2415       }
2416    }
2417
2418    /* a = a + q */
2419    if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
2420       goto ERR;
2421    }
2422
2423    if (mp_cmp_mag(a, n) != MP_LT) {
2424       s_mp_sub(a, n, a);
2425       goto top;
2426    }
2427
2428 ERR:
2429    mp_clear(&q);
2430    return res;
2431 }
2432
2433
2434 /* determines the setup value */
2435 int mp_reduce_2k_setup(mp_int *a, mp_digit *d)
2436 {
2437    int res, p;
2438    mp_int tmp;
2439
2440    if ((res = mp_init(&tmp)) != MP_OKAY) {
2441       return res;
2442    }
2443
2444    p = mp_count_bits(a);
2445    if ((res = mp_2expt(&tmp, p)) != MP_OKAY) {
2446       mp_clear(&tmp);
2447       return res;
2448    }
2449
2450    if ((res = s_mp_sub(&tmp, a, &tmp)) != MP_OKAY) {
2451       mp_clear(&tmp);
2452       return res;
2453    }
2454
2455    *d = tmp.dp[0];
2456    mp_clear(&tmp);
2457    return MP_OKAY;
2458 }
2459
2460
2461 /* computes a = 2**b
2462  *
2463  * Simple algorithm which zeroes the int, grows it then just sets one bit
2464  * as required.
2465  */
2466 int
2467 mp_2expt (mp_int * a, int b)
2468 {
2469   int     res;
2470
2471   /* zero a as per default */
2472   mp_zero (a);
2473
2474   /* grow a to accomodate the single bit */
2475   if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
2476     return res;
2477   }
2478
2479   /* set the used count of where the bit will go */
2480   a->used = b / DIGIT_BIT + 1;
2481
2482   /* put the single bit in its place */
2483   a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT);
2484
2485   return MP_OKAY;
2486 }
2487
2488
2489 /* multiply by a digit */
2490 int
2491 mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
2492 {
2493   mp_digit u, *tmpa, *tmpc;
2494   mp_word  r;
2495   int      ix, res, olduse;
2496
2497   /* make sure c is big enough to hold a*b */
2498   if (c->alloc < a->used + 1) {
2499     if ((res = mp_grow (c, a->used + 1)) != MP_OKAY) {
2500       return res;
2501     }
2502   }
2503
2504   /* get the original destinations used count */
2505   olduse = c->used;
2506
2507   /* set the sign */
2508   c->sign = a->sign;
2509
2510   /* alias for a->dp [source] */
2511   tmpa = a->dp;
2512
2513   /* alias for c->dp [dest] */
2514   tmpc = c->dp;
2515
2516   /* zero carry */
2517   u = 0;
2518
2519   /* compute columns */
2520   for (ix = 0; ix < a->used; ix++) {
2521     /* compute product and carry sum for this term */
2522     r       = ((mp_word) u) + ((mp_word)*tmpa++) * ((mp_word)b);
2523
2524     /* mask off higher bits to get a single digit */
2525     *tmpc++ = (mp_digit) (r & ((mp_word) MP_MASK));
2526
2527     /* send carry into next iteration */
2528     u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
2529   }
2530
2531   /* store final carry [if any] and increment ix offset  */
2532   *tmpc++ = u;
2533   ++ix;
2534
2535   /* now zero digits above the top */
2536   while (ix++ < olduse) {
2537      *tmpc++ = 0;
2538   }
2539
2540   /* set used count */
2541   c->used = a->used + 1;
2542   mp_clamp(c);
2543
2544   return MP_OKAY;
2545 }
2546
2547
2548 /* d = a * b (mod c) */
2549 int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
2550 {
2551   int     res;
2552   mp_int  t;
2553
2554   if ((res = mp_init (&t)) != MP_OKAY) {
2555     return res;
2556   }
2557
2558   if ((res = mp_mul (a, b, &t)) != MP_OKAY) {
2559     mp_clear (&t);
2560     return res;
2561   }
2562   res = mp_mod (&t, c, d);
2563   mp_clear (&t);
2564   return res;
2565 }
2566
2567
2568 /* computes b = a*a */
2569 int
2570 mp_sqr (mp_int * a, mp_int * b)
2571 {
2572   int     res;
2573
2574   {
2575 #ifdef BN_FAST_S_MP_SQR_C
2576     /* can we use the fast comba multiplier? */
2577     if ((a->used * 2 + 1) < MP_WARRAY &&
2578          a->used <
2579          (1 << (sizeof(mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) {
2580       res = fast_s_mp_sqr (a, b);
2581     } else
2582 #endif
2583 #ifdef BN_S_MP_SQR_C
2584       res = s_mp_sqr (a, b);
2585 #else
2586       res = MP_VAL;
2587 #endif
2588   }
2589   b->sign = MP_ZPOS;
2590   return res;
2591 }
2592
2593
2594 /* high level multiplication (handles sign) */
2595 int mp_mul (mp_int * a, mp_int * b, mp_int * c)
2596 {
2597   int     res, neg;
2598   neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
2599
2600   {
2601     /* can we use the fast multiplier?
2602      *
2603      * The fast multiplier can be used if the output will
2604      * have less than MP_WARRAY digits and the number of
2605      * digits won't affect carry propagation
2606      */
2607     int     digs = a->used + b->used + 1;
2608
2609 #ifdef BN_FAST_S_MP_MUL_DIGS_C
2610     if ((digs < MP_WARRAY) &&
2611         MIN(a->used, b->used) <=
2612         (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
2613       res = fast_s_mp_mul_digs (a, b, c, digs);
2614     } else
2615 #endif
2616 #ifdef BN_S_MP_MUL_DIGS_C
2617       res = s_mp_mul (a, b, c); /* uses s_mp_mul_digs */
2618 #else
2619       res = MP_VAL;
2620 #endif
2621
2622   }
2623   c->sign = (c->used > 0) ? neg : MP_ZPOS;
2624   return res;
2625 }
2626
2627
2628 /* b = a*2 */
2629 int mp_mul_2(mp_int * a, mp_int * b)
2630 {
2631   int     x, res, oldused;
2632
2633   /* grow to accomodate result */
2634   if (b->alloc < a->used + 1) {
2635     if ((res = mp_grow (b, a->used + 1)) != MP_OKAY) {
2636       return res;
2637     }
2638   }
2639
2640   oldused = b->used;
2641   b->used = a->used;
2642
2643   {
2644     register mp_digit r, rr, *tmpa, *tmpb;
2645
2646     /* alias for source */
2647     tmpa = a->dp;
2648
2649     /* alias for dest */
2650     tmpb = b->dp;
2651
2652     /* carry */
2653     r = 0;
2654     for (x = 0; x < a->used; x++) {
2655
2656       /* get what will be the *next* carry bit from the
2657        * MSB of the current digit
2658        */
2659       rr = *tmpa >> ((mp_digit)(DIGIT_BIT - 1));
2660
2661       /* now shift up this digit, add in the carry [from the previous] */
2662       *tmpb++ = ((*tmpa++ << ((mp_digit)1)) | r) & MP_MASK;
2663
2664       /* copy the carry that would be from the source
2665        * digit into the next iteration
2666        */
2667       r = rr;
2668     }
2669
2670     /* new leading digit? */
2671     if (r != 0) {
2672       /* add a MSB which is always 1 at this point */
2673       *tmpb = 1;
2674       ++(b->used);
2675     }
2676
2677     /* now zero any excess digits on the destination
2678      * that we didn't write to
2679      */
2680     tmpb = b->dp + b->used;
2681     for (x = b->used; x < oldused; x++) {
2682       *tmpb++ = 0;
2683     }
2684   }
2685   b->sign = a->sign;
2686   return MP_OKAY;
2687 }
2688
2689
2690 /* divide by three (based on routine from MPI and the GMP manual) */
2691 int
2692 mp_div_3 (mp_int * a, mp_int *c, mp_digit * d)
2693 {
2694   mp_int   q;
2695   mp_word  w, t;
2696   mp_digit b;
2697   int      res, ix;
2698
2699   /* b = 2**DIGIT_BIT / 3 */
2700   b = (((mp_word)1) << ((mp_word)DIGIT_BIT)) / ((mp_word)3);
2701
2702   if ((res = mp_init_size(&q, a->used)) != MP_OKAY) {
2703      return res;
2704   }
2705
2706   q.used = a->used;
2707   q.sign = a->sign;
2708   w = 0;
2709   for (ix = a->used - 1; ix >= 0; ix--) {
2710      w = (w << ((mp_word)DIGIT_BIT)) | ((mp_word)a->dp[ix]);
2711
2712      if (w >= 3) {
2713         /* multiply w by [1/3] */
2714         t = (w * ((mp_word)b)) >> ((mp_word)DIGIT_BIT);
2715
2716         /* now subtract 3 * [w/3] from w, to get the remainder */
2717         w -= t+t+t;
2718
2719         /* fixup the remainder as required since
2720          * the optimization is not exact.
2721          */
2722         while (w >= 3) {
2723            t += 1;
2724            w -= 3;
2725         }
2726       } else {
2727         t = 0;
2728       }
2729       q.dp[ix] = (mp_digit)t;
2730   }
2731
2732   /* [optional] store the remainder */
2733   if (d != NULL) {
2734      *d = (mp_digit)w;
2735   }
2736
2737   /* [optional] store the quotient */
2738   if (c != NULL) {
2739      mp_clamp(&q);
2740      mp_exch(&q, c);
2741   }
2742   mp_clear(&q);
2743
2744   return res;
2745 }
2746
2747
2748 /* init an mp_init for a given size */
2749 int mp_init_size (mp_int * a, int size)
2750 {
2751   int x;
2752
2753   /* pad size so there are always extra digits */
2754   size += (MP_PREC * 2) - (size % MP_PREC);
2755
2756   /* alloc mem */
2757   a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * size, 0,
2758                                       DYNAMIC_TYPE_BIGINT);
2759   if (a->dp == NULL) {
2760     return MP_MEM;
2761   }
2762
2763   /* set the members */
2764   a->used  = 0;
2765   a->alloc = size;
2766   a->sign  = MP_ZPOS;
2767
2768   /* zero the digits */
2769   for (x = 0; x < size; x++) {
2770       a->dp[x] = 0;
2771   }
2772
2773   return MP_OKAY;
2774 }
2775
2776
2777 /* the jist of squaring...
2778  * you do like mult except the offset of the tmpx [one that
2779  * starts closer to zero] can't equal the offset of tmpy.
2780  * So basically you set up iy like before then you min it with
2781  * (ty-tx) so that it never happens.  You double all those
2782  * you add in the inner loop
2783
2784 After that loop you do the squares and add them in.
2785 */
2786
2787 int fast_s_mp_sqr (mp_int * a, mp_int * b)
2788 {
2789   int       olduse, res, pa, ix, iz;
2790 #ifdef CYASSL_SMALL_STACK
2791   mp_digit* W;    /* uses dynamic memory and slower */
2792 #else
2793   mp_digit W[MP_WARRAY];
2794 #endif
2795   mp_digit  *tmpx;
2796   mp_word   W1;
2797
2798   /* grow the destination as required */
2799   pa = a->used + a->used;
2800   if (b->alloc < pa) {
2801     if ((res = mp_grow (b, pa)) != MP_OKAY) {
2802       return res;
2803     }
2804   }
2805
2806   if (pa > MP_WARRAY)
2807     return MP_RANGE;  /* TAO range check */
2808
2809 #ifdef CYASSL_SMALL_STACK
2810   W = (mp_digit*)XMALLOC(sizeof(mp_digit) * MP_WARRAY, 0, DYNAMIC_TYPE_BIGINT);
2811   if (W == NULL)
2812     return MP_MEM;
2813 #endif
2814
2815   /* number of output digits to produce */
2816   W1 = 0;
2817   for (ix = 0; ix < pa; ix++) {
2818       int      tx, ty, iy;
2819       mp_word  _W;
2820       mp_digit *tmpy;
2821
2822       /* clear counter */
2823       _W = 0;
2824
2825       /* get offsets into the two bignums */
2826       ty = MIN(a->used-1, ix);
2827       tx = ix - ty;
2828
2829       /* setup temp aliases */
2830       tmpx = a->dp + tx;
2831       tmpy = a->dp + ty;
2832
2833       /* this is the number of times the loop will iterrate, essentially
2834          while (tx++ < a->used && ty-- >= 0) { ... }
2835        */
2836       iy = MIN(a->used-tx, ty+1);
2837
2838       /* now for squaring tx can never equal ty
2839        * we halve the distance since they approach at a rate of 2x
2840        * and we have to round because odd cases need to be executed
2841        */
2842       iy = MIN(iy, (ty-tx+1)>>1);
2843
2844       /* execute loop */
2845       for (iz = 0; iz < iy; iz++) {
2846          _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
2847       }
2848
2849       /* double the inner product and add carry */
2850       _W = _W + _W + W1;
2851
2852       /* even columns have the square term in them */
2853       if ((ix&1) == 0) {
2854          _W += ((mp_word)a->dp[ix>>1])*((mp_word)a->dp[ix>>1]);
2855       }
2856
2857       /* store it */
2858       W[ix] = (mp_digit)(_W & MP_MASK);
2859
2860       /* make next carry */
2861       W1 = _W >> ((mp_word)DIGIT_BIT);
2862   }
2863
2864   /* setup dest */
2865   olduse  = b->used;
2866   b->used = a->used+a->used;
2867
2868   {
2869     mp_digit *tmpb;
2870     tmpb = b->dp;
2871     for (ix = 0; ix < pa; ix++) {
2872       *tmpb++ = W[ix] & MP_MASK;
2873     }
2874
2875     /* clear unused digits [that existed in the old copy of c] */
2876     for (; ix < olduse; ix++) {
2877       *tmpb++ = 0;
2878     }
2879   }
2880   mp_clamp (b);
2881
2882 #ifdef CYASSL_SMALL_STACK
2883   XFREE(W, 0, DYNAMIC_TYPE_BIGINT);
2884 #endif
2885
2886   return MP_OKAY;
2887 }
2888
2889
2890 /* Fast (comba) multiplier
2891  *
2892  * This is the fast column-array [comba] multiplier.  It is
2893  * designed to compute the columns of the product first
2894  * then handle the carries afterwards.  This has the effect
2895  * of making the nested loops that compute the columns very
2896  * simple and schedulable on super-scalar processors.
2897  *
2898  * This has been modified to produce a variable number of
2899  * digits of output so if say only a half-product is required
2900  * you don't have to compute the upper half (a feature
2901  * required for fast Barrett reduction).
2902  *
2903  * Based on Algorithm 14.12 on pp.595 of HAC.
2904  *
2905  */
2906 int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
2907 {
2908   int     olduse, res, pa, ix, iz;
2909 #ifdef CYASSL_SMALL_STACK
2910   mp_digit* W;    /* uses dynamic memory and slower */
2911 #else
2912   mp_digit W[MP_WARRAY];
2913 #endif
2914   register mp_word  _W;
2915
2916   /* grow the destination as required */
2917   if (c->alloc < digs) {
2918     if ((res = mp_grow (c, digs)) != MP_OKAY) {
2919       return res;
2920     }
2921   }
2922
2923   /* number of output digits to produce */
2924   pa = MIN(digs, a->used + b->used);
2925   if (pa > MP_WARRAY)
2926     return MP_RANGE;  /* TAO range check */
2927
2928 #ifdef CYASSL_SMALL_STACK
2929   W = (mp_digit*)XMALLOC(sizeof(mp_digit) * MP_WARRAY, 0, DYNAMIC_TYPE_BIGINT);
2930   if (W == NULL)
2931     return MP_MEM;
2932 #endif
2933
2934   /* clear the carry */
2935   _W = 0;
2936   for (ix = 0; ix < pa; ix++) {
2937       int      tx, ty;
2938       int      iy;
2939       mp_digit *tmpx, *tmpy;
2940
2941       /* get offsets into the two bignums */
2942       ty = MIN(b->used-1, ix);
2943       tx = ix - ty;
2944
2945       /* setup temp aliases */
2946       tmpx = a->dp + tx;
2947       tmpy = b->dp + ty;
2948
2949       /* this is the number of times the loop will iterrate, essentially
2950          while (tx++ < a->used && ty-- >= 0) { ... }
2951        */
2952       iy = MIN(a->used-tx, ty+1);
2953
2954       /* execute loop */
2955       for (iz = 0; iz < iy; ++iz) {
2956          _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
2957
2958       }
2959
2960       /* store term */
2961       W[ix] = ((mp_digit)_W) & MP_MASK;
2962
2963       /* make next carry */
2964       _W = _W >> ((mp_word)DIGIT_BIT);
2965  }
2966
2967   /* setup dest */
2968   olduse  = c->used;
2969   c->used = pa;
2970
2971   {
2972     register mp_digit *tmpc;
2973     tmpc = c->dp;
2974     for (ix = 0; ix < pa+1; ix++) {
2975       /* now extract the previous digit [below the carry] */
2976       *tmpc++ = W[ix];
2977     }
2978
2979     /* clear unused digits [that existed in the old copy of c] */
2980     for (; ix < olduse; ix++) {
2981       *tmpc++ = 0;
2982     }
2983   }
2984   mp_clamp (c);
2985
2986 #ifdef CYASSL_SMALL_STACK
2987   XFREE(W, 0, DYNAMIC_TYPE_BIGINT);
2988 #endif
2989
2990   return MP_OKAY;
2991 }
2992
2993
2994 /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
2995 int s_mp_sqr (mp_int * a, mp_int * b)
2996 {
2997   mp_int  t;
2998   int     res, ix, iy, pa;
2999   mp_word r;
3000   mp_digit u, tmpx, *tmpt;
3001
3002   pa = a->used;
3003   if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
3004     return res;
3005   }
3006
3007   /* default used is maximum possible size */
3008   t.used = 2*pa + 1;
3009
3010   for (ix = 0; ix < pa; ix++) {
3011     /* first calculate the digit at 2*ix */
3012     /* calculate double precision result */
3013     r = ((mp_word) t.dp[2*ix]) +
3014         ((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);
3015
3016     /* store lower part in result */
3017     t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));
3018
3019     /* get the carry */
3020     u           = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
3021
3022     /* left hand side of A[ix] * A[iy] */
3023     tmpx        = a->dp[ix];
3024
3025     /* alias for where to store the results */
3026     tmpt        = t.dp + (2*ix + 1);
3027
3028     for (iy = ix + 1; iy < pa; iy++) {
3029       /* first calculate the product */
3030       r       = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
3031
3032       /* now calculate the double precision result, note we use
3033        * addition instead of *2 since it's easier to optimize
3034        */
3035       r       = ((mp_word) *tmpt) + r + r + ((mp_word) u);
3036
3037       /* store lower part */
3038       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
3039
3040       /* get carry */
3041       u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
3042     }
3043     /* propagate upwards */
3044     while (u != ((mp_digit) 0)) {
3045       r       = ((mp_word) *tmpt) + ((mp_word) u);
3046       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
3047       u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
3048     }
3049   }
3050
3051   mp_clamp (&t);
3052   mp_exch (&t, b);
3053   mp_clear (&t);
3054   return MP_OKAY;
3055 }
3056
3057
3058 /* multiplies |a| * |b| and only computes upto digs digits of result
3059  * HAC pp. 595, Algorithm 14.12  Modified so you can control how
3060  * many digits of output are created.
3061  */
3062 int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
3063 {
3064   mp_int  t;
3065   int     res, pa, pb, ix, iy;
3066   mp_digit u;
3067   mp_word r;
3068   mp_digit tmpx, *tmpt, *tmpy;
3069
3070   /* can we use the fast multiplier? */
3071   if (((digs) < MP_WARRAY) &&
3072       MIN (a->used, b->used) <
3073           (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
3074     return fast_s_mp_mul_digs (a, b, c, digs);
3075   }
3076
3077   if ((res = mp_init_size (&t, digs)) != MP_OKAY) {
3078     return res;
3079   }
3080   t.used = digs;
3081
3082   /* compute the digits of the product directly */
3083   pa = a->used;
3084   for (ix = 0; ix < pa; ix++) {
3085     /* set the carry to zero */
3086     u = 0;
3087
3088     /* limit ourselves to making digs digits of output */
3089     pb = MIN (b->used, digs - ix);
3090
3091     /* setup some aliases */
3092     /* copy of the digit from a used within the nested loop */
3093     tmpx = a->dp[ix];
3094
3095     /* an alias for the destination shifted ix places */
3096     tmpt = t.dp + ix;
3097
3098     /* an alias for the digits of b */
3099     tmpy = b->dp;
3100
3101     /* compute the columns of the output and propagate the carry */
3102     for (iy = 0; iy < pb; iy++) {
3103       /* compute the column as a mp_word */
3104       r       = ((mp_word)*tmpt) +
3105                 ((mp_word)tmpx) * ((mp_word)*tmpy++) +
3106                 ((mp_word) u);
3107
3108       /* the new column is the lower part of the result */
3109       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
3110
3111       /* get the carry word from the result */
3112       u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
3113     }
3114     /* set carry if it is placed below digs */
3115     if (ix + iy < digs) {
3116       *tmpt = u;
3117     }
3118   }
3119
3120   mp_clamp (&t);
3121   mp_exch (&t, c);
3122
3123   mp_clear (&t);
3124   return MP_OKAY;
3125 }
3126
3127
3128 /*
3129  * shifts with subtractions when the result is greater than b.
3130  *
3131  * The method is slightly modified to shift B unconditionally upto just under
3132  * the leading bit of b.  This saves alot of multiple precision shifting.
3133  */
3134 int mp_montgomery_calc_normalization (mp_int * a, mp_int * b)
3135 {
3136   int     x, bits, res;
3137
3138   /* how many bits of last digit does b use */
3139   bits = mp_count_bits (b) % DIGIT_BIT;
3140
3141   if (b->used > 1) {
3142      if ((res = mp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) != MP_OKAY) {
3143         return res;
3144      }
3145   } else {
3146      mp_set(a, 1);
3147      bits = 1;
3148   }
3149
3150
3151   /* now compute C = A * B mod b */
3152   for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
3153     if ((res = mp_mul_2 (a, a)) != MP_OKAY) {
3154       return res;
3155     }
3156     if (mp_cmp_mag (a, b) != MP_LT) {
3157       if ((res = s_mp_sub (a, b, a)) != MP_OKAY) {
3158         return res;
3159       }
3160     }
3161   }
3162
3163   return MP_OKAY;
3164 }
3165
3166
3167 #ifdef MP_LOW_MEM
3168    #define TAB_SIZE 32
3169 #else
3170    #define TAB_SIZE 256
3171 #endif
3172
3173 int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
3174 {
3175   mp_int  M[TAB_SIZE], res, mu;
3176   mp_digit buf;
3177   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
3178   int (*redux)(mp_int*,mp_int*,mp_int*);
3179
3180   /* find window size */
3181   x = mp_count_bits (X);
3182   if (x <= 7) {
3183     winsize = 2;
3184   } else if (x <= 36) {
3185     winsize = 3;
3186   } else if (x <= 140) {
3187     winsize = 4;
3188   } else if (x <= 450) {
3189     winsize = 5;
3190   } else if (x <= 1303) {
3191     winsize = 6;
3192   } else if (x <= 3529) {
3193     winsize = 7;
3194   } else {
3195     winsize = 8;
3196   }
3197
3198 #ifdef MP_LOW_MEM
3199     if (winsize > 5) {
3200        winsize = 5;
3201     }
3202 #endif
3203
3204   /* init M array */
3205   /* init first cell */
3206   if ((err = mp_init(&M[1])) != MP_OKAY) {
3207      return err;
3208   }
3209
3210   /* now init the second half of the array */
3211   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
3212     if ((err = mp_init(&M[x])) != MP_OKAY) {
3213       for (y = 1<<(winsize-1); y < x; y++) {
3214         mp_clear (&M[y]);
3215       }
3216       mp_clear(&M[1]);
3217       return err;
3218     }
3219   }
3220
3221   /* create mu, used for Barrett reduction */
3222   if ((err = mp_init (&mu)) != MP_OKAY) {
3223     goto LBL_M;
3224   }
3225
3226   if (redmode == 0) {
3227      if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
3228         goto LBL_MU;
3229      }
3230      redux = mp_reduce;
3231   } else {
3232      if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
3233         goto LBL_MU;
3234      }
3235      redux = mp_reduce_2k_l;
3236   }
3237
3238   /* create M table
3239    *
3240    * The M table contains powers of the base,
3241    * e.g. M[x] = G**x mod P
3242    *
3243    * The first half of the table is not
3244    * computed though accept for M[0] and M[1]
3245    */
3246   if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
3247     goto LBL_MU;
3248   }
3249
3250   /* compute the value at M[1<<(winsize-1)] by squaring
3251    * M[1] (winsize-1) times
3252    */
3253   if ((err = mp_copy (&M[1], &M[(mp_digit)(1 << (winsize - 1))])) != MP_OKAY) {
3254     goto LBL_MU;
3255   }
3256
3257   for (x = 0; x < (winsize - 1); x++) {
3258     /* square it */
3259     if ((err = mp_sqr (&M[(mp_digit)(1 << (winsize - 1))],
3260                        &M[(mp_digit)(1 << (winsize - 1))])) != MP_OKAY) {
3261       goto LBL_MU;
3262     }
3263
3264     /* reduce modulo P */
3265     if ((err = redux (&M[(mp_digit)(1 << (winsize - 1))], P, &mu)) != MP_OKAY) {
3266       goto LBL_MU;
3267     }
3268   }
3269
3270   /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
3271    * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
3272    */
3273   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
3274     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
3275       goto LBL_MU;
3276     }
3277     if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
3278       goto LBL_MU;
3279     }
3280   }
3281
3282   /* setup result */
3283   if ((err = mp_init (&res)) != MP_OKAY) {
3284     goto LBL_MU;
3285   }
3286   mp_set (&res, 1);
3287
3288   /* set initial mode and bit cnt */
3289   mode   = 0;
3290   bitcnt = 1;
3291   buf    = 0;
3292   digidx = X->used - 1;
3293   bitcpy = 0;
3294   bitbuf = 0;
3295
3296   for (;;) {
3297     /* grab next digit as required */
3298     if (--bitcnt == 0) {
3299       /* if digidx == -1 we are out of digits */
3300       if (digidx == -1) {
3301         break;
3302       }
3303       /* read next digit and reset the bitcnt */
3304       buf    = X->dp[digidx--];
3305       bitcnt = (int) DIGIT_BIT;
3306     }
3307
3308     /* grab the next msb from the exponent */
3309     y     = (int)(buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
3310     buf <<= (mp_digit)1;
3311
3312     /* if the bit is zero and mode == 0 then we ignore it
3313      * These represent the leading zero bits before the first 1 bit
3314      * in the exponent.  Technically this opt is not required but it
3315      * does lower the # of trivial squaring/reductions used
3316      */
3317     if (mode == 0 && y == 0) {
3318       continue;
3319     }
3320
3321     /* if the bit is zero and mode == 1 then we square */
3322     if (mode == 1 && y == 0) {
3323       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
3324         goto LBL_RES;
3325       }
3326       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
3327         goto LBL_RES;
3328       }
3329       continue;
3330     }
3331
3332     /* else we add it to the window */
3333     bitbuf |= (y << (winsize - ++bitcpy));
3334     mode    = 2;
3335
3336     if (bitcpy == winsize) {
3337       /* ok window is filled so square as required and multiply  */
3338       /* square first */
3339       for (x = 0; x < winsize; x++) {
3340         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
3341           goto LBL_RES;
3342         }
3343         if ((err = redux (&res, P, &mu)) != MP_OKAY) {
3344           goto LBL_RES;
3345         }
3346       }
3347
3348       /* then multiply */
3349       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
3350         goto LBL_RES;
3351       }
3352       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
3353         goto LBL_RES;
3354       }
3355
3356       /* empty window and reset */
3357       bitcpy = 0;
3358       bitbuf = 0;
3359       mode   = 1;
3360     }
3361   }
3362
3363   /* if bits remain then square/multiply */
3364   if (mode == 2 && bitcpy > 0) {
3365     /* square then multiply if the bit is set */
3366     for (x = 0; x < bitcpy; x++) {
3367       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
3368         goto LBL_RES;
3369       }
3370       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
3371         goto LBL_RES;
3372       }
3373
3374       bitbuf <<= 1;
3375       if ((bitbuf & (1 << winsize)) != 0) {
3376         /* then multiply */
3377         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
3378           goto LBL_RES;
3379         }
3380         if ((err = redux (&res, P, &mu)) != MP_OKAY) {
3381           goto LBL_RES;
3382         }
3383       }
3384     }
3385   }
3386
3387   mp_exch (&res, Y);
3388   err = MP_OKAY;
3389 LBL_RES:mp_clear (&res);
3390 LBL_MU:mp_clear (&mu);
3391 LBL_M:
3392   mp_clear(&M[1]);
3393   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
3394     mp_clear (&M[x]);
3395   }
3396   return err;
3397 }
3398
3399
3400 /* pre-calculate the value required for Barrett reduction
3401  * For a given modulus "b" it calulates the value required in "a"
3402  */
3403 int mp_reduce_setup (mp_int * a, mp_int * b)
3404 {
3405   int     res;
3406
3407   if ((res = mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) {
3408     return res;
3409   }
3410   return mp_div (a, b, a, NULL);
3411 }
3412
3413
3414 /* reduces x mod m, assumes 0 < x < m**2, mu is
3415  * precomputed via mp_reduce_setup.
3416  * From HAC pp.604 Algorithm 14.42
3417  */
3418 int mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
3419 {
3420   mp_int  q;
3421   int     res, um = m->used;
3422
3423   /* q = x */
3424   if ((res = mp_init_copy (&q, x)) != MP_OKAY) {
3425     return res;
3426   }
3427
3428   /* q1 = x / b**(k-1)  */
3429   mp_rshd (&q, um - 1);
3430
3431   /* according to HAC this optimization is ok */
3432   if (((mp_word) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
3433     if ((res = mp_mul (&q, mu, &q)) != MP_OKAY) {
3434       goto CLEANUP;
3435     }
3436   } else {
3437 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
3438     if ((res = s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
3439       goto CLEANUP;
3440     }
3441 #elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C)
3442     if ((res = fast_s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
3443       goto CLEANUP;
3444     }
3445 #else
3446     {
3447       res = MP_VAL;
3448       goto CLEANUP;
3449     }
3450 #endif
3451   }
3452
3453   /* q3 = q2 / b**(k+1) */
3454   mp_rshd (&q, um + 1);
3455
3456   /* x = x mod b**(k+1), quick (no division) */
3457   if ((res = mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
3458     goto CLEANUP;
3459   }
3460
3461   /* q = q * m mod b**(k+1), quick (no division) */
3462   if ((res = s_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) {
3463     goto CLEANUP;
3464   }
3465
3466   /* x = x - q */
3467   if ((res = mp_sub (x, &q, x)) != MP_OKAY) {
3468     goto CLEANUP;
3469   }
3470
3471   /* If x < 0, add b**(k+1) to it */
3472   if (mp_cmp_d (x, 0) == MP_LT) {
3473     mp_set (&q, 1);
3474     if ((res = mp_lshd (&q, um + 1)) != MP_OKAY)
3475       goto CLEANUP;
3476     if ((res = mp_add (x, &q, x)) != MP_OKAY)
3477       goto CLEANUP;
3478   }
3479
3480   /* Back off if it's too big */
3481   while (mp_cmp (x, m) != MP_LT) {
3482     if ((res = s_mp_sub (x, m, x)) != MP_OKAY) {
3483       goto CLEANUP;
3484     }
3485   }
3486
3487 CLEANUP:
3488   mp_clear (&q);
3489
3490   return res;
3491 }
3492
3493
3494 /* reduces a modulo n where n is of the form 2**p - d
3495    This differs from reduce_2k since "d" can be larger
3496    than a single digit.
3497 */
3498 int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d)
3499 {
3500    mp_int q;
3501    int    p, res;
3502
3503    if ((res = mp_init(&q)) != MP_OKAY) {
3504       return res;
3505    }
3506
3507    p = mp_count_bits(n);
3508 top:
3509    /* q = a/2**p, a = a mod 2**p */
3510    if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
3511       goto ERR;
3512    }
3513
3514    /* q = q * d */
3515    if ((res = mp_mul(&q, d, &q)) != MP_OKAY) {
3516       goto ERR;
3517    }
3518
3519    /* a = a + q */
3520    if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
3521       goto ERR;
3522    }
3523
3524    if (mp_cmp_mag(a, n) != MP_LT) {
3525       s_mp_sub(a, n, a);
3526       goto top;
3527    }
3528
3529 ERR:
3530    mp_clear(&q);
3531    return res;
3532 }
3533
3534
3535 /* determines the setup value */
3536 int mp_reduce_2k_setup_l(mp_int *a, mp_int *d)
3537 {
3538    int    res;
3539    mp_int tmp;
3540
3541    if ((res = mp_init(&tmp)) != MP_OKAY) {
3542       return res;
3543    }
3544
3545    if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
3546       goto ERR;
3547    }
3548
3549    if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) {
3550       goto ERR;
3551    }
3552
3553 ERR:
3554    mp_clear(&tmp);
3555    return res;
3556 }
3557
3558
3559 /* multiplies |a| * |b| and does not compute the lower digs digits
3560  * [meant to get the higher part of the product]
3561  */
3562 int
3563 s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
3564 {
3565   mp_int  t;
3566   int     res, pa, pb, ix, iy;
3567   mp_digit u;
3568   mp_word r;
3569   mp_digit tmpx, *tmpt, *tmpy;
3570
3571   /* can we use the fast multiplier? */
3572 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
3573   if (((a->used + b->used + 1) < MP_WARRAY)
3574       && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
3575     return fast_s_mp_mul_high_digs (a, b, c, digs);
3576   }
3577 #endif
3578
3579   if ((res = mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) {
3580     return res;
3581   }
3582   t.used = a->used + b->used + 1;
3583
3584   pa = a->used;
3585   pb = b->used;
3586   for (ix = 0; ix < pa; ix++) {
3587     /* clear the carry */
3588     u = 0;
3589
3590     /* left hand side of A[ix] * B[iy] */
3591     tmpx = a->dp[ix];
3592
3593     /* alias to the address of where the digits will be stored */
3594     tmpt = &(t.dp[digs]);
3595
3596     /* alias for where to read the right hand side from */
3597     tmpy = b->dp + (digs - ix);
3598
3599     for (iy = digs - ix; iy < pb; iy++) {
3600       /* calculate the double precision result */
3601       r       = ((mp_word)*tmpt) +
3602                 ((mp_word)tmpx) * ((mp_word)*tmpy++) +
3603                 ((mp_word) u);
3604
3605       /* get the lower part */
3606       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
3607
3608       /* carry the carry */
3609       u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
3610     }
3611     *tmpt = u;
3612   }
3613   mp_clamp (&t);
3614   mp_exch (&t, c);
3615   mp_clear (&t);
3616   return MP_OKAY;
3617 }
3618
3619
3620 /* this is a modified version of fast_s_mul_digs that only produces
3621  * output digits *above* digs.  See the comments for fast_s_mul_digs
3622  * to see how it works.
3623  *
3624  * This is used in the Barrett reduction since for one of the multiplications
3625  * only the higher digits were needed.  This essentially halves the work.
3626  *
3627  * Based on Algorithm 14.12 on pp.595 of HAC.
3628  */
3629 int fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
3630 {
3631   int     olduse, res, pa, ix, iz;
3632 #ifdef CYASSL_SMALL_STACK
3633   mp_digit* W;    /* uses dynamic memory and slower */
3634 #else
3635   mp_digit W[MP_WARRAY];
3636 #endif
3637   mp_word  _W;
3638
3639   /* grow the destination as required */
3640   pa = a->used + b->used;
3641   if (c->alloc < pa) {
3642     if ((res = mp_grow (c, pa)) != MP_OKAY) {
3643       return res;
3644     }
3645   }
3646
3647   if (pa > MP_WARRAY)
3648     return MP_RANGE;  /* TAO range check */
3649
3650 #ifdef CYASSL_SMALL_STACK
3651   W = (mp_digit*)XMALLOC(sizeof(mp_digit) * MP_WARRAY, 0, DYNAMIC_TYPE_BIGINT);
3652   if (W == NULL)
3653     return MP_MEM;
3654 #endif
3655
3656   /* number of output digits to produce */
3657   pa = a->used + b->used;
3658   _W = 0;
3659   for (ix = digs; ix < pa; ix++) {
3660       int      tx, ty, iy;
3661       mp_digit *tmpx, *tmpy;
3662
3663       /* get offsets into the two bignums */
3664       ty = MIN(b->used-1, ix);
3665       tx = ix - ty;
3666
3667       /* setup temp aliases */
3668       tmpx = a->dp + tx;
3669       tmpy = b->dp + ty;
3670
3671       /* this is the number of times the loop will iterrate, essentially its
3672          while (tx++ < a->used && ty-- >= 0) { ... }
3673        */
3674       iy = MIN(a->used-tx, ty+1);
3675
3676       /* execute loop */
3677       for (iz = 0; iz < iy; iz++) {
3678          _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
3679       }
3680
3681       /* store term */
3682       W[ix] = ((mp_digit)_W) & MP_MASK;
3683
3684       /* make next carry */
3685       _W = _W >> ((mp_word)DIGIT_BIT);
3686   }
3687
3688   /* setup dest */
3689   olduse  = c->used;
3690   c->used = pa;
3691
3692   {
3693     register mp_digit *tmpc;
3694
3695     tmpc = c->dp + digs;
3696     for (ix = digs; ix <= pa; ix++) {
3697       /* now extract the previous digit [below the carry] */
3698       *tmpc++ = W[ix];
3699     }
3700
3701     /* clear unused digits [that existed in the old copy of c] */
3702     for (; ix < olduse; ix++) {
3703       *tmpc++ = 0;
3704     }
3705   }
3706   mp_clamp (c);
3707
3708 #ifdef CYASSL_SMALL_STACK
3709   XFREE(W, 0, DYNAMIC_TYPE_BIGINT);
3710 #endif
3711
3712   return MP_OKAY;
3713 }
3714
3715
3716 /* set a 32-bit const */
3717 int mp_set_int (mp_int * a, unsigned long b)
3718 {
3719   int     x, res;
3720
3721   mp_zero (a);
3722
3723   /* set four bits at a time */
3724   for (x = 0; x < 8; x++) {
3725     /* shift the number up four bits */
3726     if ((res = mp_mul_2d (a, 4, a)) != MP_OKAY) {
3727       return res;
3728     }
3729
3730     /* OR in the top four bits of the source */
3731     a->dp[0] |= (b >> 28) & 15;
3732
3733     /* shift the source up to the next four bits */
3734     b <<= 4;
3735
3736     /* ensure that digits are not clamped off */
3737     a->used += 1;
3738   }
3739   mp_clamp (a);
3740   return MP_OKAY;
3741 }
3742
3743
3744 #if defined(CYASSL_KEY_GEN) || defined(HAVE_ECC)
3745
3746 /* c = a * a (mod b) */
3747 int mp_sqrmod (mp_int * a, mp_int * b, mp_int * c)
3748 {
3749   int     res;
3750   mp_int  t;
3751
3752   if ((res = mp_init (&t)) != MP_OKAY) {
3753     return res;
3754   }
3755
3756   if ((res = mp_sqr (a, &t)) != MP_OKAY) {
3757     mp_clear (&t);
3758     return res;
3759   }
3760   res = mp_mod (&t, b, c);
3761   mp_clear (&t);
3762   return res;
3763 }
3764
3765 #endif
3766
3767
3768 #if defined(HAVE_ECC) || !defined(NO_PWDBASED) || defined(CYASSL_SNIFFER) || defined(CYASSL_HAVE_WOLFSCEP) || defined(CYASSL_KEY_GEN)
3769
3770 /* single digit addition */
3771 int mp_add_d (mp_int* a, mp_digit b, mp_int* c)
3772 {
3773   int     res, ix, oldused;
3774   mp_digit *tmpa, *tmpc, mu;
3775
3776   /* grow c as required */
3777   if (c->alloc < a->used + 1) {
3778      if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) {
3779         return res;
3780      }
3781   }
3782
3783   /* if a is negative and |a| >= b, call c = |a| - b */
3784   if (a->sign == MP_NEG && (a->used > 1 || a->dp[0] >= b)) {
3785      /* temporarily fix sign of a */
3786      a->sign = MP_ZPOS;
3787
3788      /* c = |a| - b */
3789      res = mp_sub_d(a, b, c);
3790
3791      /* fix sign  */
3792      a->sign = c->sign = MP_NEG;
3793
3794      /* clamp */
3795      mp_clamp(c);
3796
3797      return res;
3798   }
3799
3800   /* old number of used digits in c */
3801   oldused = c->used;
3802
3803   /* sign always positive */
3804   c->sign = MP_ZPOS;
3805
3806   /* source alias */
3807   tmpa    = a->dp;
3808
3809   /* destination alias */
3810   tmpc    = c->dp;
3811
3812   /* if a is positive */
3813   if (a->sign == MP_ZPOS) {
3814      /* add digit, after this we're propagating
3815       * the carry.
3816       */
3817      *tmpc   = *tmpa++ + b;
3818      mu      = *tmpc >> DIGIT_BIT;
3819      *tmpc++ &= MP_MASK;
3820
3821      /* now handle rest of the digits */
3822      for (ix = 1; ix < a->used; ix++) {
3823         *tmpc   = *tmpa++ + mu;
3824         mu      = *tmpc >> DIGIT_BIT;
3825         *tmpc++ &= MP_MASK;
3826      }
3827      /* set final carry */
3828      if (mu != 0 && ix < c->alloc) {
3829         ix++;
3830         *tmpc++  = mu;
3831      }
3832
3833      /* setup size */
3834      c->used = a->used + 1;
3835   } else {
3836      /* a was negative and |a| < b */
3837      c->used  = 1;
3838
3839      /* the result is a single digit */
3840      if (a->used == 1) {
3841         *tmpc++  =  b - a->dp[0];
3842      } else {
3843         *tmpc++  =  b;
3844      }
3845
3846      /* setup count so the clearing of oldused
3847       * can fall through correctly
3848       */
3849      ix       = 1;
3850   }
3851
3852   /* now zero to oldused */
3853   while (ix++ < oldused) {
3854      *tmpc++ = 0;
3855   }
3856   mp_clamp(c);
3857
3858   return MP_OKAY;
3859 }
3860
3861
3862 /* single digit subtraction */
3863 int mp_sub_d (mp_int * a, mp_digit b, mp_int * c)
3864 {
3865   mp_digit *tmpa, *tmpc, mu;
3866   int       res, ix, oldused;
3867
3868   /* grow c as required */
3869   if (c->alloc < a->used + 1) {
3870      if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) {
3871         return res;
3872      }
3873   }
3874
3875   /* if a is negative just do an unsigned
3876    * addition [with fudged signs]
3877    */
3878   if (a->sign == MP_NEG) {
3879      a->sign = MP_ZPOS;
3880      res     = mp_add_d(a, b, c);
3881      a->sign = c->sign = MP_NEG;
3882
3883      /* clamp */
3884      mp_clamp(c);
3885
3886      return res;
3887   }
3888
3889   /* setup regs */
3890   oldused = c->used;
3891   tmpa    = a->dp;
3892   tmpc    = c->dp;
3893
3894   /* if a <= b simply fix the single digit */
3895   if ((a->used == 1 && a->dp[0] <= b) || a->used == 0) {
3896      if (a->used == 1) {
3897         *tmpc++ = b - *tmpa;
3898      } else {
3899         *tmpc++ = b;
3900      }
3901      ix      = 1;
3902
3903      /* negative/1digit */
3904      c->sign = MP_NEG;
3905      c->used = 1;
3906   } else {
3907      /* positive/size */
3908      c->sign = MP_ZPOS;
3909      c->used = a->used;
3910
3911      /* subtract first digit */
3912      *tmpc    = *tmpa++ - b;
3913      mu       = *tmpc >> (sizeof(mp_digit) * CHAR_BIT - 1);
3914      *tmpc++ &= MP_MASK;
3915
3916      /* handle rest of the digits */
3917      for (ix = 1; ix < a->used; ix++) {
3918         *tmpc    = *tmpa++ - mu;
3919         mu       = *tmpc >> (sizeof(mp_digit) * CHAR_BIT - 1);
3920         *tmpc++ &= MP_MASK;
3921      }
3922   }
3923
3924   /* zero excess digits */
3925   while (ix++ < oldused) {
3926      *tmpc++ = 0;
3927   }
3928   mp_clamp(c);
3929   return MP_OKAY;
3930 }
3931
3932 #endif /* defined(HAVE_ECC) || !defined(NO_PWDBASED) */
3933
3934
3935 #ifdef CYASSL_KEY_GEN
3936
3937 int mp_cnt_lsb(mp_int *a);
3938
3939 static int s_is_power_of_two(mp_digit b, int *p)
3940 {
3941    int x;
3942
3943    /* fast return if no power of two */
3944    if ((b==0) || (b & (b-1))) {
3945       return 0;
3946    }
3947
3948    for (x = 0; x < DIGIT_BIT; x++) {
3949       if (b == (((mp_digit)1)<<x)) {
3950          *p = x;
3951          return 1;
3952       }
3953    }
3954    return 0;
3955 }
3956
3957 /* single digit division (based on routine from MPI) */
3958 static int mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d)
3959 {
3960   mp_int  q;
3961   mp_word w;
3962   mp_digit t;
3963   int     res, ix;
3964
3965   /* cannot divide by zero */
3966   if (b == 0) {
3967      return MP_VAL;
3968   }
3969
3970   /* quick outs */
3971   if (b == 1 || mp_iszero(a) == 1) {
3972      if (d != NULL) {
3973         *d = 0;
3974      }
3975      if (c != NULL) {
3976         return mp_copy(a, c);
3977      }
3978      return MP_OKAY;
3979   }
3980
3981   /* power of two ? */
3982   if (s_is_power_of_two(b, &ix) == 1) {
3983      if (d != NULL) {
3984         *d = a->dp[0] & ((((mp_digit)1)<<ix) - 1);
3985      }
3986      if (c != NULL) {
3987         return mp_div_2d(a, ix, c, NULL);
3988      }
3989      return MP_OKAY;
3990   }
3991
3992 #ifdef BN_MP_DIV_3_C
3993   /* three? */
3994   if (b == 3) {
3995      return mp_div_3(a, c, d);
3996   }
3997 #endif
3998
3999   /* no easy answer [c'est la vie].  Just division */
4000   if ((res = mp_init_size(&q, a->used)) != MP_OKAY) {
4001      return res;
4002   }
4003
4004   q.used = a->used;
4005   q.sign = a->sign;
4006   w = 0;
4007   for (ix = a->used - 1; ix >= 0; ix--) {
4008      w = (w << ((mp_word)DIGIT_BIT)) | ((mp_word)a->dp[ix]);
4009
4010      if (w >= b) {
4011         t = (mp_digit)(w / b);
4012         w -= ((mp_word)t) * ((mp_word)b);
4013       } else {
4014         t = 0;
4015       }
4016       q.dp[ix] = (mp_digit)t;
4017   }
4018
4019   if (d != NULL) {
4020      *d = (mp_digit)w;
4021   }
4022
4023   if (c != NULL) {
4024      mp_clamp(&q);
4025      mp_exch(&q, c);
4026   }
4027   mp_clear(&q);
4028
4029   return res;
4030 }
4031
4032
4033 static int mp_mod_d (mp_int * a, mp_digit b, mp_digit * c)
4034 {
4035   return mp_div_d(a, b, NULL, c);
4036 }
4037
4038
4039 const mp_digit ltm_prime_tab[] = {
4040   0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
4041   0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
4042   0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
4043   0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F,
4044 #ifndef MP_8BIT
4045   0x0083,
4046   0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
4047   0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
4048   0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
4049   0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
4050
4051   0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
4052   0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
4053   0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
4054   0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
4055   0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
4056   0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
4057   0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
4058   0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
4059
4060   0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
4061   0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
4062   0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
4063   0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
4064   0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
4065   0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
4066   0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
4067   0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
4068
4069   0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
4070   0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
4071   0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
4072   0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
4073   0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
4074   0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
4075   0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
4076   0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
4077 #endif
4078 };
4079
4080
4081 /* Miller-Rabin test of "a" to the base of "b" as described in
4082  * HAC pp. 139 Algorithm 4.24
4083  *
4084  * Sets result to 0 if definitely composite or 1 if probably prime.
4085  * Randomly the chance of error is no more than 1/4 and often
4086  * very much lower.
4087  */
4088 static int mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
4089 {
4090   mp_int  n1, y, r;
4091   int     s, j, err;
4092
4093   /* default */
4094   *result = MP_NO;
4095
4096   /* ensure b > 1 */
4097   if (mp_cmp_d(b, 1) != MP_GT) {
4098      return MP_VAL;
4099   }
4100
4101   /* get n1 = a - 1 */
4102   if ((err = mp_init_copy (&n1, a)) != MP_OKAY) {
4103     return err;
4104   }
4105   if ((err = mp_sub_d (&n1, 1, &n1)) != MP_OKAY) {
4106     goto LBL_N1;
4107   }
4108
4109   /* set 2**s * r = n1 */
4110   if ((err = mp_init_copy (&r, &n1)) != MP_OKAY) {
4111     goto LBL_N1;
4112   }
4113
4114   /* count the number of least significant bits
4115    * which are zero
4116    */
4117   s = mp_cnt_lsb(&r);
4118
4119   /* now divide n - 1 by 2**s */
4120   if ((err = mp_div_2d (&r, s, &r, NULL)) != MP_OKAY) {
4121     goto LBL_R;
4122   }
4123
4124   /* compute y = b**r mod a */
4125   if ((err = mp_init (&y)) != MP_OKAY) {
4126     goto LBL_R;
4127   }
4128   if ((err = mp_exptmod (b, &r, a, &y)) != MP_OKAY) {
4129     goto LBL_Y;
4130   }
4131
4132   /* if y != 1 and y != n1 do */
4133   if (mp_cmp_d (&y, 1) != MP_EQ && mp_cmp (&y, &n1) != MP_EQ) {
4134     j = 1;
4135     /* while j <= s-1 and y != n1 */
4136     while ((j <= (s - 1)) && mp_cmp (&y, &n1) != MP_EQ) {
4137       if ((err = mp_sqrmod (&y, a, &y)) != MP_OKAY) {
4138          goto LBL_Y;
4139       }
4140
4141       /* if y == 1 then composite */
4142       if (mp_cmp_d (&y, 1) == MP_EQ) {
4143          goto LBL_Y;
4144       }
4145
4146       ++j;
4147     }
4148
4149     /* if y != n1 then composite */
4150     if (mp_cmp (&y, &n1) != MP_EQ) {
4151       goto LBL_Y;
4152     }
4153   }
4154
4155   /* probably prime now */
4156   *result = MP_YES;
4157 LBL_Y:mp_clear (&y);
4158 LBL_R:mp_clear (&r);
4159 LBL_N1:mp_clear (&n1);
4160   return err;
4161 }
4162
4163
4164 /* determines if an integers is divisible by one
4165  * of the first PRIME_SIZE primes or not
4166  *
4167  * sets result to 0 if not, 1 if yes
4168  */
4169 static int mp_prime_is_divisible (mp_int * a, int *result)
4170 {
4171   int     err, ix;
4172   mp_digit res;
4173
4174   /* default to not */
4175   *result = MP_NO;
4176
4177   for (ix = 0; ix < PRIME_SIZE; ix++) {
4178     /* what is a mod LBL_prime_tab[ix] */
4179     if ((err = mp_mod_d (a, ltm_prime_tab[ix], &res)) != MP_OKAY) {
4180       return err;
4181     }
4182
4183     /* is the residue zero? */
4184     if (res == 0) {
4185       *result = MP_YES;
4186       return MP_OKAY;
4187     }
4188   }
4189
4190   return MP_OKAY;
4191 }
4192
4193
4194 /*
4195  * Sets result to 1 if probably prime, 0 otherwise
4196  */
4197 int mp_prime_is_prime (mp_int * a, int t, int *result)
4198 {
4199   mp_int  b;
4200   int     ix, err, res;
4201
4202   /* default to no */
4203   *result = MP_NO;
4204
4205   /* valid value of t? */
4206   if (t <= 0 || t > PRIME_SIZE) {
4207     return MP_VAL;
4208   }
4209
4210   /* is the input equal to one of the primes in the table? */
4211   for (ix = 0; ix < PRIME_SIZE; ix++) {
4212       if (mp_cmp_d(a, ltm_prime_tab[ix]) == MP_EQ) {
4213          *result = 1;
4214          return MP_OKAY;
4215       }
4216   }
4217
4218   /* first perform trial division */
4219   if ((err = mp_prime_is_divisible (a, &res)) != MP_OKAY) {
4220     return err;
4221   }
4222
4223   /* return if it was trivially divisible */
4224   if (res == MP_YES) {
4225     return MP_OKAY;
4226   }
4227
4228   /* now perform the miller-rabin rounds */
4229   if ((err = mp_init (&b)) != MP_OKAY) {
4230     return err;
4231   }
4232
4233   for (ix = 0; ix < t; ix++) {
4234     /* set the prime */
4235     mp_set (&b, ltm_prime_tab[ix]);
4236
4237     if ((err = mp_prime_miller_rabin (a, &b, &res)) != MP_OKAY) {
4238       goto LBL_B;
4239     }
4240
4241     if (res == MP_NO) {
4242       goto LBL_B;
4243     }
4244   }
4245
4246   /* passed the test */
4247   *result = MP_YES;
4248 LBL_B:mp_clear (&b);
4249   return err;
4250 }
4251
4252
4253 /* computes least common multiple as |a*b|/(a, b) */
4254 int mp_lcm (mp_int * a, mp_int * b, mp_int * c)
4255 {
4256   int     res;
4257   mp_int  t1, t2;
4258
4259
4260   if ((res = mp_init_multi (&t1, &t2, NULL, NULL, NULL, NULL)) != MP_OKAY) {
4261     return res;
4262   }
4263
4264   /* t1 = get the GCD of the two inputs */
4265   if ((res = mp_gcd (a, b, &t1)) != MP_OKAY) {
4266     goto LBL_T;
4267   }
4268
4269   /* divide the smallest by the GCD */
4270   if (mp_cmp_mag(a, b) == MP_LT) {
4271      /* store quotient in t2 such that t2 * b is the LCM */
4272      if ((res = mp_div(a, &t1, &t2, NULL)) != MP_OKAY) {
4273         goto LBL_T;
4274      }
4275      res = mp_mul(b, &t2, c);
4276   } else {
4277      /* store quotient in t2 such that t2 * a is the LCM */
4278      if ((res = mp_div(b, &t1, &t2, NULL)) != MP_OKAY) {
4279         goto LBL_T;
4280      }
4281      res = mp_mul(a, &t2, c);
4282   }
4283
4284   /* fix the sign to positive */
4285   c->sign = MP_ZPOS;
4286
4287 LBL_T:
4288   mp_clear(&t1);
4289   mp_clear(&t2);
4290   return res;
4291 }
4292
4293
4294 static const int lnz[16] = {
4295    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
4296 };
4297
4298 /* Counts the number of lsbs which are zero before the first zero bit */
4299 int mp_cnt_lsb(mp_int *a)
4300 {
4301     int x;
4302     mp_digit q, qq;
4303
4304     /* easy out */
4305     if (mp_iszero(a) == 1) {
4306         return 0;
4307     }
4308
4309     /* scan lower digits until non-zero */
4310     for (x = 0; x < a->used && a->dp[x] == 0; x++);
4311     q = a->dp[x];
4312     x *= DIGIT_BIT;
4313
4314     /* now scan this digit until a 1 is found */
4315     if ((q & 1) == 0) {
4316         do {
4317             qq  = q & 15;
4318             x  += lnz[qq];
4319             q >>= 4;
4320         } while (qq == 0);
4321     }
4322     return x;
4323 }
4324
4325
4326 /* Greatest Common Divisor using the binary method */
4327 int mp_gcd (mp_int * a, mp_int * b, mp_int * c)
4328 {
4329     mp_int  u, v;
4330     int     k, u_lsb, v_lsb, res;
4331
4332     /* either zero than gcd is the largest */
4333     if (mp_iszero (a) == MP_YES) {
4334         return mp_abs (b, c);
4335     }
4336     if (mp_iszero (b) == MP_YES) {
4337         return mp_abs (a, c);
4338     }
4339
4340     /* get copies of a and b we can modify */
4341     if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
4342         return res;
4343     }
4344
4345     if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
4346         goto LBL_U;
4347     }
4348
4349     /* must be positive for the remainder of the algorithm */
4350     u.sign = v.sign = MP_ZPOS;
4351
4352     /* B1.  Find the common power of two for u and v */
4353     u_lsb = mp_cnt_lsb(&u);
4354     v_lsb = mp_cnt_lsb(&v);
4355     k     = MIN(u_lsb, v_lsb);
4356
4357     if (k > 0) {
4358         /* divide the power of two out */
4359         if ((res = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
4360             goto LBL_V;
4361         }
4362
4363         if ((res = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
4364             goto LBL_V;
4365         }
4366     }
4367
4368     /* divide any remaining factors of two out */
4369     if (u_lsb != k) {
4370         if ((res = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
4371             goto LBL_V;
4372         }
4373     }
4374
4375     if (v_lsb != k) {
4376         if ((res = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
4377             goto LBL_V;
4378         }
4379     }
4380
4381     while (mp_iszero(&v) == 0) {
4382         /* make sure v is the largest */
4383         if (mp_cmp_mag(&u, &v) == MP_GT) {
4384             /* swap u and v to make sure v is >= u */
4385             mp_exch(&u, &v);
4386         }
4387
4388         /* subtract smallest from largest */
4389         if ((res = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
4390             goto LBL_V;
4391         }
4392
4393         /* Divide out all factors of two */
4394         if ((res = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
4395             goto LBL_V;
4396         }
4397     }
4398
4399     /* multiply by 2**k which we divided out at the beginning */
4400     if ((res = mp_mul_2d (&u, k, c)) != MP_OKAY) {
4401         goto LBL_V;
4402     }
4403     c->sign = MP_ZPOS;
4404     res = MP_OKAY;
4405 LBL_V:mp_clear (&u);
4406 LBL_U:mp_clear (&v);
4407     return res;
4408 }
4409
4410
4411
4412 #endif /* CYASSL_KEY_GEN */
4413
4414
4415 #ifdef HAVE_ECC
4416
4417 /* chars used in radix conversions */
4418 const char *mp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
4419
4420 /* read a string [ASCII] in a given radix */
4421 int mp_read_radix (mp_int * a, const char *str, int radix)
4422 {
4423   int     y, res, neg;
4424   char    ch;
4425
4426   /* zero the digit bignum */
4427   mp_zero(a);
4428
4429   /* make sure the radix is ok */
4430   if (radix < 2 || radix > 64) {
4431     return MP_VAL;
4432   }
4433
4434   /* if the leading digit is a
4435    * minus set the sign to negative.
4436    */
4437   if (*str == '-') {
4438     ++str;
4439     neg = MP_NEG;
4440   } else {
4441     neg = MP_ZPOS;
4442   }
4443
4444   /* set the integer to the default of zero */
4445   mp_zero (a);
4446
4447   /* process each digit of the string */
4448   while (*str) {
4449     /* if the radix < 36 the conversion is case insensitive
4450      * this allows numbers like 1AB and 1ab to represent the same  value
4451      * [e.g. in hex]
4452      */
4453     ch = (char) ((radix < 36) ? XTOUPPER(*str) : *str);
4454     for (y = 0; y < 64; y++) {
4455       if (ch == mp_s_rmap[y]) {
4456          break;
4457       }
4458     }
4459
4460     /* if the char was found in the map
4461      * and is less than the given radix add it
4462      * to the number, otherwise exit the loop.
4463      */
4464     if (y < radix) {
4465       if ((res = mp_mul_d (a, (mp_digit) radix, a)) != MP_OKAY) {
4466          return res;
4467       }
4468       if ((res = mp_add_d (a, (mp_digit) y, a)) != MP_OKAY) {
4469          return res;
4470       }
4471     } else {
4472       break;
4473     }
4474     ++str;
4475   }
4476
4477   /* set the sign only if a != 0 */
4478   if (mp_iszero(a) != 1) {
4479      a->sign = neg;
4480   }
4481   return MP_OKAY;
4482 }
4483
4484 #endif /* HAVE_ECC */
4485
4486 #endif /* USE_FAST_MATH */
4487
4488 #endif /* NO_BIG_INT */