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

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