src/tls/libtommath.c

   1 /*
   2  * Minimal code for RSA support from LibTomMath 0.3.9
   3  * http://math.libtomcrypt.com/
   4  * http://math.libtomcrypt.com/files/ltm-0.39.tar.bz2
   5  * This library was released in public domain by Tom St Denis.
   6  *
   7  * The combination in this file may not use all of the optimized algorithms
   8  * from LibTomMath and may be considerable slower than the LibTomMath with its
   9  * default settings. The main purpose of having this version here is to make it
  10  * easier to build bignum.c wrapper without having to install and build an
  11  * external library.
  12  *
  13  * If CONFIG_INTERNAL_LIBTOMMATH is defined, bignum.c includes this
  14  * libtommath.c file instead of using the external LibTomMath library.
  15  */
  16
  17 #ifndef CHAR_BIT
  18 #define CHAR_BIT 8
  19 #endif
  20
  21 #define BN_MP_INVMOD_C
  22 #define BN_S_MP_EXPTMOD_C /* Note: #undef in tommath_superclass.h; this would
  23                            * require BN_MP_EXPTMOD_FAST_C instead */
  24 #define BN_S_MP_MUL_DIGS_C
  25 #define BN_MP_INVMOD_SLOW_C
  26 #define BN_S_MP_SQR_C
  27 #define BN_S_MP_MUL_HIGH_DIGS_C /* Note: #undef in tommath_superclass.h; this
  28                                  * would require other than mp_reduce */
  29
  30 #ifdef LTM_FAST_EXPTMOD
  31 /* Include faster exptmod (Montgomery) at the cost of about 2.5 kB in code */
  32 #define BN_MP_EXPTMOD_FAST_C
  33 #define BN_MP_MONTGOMERY_SETUP_C
  34 #define BN_FAST_MP_MONTGOMERY_REDUCE_C
  35 #define BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
  36 #define BN_MP_MUL_2_C
  37 #endif /* LTM_FAST_EXPTMOD */
  38
  39 #ifdef LTM_FAST_SQR
  40 /* Include faster sqr at the cost of about 0.5 kB in code */
  41 #define BN_FAST_S_MP_SQR_C
  42 #endif /* LTM_FAST_SQR */
  43
  44 /* Current uses do not require support for negative exponent in exptmod, so we
  45  * can save about 1.5 kB in leaving out invmod. */
  46 #define LTM_NO_NEG_EXP
  47
  48 /* from tommath.h */
  49
  50 #ifndef MIN
  51    #define MIN(x,y) ((x)<(y)?(x):(y))
  52 #endif
  53
  54 #ifndef MAX
  55    #define MAX(x,y) ((x)>(y)?(x):(y))
  56 #endif
  57
  58 #define  OPT_CAST(x)
  59
  60 typedef unsigned long mp_digit;
  61 typedef u64 mp_word;
  62
  63 #define DIGIT_BIT          28
  64 #define MP_28BIT
  65
  66
  67 #define XMALLOC  os_malloc
  68 #define XFREE    os_free
  69 #define XREALLOC os_realloc
  70
  71
  72 #define MP_MASK          ((((mp_digit)1)<<((mp_digit)DIGIT_BIT))-((mp_digit)1))
  73
  74 #define MP_LT        -1   /* less than */
  75 #define MP_EQ         0   /* equal to */
  76 #define MP_GT         1   /* greater than */
  77
  78 #define MP_ZPOS       0   /* positive integer */
  79 #define MP_NEG        1   /* negative */
  80
  81 #define MP_OKAY       0   /* ok result */
  82 #define MP_MEM        -2  /* out of mem */
  83 #define MP_VAL        -3  /* invalid input */
  84
  85 #define MP_YES        1   /* yes response */
  86 #define MP_NO         0   /* no response */
  87
  88 typedef int           mp_err;
  89
  90 /* define this to use lower memory usage routines (exptmods mostly) */
  91 #define MP_LOW_MEM
  92
  93 /* default precision */
  94 #ifndef MP_PREC
  95    #ifndef MP_LOW_MEM
  96       #define MP_PREC                 32     /* default digits of precision */
  97    #else
  98       #define MP_PREC                 8      /* default digits of precision */
  99    #endif
 100 #endif
 101
 102 /* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */
 103 #define MP_WARRAY               (1 << (sizeof(mp_word) * CHAR_BIT - 2 * DIGIT_BIT + 1))
 104
 105 /* the infamous mp_int structure */
 106 typedef struct  {
 107     int used, alloc, sign;
 108     mp_digit *dp;
 109 } mp_int;
 110
 111
 112 /* ---> Basic Manipulations <--- */
 113 #define mp_iszero(a) (((a)->used == 0) ? MP_YES : MP_NO)
 114 #define mp_iseven(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 0)) ? MP_YES : MP_NO)
 115 #define mp_isodd(a)  (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? MP_YES : MP_NO)
 116
 117
 118 /* prototypes for copied functions */
 119 #define s_mp_mul(a, b, c) s_mp_mul_digs(a, b, c, (a)->used + (b)->used + 1)
 120 static int s_mp_exptmod(mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode);
 121 static int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs);
 122 static int s_mp_sqr(mp_int * a, mp_int * b);
 123 static int s_mp_mul_high_digs(mp_int * a, mp_int * b, mp_int * c, int digs);
 124
 125 static int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs);
 126
 127 static int mp_init_multi(mp_int *mp, ...);
 128 static void mp_clear_multi(mp_int *mp, ...);
 129 static int mp_lshd(mp_int * a, int b);
 130 static void mp_set(mp_int * a, mp_digit b);
 131 static void mp_clamp(mp_int * a);
 132 static void mp_exch(mp_int * a, mp_int * b);
 133 static void mp_rshd(mp_int * a, int b);
 134 static void mp_zero(mp_int * a);
 135 static int mp_mod_2d(mp_int * a, int b, mp_int * c);
 136 static int mp_div_2d(mp_int * a, int b, mp_int * c, mp_int * d);
 137 static int mp_init_copy(mp_int * a, mp_int * b);
 138 static int mp_mul_2d(mp_int * a, int b, mp_int * c);
 139 #ifndef LTM_NO_NEG_EXP
 140 static int mp_div_2(mp_int * a, mp_int * b);
 141 static int mp_invmod(mp_int * a, mp_int * b, mp_int * c);
 142 static int mp_invmod_slow(mp_int * a, mp_int * b, mp_int * c);
 143 #endif /* LTM_NO_NEG_EXP */
 144 static int mp_copy(mp_int * a, mp_int * b);
 145 static int mp_count_bits(mp_int * a);
 146 static int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d);
 147 static int mp_mod(mp_int * a, mp_int * b, mp_int * c);
 148 static int mp_grow(mp_int * a, int size);
 149 static int mp_cmp_mag(mp_int * a, mp_int * b);
 150 static int mp_abs(mp_int * a, mp_int * b);
 151 static int mp_sqr(mp_int * a, mp_int * b);
 152 static int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d);
 153 static int mp_reduce_2k_setup_l(mp_int *a, mp_int *d);
 154 static int mp_2expt(mp_int * a, int b);
 155 static int mp_reduce_setup(mp_int * a, mp_int * b);
 156 static int mp_reduce(mp_int * x, mp_int * m, mp_int * mu);
 157 static int mp_init_size(mp_int * a, int size);
 158 #ifdef BN_MP_EXPTMOD_FAST_C
 159 static int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode);
 160 #endif /* BN_MP_EXPTMOD_FAST_C */
 161 #ifdef BN_FAST_S_MP_SQR_C
 162 static int fast_s_mp_sqr (mp_int * a, mp_int * b);
 163 #endif /* BN_FAST_S_MP_SQR_C */
 164
 165
 166
 167 /* functions from bn_<func name>.c */
 168
 169
 170 /* reverse an array, used for radix code */
 171 static void bn_reverse (unsigned char *s, int len)
 172 {
 173   int     ix, iy;
 174   unsigned char t;
 175
 176   ix = 0;
 177   iy = len - 1;
 178   while (ix < iy) {
 179     t     = s[ix];
 180     s[ix] = s[iy];
 181     s[iy] = t;
 182     ++ix;
 183     --iy;
 184   }
 185 }
 186
 187
 188 /* low level addition, based on HAC pp.594, Algorithm 14.7 */
 189 static int s_mp_add (mp_int * a, mp_int * b, mp_int * c)
 190 {
 191   mp_int *x;
 192   int     olduse, res, min, max;
 193
 194   /* find sizes, we let |a| <= |b| which means we have to sort
 195    * them.  "x" will point to the input with the most digits
 196    */
 197   if (a->used > b->used) {
 198     min = b->used;
 199     max = a->used;
 200     x = a;
 201   } else {
 202     min = a->used;
 203     max = b->used;
 204     x = b;
 205   }
 206
 207   /* init result */
 208   if (c->alloc < max + 1) {
 209     if ((res = mp_grow (c, max + 1)) != MP_OKAY) {
 210       return res;
 211     }
 212   }
 213
 214   /* get old used digit count and set new one */
 215   olduse = c->used;
 216   c->used = max + 1;
 217
 218   {
 219     register mp_digit u, *tmpa, *tmpb, *tmpc;
 220     register int i;
 221
 222     /* alias for digit pointers */
 223
 224     /* first input */
 225     tmpa = a->dp;
 226
 227     /* second input */
 228     tmpb = b->dp;
 229
 230     /* destination */
 231     tmpc = c->dp;
 232
 233     /* zero the carry */
 234     u = 0;
 235     for (i = 0; i < min; i++) {
 236       /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
 237       *tmpc = *tmpa++ + *tmpb++ + u;
 238
 239       /* U = carry bit of T[i] */
 240       u = *tmpc >> ((mp_digit)DIGIT_BIT);
 241
 242       /* take away carry bit from T[i] */
 243       *tmpc++ &= MP_MASK;
 244     }
 245
 246     /* now copy higher words if any, that is in A+B
 247      * if A or B has more digits add those in
 248      */
 249     if (min != max) {
 250       for (; i < max; i++) {
 251         /* T[i] = X[i] + U */
 252         *tmpc = x->dp[i] + u;
 253
 254         /* U = carry bit of T[i] */
 255         u = *tmpc >> ((mp_digit)DIGIT_BIT);
 256
 257         /* take away carry bit from T[i] */
 258         *tmpc++ &= MP_MASK;
 259       }
 260     }
 261
 262     /* add carry */
 263     *tmpc++ = u;
 264
 265     /* clear digits above oldused */
 266     for (i = c->used; i < olduse; i++) {
 267       *tmpc++ = 0;
 268     }
 269   }
 270
 271   mp_clamp (c);
 272   return MP_OKAY;
 273 }
 274
 275
 276 /* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
 277 static int s_mp_sub (mp_int * a, mp_int * b, mp_int * c)
 278 {
 279   int     olduse, res, min, max;
 280
 281   /* find sizes */
 282   min = b->used;
 283   max = a->used;
 284
 285   /* init result */
 286   if (c->alloc < max) {
 287     if ((res = mp_grow (c, max)) != MP_OKAY) {
 288       return res;
 289     }
 290   }
 291   olduse = c->used;
 292   c->used = max;
 293
 294   {
 295     register mp_digit u, *tmpa, *tmpb, *tmpc;
 296     register int i;
 297
 298     /* alias for digit pointers */
 299     tmpa = a->dp;
 300     tmpb = b->dp;
 301     tmpc = c->dp;
 302
 303     /* set carry to zero */
 304     u = 0;
 305     for (i = 0; i < min; i++) {
 306       /* T[i] = A[i] - B[i] - U */
 307       *tmpc = *tmpa++ - *tmpb++ - u;
 308
 309       /* U = carry bit of T[i]
 310        * Note this saves performing an AND operation since
 311        * if a carry does occur it will propagate all the way to the
 312        * MSB.  As a result a single shift is enough to get the carry
 313        */
 314       u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
 315
 316       /* Clear carry from T[i] */
 317       *tmpc++ &= MP_MASK;
 318     }
 319
 320     /* now copy higher words if any, e.g. if A has more digits than B  */
 321     for (; i < max; i++) {
 322       /* T[i] = A[i] - U */
 323       *tmpc = *tmpa++ - u;
 324
 325       /* U = carry bit of T[i] */
 326       u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
 327
 328       /* Clear carry from T[i] */
 329       *tmpc++ &= MP_MASK;
 330     }
 331
 332     /* clear digits above used (since we may not have grown result above) */
 333     for (i = c->used; i < olduse; i++) {
 334       *tmpc++ = 0;
 335     }
 336   }
 337
 338   mp_clamp (c);
 339   return MP_OKAY;
 340 }
 341
 342
 343 /* init a new mp_int */
 344 static int mp_init (mp_int * a)
 345 {
 346   int i;
 347
 348   /* allocate memory required and clear it */
 349   a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * MP_PREC);
 350   if (a->dp == NULL) {
 351     return MP_MEM;
 352   }
 353
 354   /* set the digits to zero */
 355   for (i = 0; i < MP_PREC; i++) {
 356       a->dp[i] = 0;
 357   }
 358
 359   /* set the used to zero, allocated digits to the default precision
 360    * and sign to positive */
 361   a->used  = 0;
 362   a->alloc = MP_PREC;
 363   a->sign  = MP_ZPOS;
 364
 365   return MP_OKAY;
 366 }
 367
 368
 369 /* clear one (frees)  */
 370 static void mp_clear (mp_int * a)
 371 {
 372   int i;
 373
 374   /* only do anything if a hasn't been freed previously */
 375   if (a->dp != NULL) {
 376     /* first zero the digits */
 377     for (i = 0; i < a->used; i++) {
 378         a->dp[i] = 0;
 379     }
 380
 381     /* free ram */
 382     XFREE(a->dp);
 383
 384     /* reset members to make debugging easier */
 385     a->dp    = NULL;
 386     a->alloc = a->used = 0;
 387     a->sign  = MP_ZPOS;
 388   }
 389 }
 390
 391
 392 /* high level addition (handles signs) */
 393 static int mp_add (mp_int * a, mp_int * b, mp_int * c)
 394 {
 395   int     sa, sb, res;
 396
 397   /* get sign of both inputs */
 398   sa = a->sign;
 399   sb = b->sign;
 400
 401   /* handle two cases, not four */
 402   if (sa == sb) {
 403     /* both positive or both negative */
 404     /* add their magnitudes, copy the sign */
 405     c->sign = sa;
 406     res = s_mp_add (a, b, c);
 407   } else {
 408     /* one positive, the other negative */
 409     /* subtract the one with the greater magnitude from */
 410     /* the one of the lesser magnitude.  The result gets */
 411     /* the sign of the one with the greater magnitude. */
 412     if (mp_cmp_mag (a, b) == MP_LT) {
 413       c->sign = sb;
 414       res = s_mp_sub (b, a, c);
 415     } else {
 416       c->sign = sa;
 417       res = s_mp_sub (a, b, c);
 418     }
 419   }
 420   return res;
 421 }
 422
 423
 424 /* high level subtraction (handles signs) */
 425 static int mp_sub (mp_int * a, mp_int * b, mp_int * c)
 426 {
 427   int     sa, sb, res;
 428
 429   sa = a->sign;
 430   sb = b->sign;
 431
 432   if (sa != sb) {
 433     /* subtract a negative from a positive, OR */
 434     /* subtract a positive from a negative. */
 435     /* In either case, ADD their magnitudes, */
 436     /* and use the sign of the first number. */
 437     c->sign = sa;
 438     res = s_mp_add (a, b, c);
 439   } else {
 440     /* subtract a positive from a positive, OR */
 441     /* subtract a negative from a negative. */
 442     /* First, take the difference between their */
 443     /* magnitudes, then... */
 444     if (mp_cmp_mag (a, b) != MP_LT) {
 445       /* Copy the sign from the first */
 446       c->sign = sa;
 447       /* The first has a larger or equal magnitude */
 448       res = s_mp_sub (a, b, c);
 449     } else {
 450       /* The result has the *opposite* sign from */
 451       /* the first number. */
 452       c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
 453       /* The second has a larger magnitude */
 454       res = s_mp_sub (b, a, c);
 455     }
 456   }
 457   return res;
 458 }
 459
 460
 461 /* high level multiplication (handles sign) */
 462 static int mp_mul (mp_int * a, mp_int * b, mp_int * c)
 463 {
 464   int     res, neg;
 465   neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
 466
 467   /* use Toom-Cook? */
 468 #ifdef BN_MP_TOOM_MUL_C
 469   if (MIN (a->used, b->used) >= TOOM_MUL_CUTOFF) {
 470     res = mp_toom_mul(a, b, c);
 471   } else
 472 #endif
 473 #ifdef BN_MP_KARATSUBA_MUL_C
 474   /* use Karatsuba? */
 475   if (MIN (a->used, b->used) >= KARATSUBA_MUL_CUTOFF) {
 476     res = mp_karatsuba_mul (a, b, c);
 477   } else
 478 #endif
 479   {
 480     /* can we use the fast multiplier?
 481      *
 482      * The fast multiplier can be used if the output will
 483      * have less than MP_WARRAY digits and the number of
 484      * digits won't affect carry propagation
 485      */
 486 #ifdef BN_FAST_S_MP_MUL_DIGS_C
 487     int     digs = a->used + b->used + 1;
 488
 489     if ((digs < MP_WARRAY) &&
 490         MIN(a->used, b->used) <=
 491         (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
 492       res = fast_s_mp_mul_digs (a, b, c, digs);
 493     } else
 494 #endif
 495 #ifdef BN_S_MP_MUL_DIGS_C
 496       res = s_mp_mul (a, b, c); /* uses s_mp_mul_digs */
 497 #else
 498 #error mp_mul could fail
 499       res = MP_VAL;
 500 #endif
 501
 502   }
 503   c->sign = (c->used > 0) ? neg : MP_ZPOS;
 504   return res;
 505 }
 506
 507
 508 /* d = a * b (mod c) */
 509 static int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
 510 {
 511   int     res;
 512   mp_int  t;
 513
 514   if ((res = mp_init (&t)) != MP_OKAY) {
 515     return res;
 516   }
 517
 518   if ((res = mp_mul (a, b, &t)) != MP_OKAY) {
 519     mp_clear (&t);
 520     return res;
 521   }
 522   res = mp_mod (&t, c, d);
 523   mp_clear (&t);
 524   return res;
 525 }
 526
 527
 528 /* c = a mod b, 0 <= c < b */
 529 static int mp_mod (mp_int * a, mp_int * b, mp_int * c)
 530 {
 531   mp_int  t;
 532   int     res;
 533
 534   if ((res = mp_init (&t)) != MP_OKAY) {
 535     return res;
 536   }
 537
 538   if ((res = mp_div (a, b, NULL, &t)) != MP_OKAY) {
 539     mp_clear (&t);
 540     return res;
 541   }
 542
 543   if (t.sign != b->sign) {
 544     res = mp_add (b, &t, c);
 545   } else {
 546     res = MP_OKAY;
 547     mp_exch (&t, c);
 548   }
 549
 550   mp_clear (&t);
 551   return res;
 552 }
 553
 554
 555 /* this is a shell function that calls either the normal or Montgomery
 556  * exptmod functions.  Originally the call to the montgomery code was
 557  * embedded in the normal function but that wasted alot of stack space
 558  * for nothing (since 99% of the time the Montgomery code would be called)
 559  */
 560 static int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 561 {
 562   int dr;
 563
 564   /* modulus P must be positive */
 565   if (P->sign == MP_NEG) {
 566      return MP_VAL;
 567   }
 568
 569   /* if exponent X is negative we have to recurse */
 570   if (X->sign == MP_NEG) {
 571 #ifdef LTM_NO_NEG_EXP
 572         return MP_VAL;
 573 #else /* LTM_NO_NEG_EXP */
 574 #ifdef BN_MP_INVMOD_C
 575      mp_int tmpG, tmpX;
 576      int err;
 577
 578      /* first compute 1/G mod P */
 579      if ((err = mp_init(&tmpG)) != MP_OKAY) {
 580         return err;
 581      }
 582      if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {
 583         mp_clear(&tmpG);
 584         return err;
 585      }
 586
 587      /* now get |X| */
 588      if ((err = mp_init(&tmpX)) != MP_OKAY) {
 589         mp_clear(&tmpG);
 590         return err;
 591      }
 592      if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
 593         mp_clear_multi(&tmpG, &tmpX, NULL);
 594         return err;
 595      }
 596
 597      /* and now compute (1/G)**|X| instead of G**X [X < 0] */
 598      err = mp_exptmod(&tmpG, &tmpX, P, Y);
 599      mp_clear_multi(&tmpG, &tmpX, NULL);
 600      return err;
 601 #else
 602 #error mp_exptmod would always fail
 603      /* no invmod */
 604      return MP_VAL;
 605 #endif
 606 #endif /* LTM_NO_NEG_EXP */
 607   }
 608
 609 /* modified diminished radix reduction */
 610 #if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C) && defined(BN_S_MP_EXPTMOD_C)
 611   if (mp_reduce_is_2k_l(P) == MP_YES) {
 612      return s_mp_exptmod(G, X, P, Y, 1);
 613   }
 614 #endif
 615
 616 #ifdef BN_MP_DR_IS_MODULUS_C
 617   /* is it a DR modulus? */
 618   dr = mp_dr_is_modulus(P);
 619 #else
 620   /* default to no */
 621   dr = 0;
 622 #endif
 623
 624 #ifdef BN_MP_REDUCE_IS_2K_C
 625   /* if not, is it a unrestricted DR modulus? */
 626   if (dr == 0) {
 627      dr = mp_reduce_is_2k(P) << 1;
 628   }
 629 #endif
 630
 631   /* if the modulus is odd or dr != 0 use the montgomery method */
 632 #ifdef BN_MP_EXPTMOD_FAST_C
 633   if (mp_isodd (P) == 1 || dr !=  0) {
 634     return mp_exptmod_fast (G, X, P, Y, dr);
 635   } else {
 636 #endif
 637 #ifdef BN_S_MP_EXPTMOD_C
 638     /* otherwise use the generic Barrett reduction technique */
 639     return s_mp_exptmod (G, X, P, Y, 0);
 640 #else
 641 #error mp_exptmod could fail
 642     /* no exptmod for evens */
 643     return MP_VAL;
 644 #endif
 645 #ifdef BN_MP_EXPTMOD_FAST_C
 646   }
 647 #endif
 648 }
 649
 650
 651 /* compare two ints (signed)*/
 652 static int mp_cmp (mp_int * a, mp_int * b)
 653 {
 654   /* compare based on sign */
 655   if (a->sign != b->sign) {
 656      if (a->sign == MP_NEG) {
 657         return MP_LT;
 658      } else {
 659         return MP_GT;
 660      }
 661   }
 662
 663   /* compare digits */
 664   if (a->sign == MP_NEG) {
 665      /* if negative compare opposite direction */
 666      return mp_cmp_mag(b, a);
 667   } else {
 668      return mp_cmp_mag(a, b);
 669   }
 670 }
 671
 672
 673 /* compare a digit */
 674 static int mp_cmp_d(mp_int * a, mp_digit b)
 675 {
 676   /* compare based on sign */
 677   if (a->sign == MP_NEG) {
 678     return MP_LT;
 679   }
 680
 681   /* compare based on magnitude */
 682   if (a->used > 1) {
 683     return MP_GT;
 684   }
 685
 686   /* compare the only digit of a to b */
 687   if (a->dp[0] > b) {
 688     return MP_GT;
 689   } else if (a->dp[0] < b) {
 690     return MP_LT;
 691   } else {
 692     return MP_EQ;
 693   }
 694 }
 695
 696
 697 #ifndef LTM_NO_NEG_EXP
 698 /* hac 14.61, pp608 */
 699 static int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
 700 {
 701   /* b cannot be negative */
 702   if (b->sign == MP_NEG || mp_iszero(b) == 1) {
 703     return MP_VAL;
 704   }
 705
 706 #ifdef BN_FAST_MP_INVMOD_C
 707   /* if the modulus is odd we can use a faster routine instead */
 708   if (mp_isodd (b) == 1) {
 709     return fast_mp_invmod (a, b, c);
 710   }
 711 #endif
 712
 713 #ifdef BN_MP_INVMOD_SLOW_C
 714   return mp_invmod_slow(a, b, c);
 715 #endif
 716
 717 #ifndef BN_FAST_MP_INVMOD_C
 718 #ifndef BN_MP_INVMOD_SLOW_C
 719 #error mp_invmod would always fail
 720 #endif
 721 #endif
 722   return MP_VAL;
 723 }
 724 #endif /* LTM_NO_NEG_EXP */
 725
 726
 727 /* get the size for an unsigned equivalent */
 728 static int mp_unsigned_bin_size (mp_int * a)
 729 {
 730   int     size = mp_count_bits (a);
 731   return (size / 8 + ((size & 7) != 0 ? 1 : 0));
 732 }
 733
 734
 735 #ifndef LTM_NO_NEG_EXP
 736 /* hac 14.61, pp608 */
 737 static int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c)
 738 {
 739   mp_int  x, y, u, v, A, B, C, D;
 740   int     res;
 741
 742   /* b cannot be negative */
 743   if (b->sign == MP_NEG || mp_iszero(b) == 1) {
 744     return MP_VAL;
 745   }
 746
 747   /* init temps */
 748   if ((res = mp_init_multi(&x, &y, &u, &v,
 749                            &A, &B, &C, &D, NULL)) != MP_OKAY) {
 750      return res;
 751   }
 752
 753   /* x = a, y = b */
 754   if ((res = mp_mod(a, b, &x)) != MP_OKAY) {
 755       goto LBL_ERR;
 756   }
 757   if ((res = mp_copy (b, &y)) != MP_OKAY) {
 758     goto LBL_ERR;
 759   }
 760
 761   /* 2. [modified] if x,y are both even then return an error! */
 762   if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) {
 763     res = MP_VAL;
 764     goto LBL_ERR;
 765   }
 766
 767   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
 768   if ((res = mp_copy (&x, &u)) != MP_OKAY) {
 769     goto LBL_ERR;
 770   }
 771   if ((res = mp_copy (&y, &v)) != MP_OKAY) {
 772     goto LBL_ERR;
 773   }
 774   mp_set (&A, 1);
 775   mp_set (&D, 1);
 776
 777 top:
 778   /* 4.  while u is even do */
 779   while (mp_iseven (&u) == 1) {
 780     /* 4.1 u = u/2 */
 781     if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
 782       goto LBL_ERR;
 783     }
 784     /* 4.2 if A or B is odd then */
 785     if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) {
 786       /* A = (A+y)/2, B = (B-x)/2 */
 787       if ((res = mp_add (&A, &y, &A)) != MP_OKAY) {
 788          goto LBL_ERR;
 789       }
 790       if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
 791          goto LBL_ERR;
 792       }
 793     }
 794     /* A = A/2, B = B/2 */
 795     if ((res = mp_div_2 (&A, &A)) != MP_OKAY) {
 796       goto LBL_ERR;
 797     }
 798     if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
 799       goto LBL_ERR;
 800     }
 801   }
 802
 803   /* 5.  while v is even do */
 804   while (mp_iseven (&v) == 1) {
 805     /* 5.1 v = v/2 */
 806     if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
 807       goto LBL_ERR;
 808     }
 809     /* 5.2 if C or D is odd then */
 810     if (mp_isodd (&C) == 1 || mp_isodd (&D) == 1) {
 811       /* C = (C+y)/2, D = (D-x)/2 */
 812       if ((res = mp_add (&C, &y, &C)) != MP_OKAY) {
 813          goto LBL_ERR;
 814       }
 815       if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
 816          goto LBL_ERR;
 817       }
 818     }
 819     /* C = C/2, D = D/2 */
 820     if ((res = mp_div_2 (&C, &C)) != MP_OKAY) {
 821       goto LBL_ERR;
 822     }
 823     if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
 824       goto LBL_ERR;
 825     }
 826   }
 827
 828   /* 6.  if u >= v then */
 829   if (mp_cmp (&u, &v) != MP_LT) {
 830     /* u = u - v, A = A - C, B = B - D */
 831     if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
 832       goto LBL_ERR;
 833     }
 834
 835     if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) {
 836       goto LBL_ERR;
 837     }
 838
 839     if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
 840       goto LBL_ERR;
 841     }
 842   } else {
 843     /* v - v - u, C = C - A, D = D - B */
 844     if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
 845       goto LBL_ERR;
 846     }
 847
 848     if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) {
 849       goto LBL_ERR;
 850     }
 851
 852     if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
 853       goto LBL_ERR;
 854     }
 855   }
 856
 857   /* if not zero goto step 4 */
 858   if (mp_iszero (&u) == 0)
 859     goto top;
 860
 861   /* now a = C, b = D, gcd == g*v */
 862
 863   /* if v != 1 then there is no inverse */
 864   if (mp_cmp_d (&v, 1) != MP_EQ) {
 865     res = MP_VAL;
 866     goto LBL_ERR;
 867   }
 868
 869   /* if its too low */
 870   while (mp_cmp_d(&C, 0) == MP_LT) {
 871       if ((res = mp_add(&C, b, &C)) != MP_OKAY) {
 872          goto LBL_ERR;
 873       }
 874   }
 875
 876   /* too big */
 877   while (mp_cmp_mag(&C, b) != MP_LT) {
 878       if ((res = mp_sub(&C, b, &C)) != MP_OKAY) {
 879          goto LBL_ERR;
 880       }
 881   }
 882
 883   /* C is now the inverse */
 884   mp_exch (&C, c);
 885   res = MP_OKAY;
 886 LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
 887   return res;
 888 }
 889 #endif /* LTM_NO_NEG_EXP */
 890
 891
 892 /* compare maginitude of two ints (unsigned) */
 893 static int mp_cmp_mag (mp_int * a, mp_int * b)
 894 {
 895   int     n;
 896   mp_digit *tmpa, *tmpb;
 897
 898   /* compare based on # of non-zero digits */
 899   if (a->used > b->used) {
 900     return MP_GT;
 901   }
 902
 903   if (a->used < b->used) {
 904     return MP_LT;
 905   }
 906
 907   /* alias for a */
 908   tmpa = a->dp + (a->used - 1);
 909
 910   /* alias for b */
 911   tmpb = b->dp + (a->used - 1);
 912
 913   /* compare based on digits  */
 914   for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
 915     if (*tmpa > *tmpb) {
 916       return MP_GT;
 917     }
 918
 919     if (*tmpa < *tmpb) {
 920       return MP_LT;
 921     }
 922   }
 923   return MP_EQ;
 924 }
 925
 926
 927 /* reads a unsigned char array, assumes the msb is stored first [big endian] */
 928 static int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c)
 929 {
 930   int     res;
 931
 932   /* make sure there are at least two digits */
 933   if (a->alloc < 2) {
 934      if ((res = mp_grow(a, 2)) != MP_OKAY) {
 935         return res;
 936      }
 937   }
 938
 939   /* zero the int */
 940   mp_zero (a);
 941
 942   /* read the bytes in */
 943   while (c-- > 0) {
 944     if ((res = mp_mul_2d (a, 8, a)) != MP_OKAY) {
 945       return res;
 946     }
 947
 948 #ifndef MP_8BIT
 949       a->dp[0] |= *b++;
 950       a->used += 1;
 951 #else
 952       a->dp[0] = (*b & MP_MASK);
 953       a->dp[1] |= ((*b++ >> 7U) & 1);
 954       a->used += 2;
 955 #endif
 956   }
 957   mp_clamp (a);
 958   return MP_OKAY;
 959 }
 960
 961
 962 /* store in unsigned [big endian] format */
 963 static int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
 964 {
 965   int     x, res;
 966   mp_int  t;
 967
 968   if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
 969     return res;
 970   }
 971
 972   x = 0;
 973   while (mp_iszero (&t) == 0) {
 974 #ifndef MP_8BIT
 975       b[x++] = (unsigned char) (t.dp[0] & 255);
 976 #else
 977       b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7));
 978 #endif
 979     if ((res = mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) {
 980       mp_clear (&t);
 981       return res;
 982     }
 983   }
 984   bn_reverse (b, x);
 985   mp_clear (&t);
 986   return MP_OKAY;
 987 }
 988
 989
 990 /* shift right by a certain bit count (store quotient in c, optional remainder in d) */
 991 static int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
 992 {
 993   mp_digit D, r, rr;
 994   int     x, res;
 995   mp_int  t;
 996
 997
 998   /* if the shift count is <= 0 then we do no work */
 999   if (b <= 0) {
1000     res = mp_copy (a, c);
1001     if (d != NULL) {
1002       mp_zero (d);
1003     }
1004     return res;
1005   }
1006
1007   if ((res = mp_init (&t)) != MP_OKAY) {
1008     return res;
1009   }
1010
1011   /* get the remainder */
1012   if (d != NULL) {
1013     if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) {
1014       mp_clear (&t);
1015       return res;
1016     }
1017   }
1018
1019   /* copy */
1020   if ((res = mp_copy (a, c)) != MP_OKAY) {
1021     mp_clear (&t);
1022     return res;
1023   }
1024
1025   /* shift by as many digits in the bit count */
1026   if (b >= (int)DIGIT_BIT) {
1027     mp_rshd (c, b / DIGIT_BIT);
1028   }
1029
1030   /* shift any bit count < DIGIT_BIT */
1031   D = (mp_digit) (b % DIGIT_BIT);
1032   if (D != 0) {
1033     register mp_digit *tmpc, mask, shift;
1034
1035     /* mask */
1036     mask = (((mp_digit)1) << D) - 1;
1037
1038     /* shift for lsb */
1039     shift = DIGIT_BIT - D;
1040
1041     /* alias */
1042     tmpc = c->dp + (c->used - 1);
1043
1044     /* carry */
1045     r = 0;
1046     for (x = c->used - 1; x >= 0; x--) {
1047       /* get the lower  bits of this word in a temp */
1048       rr = *tmpc & mask;
1049
1050       /* shift the current word and mix in the carry bits from the previous word */
1051       *tmpc = (*tmpc >> D) | (r << shift);
1052       --tmpc;
1053
1054       /* set the carry to the carry bits of the current word found above */
1055       r = rr;
1056     }
1057   }
1058   mp_clamp (c);
1059   if (d != NULL) {
1060     mp_exch (&t, d);
1061   }
1062   mp_clear (&t);
1063   return MP_OKAY;
1064 }
1065
1066
1067 static int mp_init_copy (mp_int * a, mp_int * b)
1068 {
1069   int     res;
1070
1071   if ((res = mp_init (a)) != MP_OKAY) {
1072     return res;
1073   }
1074   return mp_copy (b, a);
1075 }
1076
1077
1078 /* set to zero */
1079 static void mp_zero (mp_int * a)
1080 {
1081   int       n;
1082   mp_digit *tmp;
1083
1084   a->sign = MP_ZPOS;
1085   a->used = 0;
1086
1087   tmp = a->dp;
1088   for (n = 0; n < a->alloc; n++) {
1089      *tmp++ = 0;
1090   }
1091 }
1092
1093
1094 /* copy, b = a */
1095 static int mp_copy (mp_int * a, mp_int * b)
1096 {
1097   int     res, n;
1098
1099   /* if dst == src do nothing */
1100   if (a == b) {
1101     return MP_OKAY;
1102   }
1103
1104   /* grow dest */
1105   if (b->alloc < a->used) {
1106      if ((res = mp_grow (b, a->used)) != MP_OKAY) {
1107         return res;
1108      }
1109   }
1110
1111   /* zero b and copy the parameters over */
1112   {
1113     register mp_digit *tmpa, *tmpb;
1114
1115     /* pointer aliases */
1116
1117     /* source */
1118     tmpa = a->dp;
1119
1120     /* destination */
1121     tmpb = b->dp;
1122
1123     /* copy all the digits */
1124     for (n = 0; n < a->used; n++) {
1125       *tmpb++ = *tmpa++;
1126     }
1127
1128     /* clear high digits */
1129     for (; n < b->used; n++) {
1130       *tmpb++ = 0;
1131     }
1132   }
1133
1134   /* copy used count and sign */
1135   b->used = a->used;
1136   b->sign = a->sign;
1137   return MP_OKAY;
1138 }
1139
1140
1141 /* shift right a certain amount of digits */
1142 static void mp_rshd (mp_int * a, int b)
1143 {
1144   int     x;
1145
1146   /* if b <= 0 then ignore it */
1147   if (b <= 0) {
1148     return;
1149   }
1150
1151   /* if b > used then simply zero it and return */
1152   if (a->used <= b) {
1153     mp_zero (a);
1154     return;
1155   }
1156
1157   {
1158     register mp_digit *bottom, *top;
1159
1160     /* shift the digits down */
1161
1162     /* bottom */
1163     bottom = a->dp;
1164
1165     /* top [offset into digits] */
1166     top = a->dp + b;
1167
1168     /* this is implemented as a sliding window where
1169      * the window is b-digits long and digits from
1170      * the top of the window are copied to the bottom
1171      *
1172      * e.g.
1173
1174      b-2 | b-1 | b0 | b1 | b2 | ... | bb |   ---->
1175                  /\                   |      ---->
1176                   \-------------------/      ---->
1177      */
1178     for (x = 0; x < (a->used - b); x++) {
1179       *bottom++ = *top++;
1180     }
1181
1182     /* zero the top digits */
1183     for (; x < a->used; x++) {
1184       *bottom++ = 0;
1185     }
1186   }
1187
1188   /* remove excess digits */
1189   a->used -= b;
1190 }
1191
1192
1193 /* swap the elements of two integers, for cases where you can't simply swap the
1194  * mp_int pointers around
1195  */
1196 static void mp_exch (mp_int * a, mp_int * b)
1197 {
1198   mp_int  t;
1199
1200   t  = *a;
1201   *a = *b;
1202   *b = t;
1203 }
1204
1205
1206 /* trim unused digits
1207  *
1208  * This is used to ensure that leading zero digits are
1209  * trimed and the leading "used" digit will be non-zero
1210  * Typically very fast.  Also fixes the sign if there
1211  * are no more leading digits
1212  */
1213 static void mp_clamp (mp_int * a)
1214 {
1215   /* decrease used while the most significant digit is
1216    * zero.
1217    */
1218   while (a->used > 0 && a->dp[a->used - 1] == 0) {
1219     --(a->used);
1220   }
1221
1222   /* reset the sign flag if used == 0 */
1223   if (a->used == 0) {
1224     a->sign = MP_ZPOS;
1225   }
1226 }
1227
1228
1229 /* grow as required */
1230 static int mp_grow (mp_int * a, int size)
1231 {
1232   int     i;
1233   mp_digit *tmp;
1234
1235   /* if the alloc size is smaller alloc more ram */
1236   if (a->alloc < size) {
1237     /* ensure there are always at least MP_PREC digits extra on top */
1238     size += (MP_PREC * 2) - (size % MP_PREC);
1239
1240     /* reallocate the array a->dp
1241      *
1242      * We store the return in a temporary variable
1243      * in case the operation failed we don't want
1244      * to overwrite the dp member of a.
1245      */
1246     tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * size);
1247     if (tmp == NULL) {
1248       /* reallocation failed but "a" is still valid [can be freed] */
1249       return MP_MEM;
1250     }
1251
1252     /* reallocation succeeded so set a->dp */
1253     a->dp = tmp;
1254
1255     /* zero excess digits */
1256     i        = a->alloc;
1257     a->alloc = size;
1258     for (; i < a->alloc; i++) {
1259       a->dp[i] = 0;
1260     }
1261   }
1262   return MP_OKAY;
1263 }
1264
1265
1266 /* b = |a|
1267  *
1268  * Simple function copies the input and fixes the sign to positive
1269  */
1270 static int mp_abs (mp_int * a, mp_int * b)
1271 {
1272   int     res;
1273
1274   /* copy a to b */
1275   if (a != b) {
1276      if ((res = mp_copy (a, b)) != MP_OKAY) {
1277        return res;
1278      }
1279   }
1280
1281   /* force the sign of b to positive */
1282   b->sign = MP_ZPOS;
1283
1284   return MP_OKAY;
1285 }
1286
1287
1288 /* set to a digit */
1289 static void mp_set (mp_int * a, mp_digit b)
1290 {
1291   mp_zero (a);
1292   a->dp[0] = b & MP_MASK;
1293   a->used  = (a->dp[0] != 0) ? 1 : 0;
1294 }
1295
1296
1297 #ifndef LTM_NO_NEG_EXP
1298 /* b = a/2 */
1299 static int mp_div_2(mp_int * a, mp_int * b)
1300 {
1301   int     x, res, oldused;
1302
1303   /* copy */
1304   if (b->alloc < a->used) {
1305     if ((res = mp_grow (b, a->used)) != MP_OKAY) {
1306       return res;
1307     }
1308   }
1309
1310   oldused = b->used;
1311   b->used = a->used;
1312   {
1313     register mp_digit r, rr, *tmpa, *tmpb;
1314
1315     /* source alias */
1316     tmpa = a->dp + b->used - 1;
1317
1318     /* dest alias */
1319     tmpb = b->dp + b->used - 1;
1320
1321     /* carry */
1322     r = 0;
1323     for (x = b->used - 1; x >= 0; x--) {
1324       /* get the carry for the next iteration */
1325       rr = *tmpa & 1;
1326
1327       /* shift the current digit, add in carry and store */
1328       *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1329
1330       /* forward carry to next iteration */
1331       r = rr;
1332     }
1333
1334     /* zero excess digits */
1335     tmpb = b->dp + b->used;
1336     for (x = b->used; x < oldused; x++) {
1337       *tmpb++ = 0;
1338     }
1339   }
1340   b->sign = a->sign;
1341   mp_clamp (b);
1342   return MP_OKAY;
1343 }
1344 #endif /* LTM_NO_NEG_EXP */
1345
1346
1347 /* shift left by a certain bit count */
1348 static int mp_mul_2d (mp_int * a, int b, mp_int * c)
1349 {
1350   mp_digit d;
1351   int      res;
1352
1353   /* copy */
1354   if (a != c) {
1355      if ((res = mp_copy (a, c)) != MP_OKAY) {
1356        return res;
1357      }
1358   }
1359
1360   if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) {
1361      if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
1362        return res;
1363      }
1364   }
1365
1366   /* shift by as many digits in the bit count */
1367   if (b >= (int)DIGIT_BIT) {
1368     if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
1369       return res;
1370     }
1371   }
1372
1373   /* shift any bit count < DIGIT_BIT */
1374   d = (mp_digit) (b % DIGIT_BIT);
1375   if (d != 0) {
1376     register mp_digit *tmpc, shift, mask, r, rr;
1377     register int x;
1378
1379     /* bitmask for carries */
1380     mask = (((mp_digit)1) << d) - 1;
1381
1382     /* shift for msbs */
1383     shift = DIGIT_BIT - d;
1384
1385     /* alias */
1386     tmpc = c->dp;
1387
1388     /* carry */
1389     r    = 0;
1390     for (x = 0; x < c->used; x++) {
1391       /* get the higher bits of the current word */
1392       rr = (*tmpc >> shift) & mask;
1393
1394       /* shift the current word and OR in the carry */
1395       *tmpc = ((*tmpc << d) | r) & MP_MASK;
1396       ++tmpc;
1397
1398       /* set the carry to the carry bits of the current word */
1399       r = rr;
1400     }
1401
1402     /* set final carry */
1403     if (r != 0) {
1404        c->dp[(c->used)++] = r;
1405     }
1406   }
1407   mp_clamp (c);
1408   return MP_OKAY;
1409 }
1410
1411
1412 static int mp_init_multi(mp_int *mp, ...)
1413 {
1414     mp_err res = MP_OKAY;      /* Assume ok until proven otherwise */
1415     int n = 0;                 /* Number of ok inits */
1416     mp_int* cur_arg = mp;
1417     va_list args;
1418
1419     va_start(args, mp);        /* init args to next argument from caller */
1420     while (cur_arg != NULL) {
1421         if (mp_init(cur_arg) != MP_OKAY) {
1422             /* Oops - error! Back-track and mp_clear what we already
1423                succeeded in init-ing, then return error.
1424             */
1425             va_list clean_args;
1426
1427             /* end the current list */
1428             va_end(args);
1429
1430             /* now start cleaning up */
1431             cur_arg = mp;
1432             va_start(clean_args, mp);
1433             while (n--) {
1434                 mp_clear(cur_arg);
1435                 cur_arg = va_arg(clean_args, mp_int*);
1436             }
1437             va_end(clean_args);
1438             res = MP_MEM;
1439             break;
1440         }
1441         n++;
1442         cur_arg = va_arg(args, mp_int*);
1443     }
1444     va_end(args);
1445     return res;                /* Assumed ok, if error flagged above. */
1446 }
1447
1448
1449 static void mp_clear_multi(mp_int *mp, ...)
1450 {
1451     mp_int* next_mp = mp;
1452     va_list args;
1453     va_start(args, mp);
1454     while (next_mp != NULL) {
1455         mp_clear(next_mp);
1456         next_mp = va_arg(args, mp_int*);
1457     }
1458     va_end(args);
1459 }
1460
1461
1462 /* shift left a certain amount of digits */
1463 static int mp_lshd (mp_int * a, int b)
1464 {
1465   int     x, res;
1466
1467   /* if its less than zero return */
1468   if (b <= 0) {
1469     return MP_OKAY;
1470   }
1471
1472   /* grow to fit the new digits */
1473   if (a->alloc < a->used + b) {
1474      if ((res = mp_grow (a, a->used + b)) != MP_OKAY) {
1475        return res;
1476      }
1477   }
1478
1479   {
1480     register mp_digit *top, *bottom;
1481
1482     /* increment the used by the shift amount then copy upwards */
1483     a->used += b;
1484
1485     /* top */
1486     top = a->dp + a->used - 1;
1487
1488     /* base */
1489     bottom = a->dp + a->used - 1 - b;
1490
1491     /* much like mp_rshd this is implemented using a sliding window
1492      * except the window goes the otherway around.  Copying from
1493      * the bottom to the top.  see bn_mp_rshd.c for more info.
1494      */
1495     for (x = a->used - 1; x >= b; x--) {
1496       *top-- = *bottom--;
1497     }
1498
1499     /* zero the lower digits */
1500     top = a->dp;
1501     for (x = 0; x < b; x++) {
1502       *top++ = 0;
1503     }
1504   }
1505   return MP_OKAY;
1506 }
1507
1508
1509 /* returns the number of bits in an int */
1510 static int mp_count_bits (mp_int * a)
1511 {
1512   int     r;
1513   mp_digit q;
1514
1515   /* shortcut */
1516   if (a->used == 0) {
1517     return 0;
1518   }
1519
1520   /* get number of digits and add that */
1521   r = (a->used - 1) * DIGIT_BIT;
1522
1523   /* take the last digit and count the bits in it */
1524   q = a->dp[a->used - 1];
1525   while (q > ((mp_digit) 0)) {
1526     ++r;
1527     q >>= ((mp_digit) 1);
1528   }
1529   return r;
1530 }
1531
1532
1533 /* calc a value mod 2**b */
1534 static int mp_mod_2d (mp_int * a, int b, mp_int * c)
1535 {
1536   int     x, res;
1537
1538   /* if b is <= 0 then zero the int */
1539   if (b <= 0) {
1540     mp_zero (c);
1541     return MP_OKAY;
1542   }
1543
1544   /* if the modulus is larger than the value than return */
1545   if (b >= (int) (a->used * DIGIT_BIT)) {
1546     res = mp_copy (a, c);
1547     return res;
1548   }
1549
1550   /* copy */
1551   if ((res = mp_copy (a, c)) != MP_OKAY) {
1552     return res;
1553   }
1554
1555   /* zero digits above the last digit of the modulus */
1556   for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
1557     c->dp[x] = 0;
1558   }
1559   /* clear the digit that is not completely outside/inside the modulus */
1560   c->dp[b / DIGIT_BIT] &=
1561     (mp_digit) ((((mp_digit) 1) << (((mp_digit) b) % DIGIT_BIT)) - ((mp_digit) 1));
1562   mp_clamp (c);
1563   return MP_OKAY;
1564 }
1565
1566
1567 /* slower bit-bang division... also smaller */
1568 static int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d)
1569 {
1570    mp_int ta, tb, tq, q;
1571    int    res, n, n2;
1572
1573   /* is divisor zero ? */
1574   if (mp_iszero (b) == 1) {
1575     return MP_VAL;
1576   }
1577
1578   /* if a < b then q=0, r = a */
1579   if (mp_cmp_mag (a, b) == MP_LT) {
1580     if (d != NULL) {
1581       res = mp_copy (a, d);
1582     } else {
1583       res = MP_OKAY;
1584     }
1585     if (c != NULL) {
1586       mp_zero (c);
1587     }
1588     return res;
1589   }
1590
1591   /* init our temps */
1592   if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) {
1593      return res;
1594   }
1595
1596
1597   mp_set(&tq, 1);
1598   n = mp_count_bits(a) - mp_count_bits(b);
1599   if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
1600       ((res = mp_abs(b, &tb)) != MP_OKAY) ||
1601       ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
1602       ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
1603       goto LBL_ERR;
1604   }
1605
1606   while (n-- >= 0) {
1607      if (mp_cmp(&tb, &ta) != MP_GT) {
1608         if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
1609             ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
1610            goto LBL_ERR;
1611         }
1612      }
1613      if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
1614          ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
1615            goto LBL_ERR;
1616      }
1617   }
1618
1619   /* now q == quotient and ta == remainder */
1620   n  = a->sign;
1621   n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
1622   if (c != NULL) {
1623      mp_exch(c, &q);
1624      c->sign  = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
1625   }
1626   if (d != NULL) {
1627      mp_exch(d, &ta);
1628      d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
1629   }
1630 LBL_ERR:
1631    mp_clear_multi(&ta, &tb, &tq, &q, NULL);
1632    return res;
1633 }
1634
1635
1636 #ifdef MP_LOW_MEM
1637    #define TAB_SIZE 32
1638 #else
1639    #define TAB_SIZE 256
1640 #endif
1641
1642 static int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
1643 {
1644   mp_int  M[TAB_SIZE], res, mu;
1645   mp_digit buf;
1646   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
1647   int (*redux)(mp_int*,mp_int*,mp_int*);
1648
1649   /* find window size */
1650   x = mp_count_bits (X);
1651   if (x <= 7) {
1652     winsize = 2;
1653   } else if (x <= 36) {
1654     winsize = 3;
1655   } else if (x <= 140) {
1656     winsize = 4;
1657   } else if (x <= 450) {
1658     winsize = 5;
1659   } else if (x <= 1303) {
1660     winsize = 6;
1661   } else if (x <= 3529) {
1662     winsize = 7;
1663   } else {
1664     winsize = 8;
1665   }
1666
1667 #ifdef MP_LOW_MEM
1668     if (winsize > 5) {
1669        winsize = 5;
1670     }
1671 #endif
1672
1673   /* init M array */
1674   /* init first cell */
1675   if ((err = mp_init(&M[1])) != MP_OKAY) {
1676      return err;
1677   }
1678
1679   /* now init the second half of the array */
1680   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
1681     if ((err = mp_init(&M[x])) != MP_OKAY) {
1682       for (y = 1<<(winsize-1); y < x; y++) {
1683         mp_clear (&M[y]);
1684       }
1685       mp_clear(&M[1]);
1686       return err;
1687     }
1688   }
1689
1690   /* create mu, used for Barrett reduction */
1691   if ((err = mp_init (&mu)) != MP_OKAY) {
1692     goto LBL_M;
1693   }
1694
1695   if (redmode == 0) {
1696      if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
1697         goto LBL_MU;
1698      }
1699      redux = mp_reduce;
1700   } else {
1701      if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
1702         goto LBL_MU;
1703      }
1704      redux = mp_reduce_2k_l;
1705   }
1706
1707   /* create M table
1708    *
1709    * The M table contains powers of the base,
1710    * e.g. M[x] = G**x mod P
1711    *
1712    * The first half of the table is not
1713    * computed though accept for M[0] and M[1]
1714    */
1715   if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
1716     goto LBL_MU;
1717   }
1718
1719   /* compute the value at M[1<<(winsize-1)] by squaring
1720    * M[1] (winsize-1) times
1721    */
1722   if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
1723     goto LBL_MU;
1724   }
1725
1726   for (x = 0; x < (winsize - 1); x++) {
1727     /* square it */
1728     if ((err = mp_sqr (&M[1 << (winsize - 1)],
1729                        &M[1 << (winsize - 1)])) != MP_OKAY) {
1730       goto LBL_MU;
1731     }
1732
1733     /* reduce modulo P */
1734     if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
1735       goto LBL_MU;
1736     }
1737   }
1738
1739   /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
1740    * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
1741    */
1742   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
1743     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
1744       goto LBL_MU;
1745     }
1746     if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
1747       goto LBL_MU;
1748     }
1749   }
1750
1751   /* setup result */
1752   if ((err = mp_init (&res)) != MP_OKAY) {
1753     goto LBL_MU;
1754   }
1755   mp_set (&res, 1);
1756
1757   /* set initial mode and bit cnt */
1758   mode   = 0;
1759   bitcnt = 1;
1760   buf    = 0;
1761   digidx = X->used - 1;
1762   bitcpy = 0;
1763   bitbuf = 0;
1764
1765   for (;;) {
1766     /* grab next digit as required */
1767     if (--bitcnt == 0) {
1768       /* if digidx == -1 we are out of digits */
1769       if (digidx == -1) {
1770         break;
1771       }
1772       /* read next digit and reset the bitcnt */
1773       buf    = X->dp[digidx--];
1774       bitcnt = (int) DIGIT_BIT;
1775     }
1776
1777     /* grab the next msb from the exponent */
1778     y     = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
1779     buf <<= (mp_digit)1;
1780
1781     /* if the bit is zero and mode == 0 then we ignore it
1782      * These represent the leading zero bits before the first 1 bit
1783      * in the exponent.  Technically this opt is not required but it
1784      * does lower the # of trivial squaring/reductions used
1785      */
1786     if (mode == 0 && y == 0) {
1787       continue;
1788     }
1789
1790     /* if the bit is zero and mode == 1 then we square */
1791     if (mode == 1 && y == 0) {
1792       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
1793         goto LBL_RES;
1794       }
1795       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
1796         goto LBL_RES;
1797       }
1798       continue;
1799     }
1800
1801     /* else we add it to the window */
1802     bitbuf |= (y << (winsize - ++bitcpy));
1803     mode    = 2;
1804
1805     if (bitcpy == winsize) {
1806       /* ok window is filled so square as required and multiply  */
1807       /* square first */
1808       for (x = 0; x < winsize; x++) {
1809         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
1810           goto LBL_RES;
1811         }
1812         if ((err = redux (&res, P, &mu)) != MP_OKAY) {
1813           goto LBL_RES;
1814         }
1815       }
1816
1817       /* then multiply */
1818       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
1819         goto LBL_RES;
1820       }
1821       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
1822         goto LBL_RES;
1823       }
1824
1825       /* empty window and reset */
1826       bitcpy = 0;
1827       bitbuf = 0;
1828       mode   = 1;
1829     }
1830   }
1831
1832   /* if bits remain then square/multiply */
1833   if (mode == 2 && bitcpy > 0) {
1834     /* square then multiply if the bit is set */
1835     for (x = 0; x < bitcpy; x++) {
1836       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
1837         goto LBL_RES;
1838       }
1839       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
1840         goto LBL_RES;
1841       }
1842
1843       bitbuf <<= 1;
1844       if ((bitbuf & (1 << winsize)) != 0) {
1845         /* then multiply */
1846         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
1847           goto LBL_RES;
1848         }
1849         if ((err = redux (&res, P, &mu)) != MP_OKAY) {
1850           goto LBL_RES;
1851         }
1852       }
1853     }
1854   }
1855
1856   mp_exch (&res, Y);
1857   err = MP_OKAY;
1858 LBL_RES:mp_clear (&res);
1859 LBL_MU:mp_clear (&mu);
1860 LBL_M:
1861   mp_clear(&M[1]);
1862   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
1863     mp_clear (&M[x]);
1864   }
1865   return err;
1866 }
1867
1868
1869 /* computes b = a*a */
1870 static int mp_sqr (mp_int * a, mp_int * b)
1871 {
1872   int     res;
1873
1874 #ifdef BN_MP_TOOM_SQR_C
1875   /* use Toom-Cook? */
1876   if (a->used >= TOOM_SQR_CUTOFF) {
1877     res = mp_toom_sqr(a, b);
1878   /* Karatsuba? */
1879   } else
1880 #endif
1881 #ifdef BN_MP_KARATSUBA_SQR_C
1882 if (a->used >= KARATSUBA_SQR_CUTOFF) {
1883     res = mp_karatsuba_sqr (a, b);
1884   } else
1885 #endif
1886   {
1887 #ifdef BN_FAST_S_MP_SQR_C
1888     /* can we use the fast comba multiplier? */
1889     if ((a->used * 2 + 1) < MP_WARRAY &&
1890          a->used <
1891          (1 << (sizeof(mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) {
1892       res = fast_s_mp_sqr (a, b);
1893     } else
1894 #endif
1895 #ifdef BN_S_MP_SQR_C
1896       res = s_mp_sqr (a, b);
1897 #else
1898 #error mp_sqr could fail
1899       res = MP_VAL;
1900 #endif
1901   }
1902   b->sign = MP_ZPOS;
1903   return res;
1904 }
1905
1906
1907 /* reduces a modulo n where n is of the form 2**p - d
1908    This differs from reduce_2k since "d" can be larger
1909    than a single digit.
1910 */
1911 static int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d)
1912 {
1913    mp_int q;
1914    int    p, res;
1915
1916    if ((res = mp_init(&q)) != MP_OKAY) {
1917       return res;
1918    }
1919
1920    p = mp_count_bits(n);
1921 top:
1922    /* q = a/2**p, a = a mod 2**p */
1923    if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
1924       goto ERR;
1925    }
1926
1927    /* q = q * d */
1928    if ((res = mp_mul(&q, d, &q)) != MP_OKAY) {
1929       goto ERR;
1930    }
1931
1932    /* a = a + q */
1933    if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
1934       goto ERR;
1935    }
1936
1937    if (mp_cmp_mag(a, n) != MP_LT) {
1938       s_mp_sub(a, n, a);
1939       goto top;
1940    }
1941
1942 ERR:
1943    mp_clear(&q);
1944    return res;
1945 }
1946
1947
1948 /* determines the setup value */
1949 static int mp_reduce_2k_setup_l(mp_int *a, mp_int *d)
1950 {
1951    int    res;
1952    mp_int tmp;
1953
1954    if ((res = mp_init(&tmp)) != MP_OKAY) {
1955       return res;
1956    }
1957
1958    if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
1959       goto ERR;
1960    }
1961
1962    if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) {
1963       goto ERR;
1964    }
1965
1966 ERR:
1967    mp_clear(&tmp);
1968    return res;
1969 }
1970
1971
1972 /* computes a = 2**b
1973  *
1974  * Simple algorithm which zeroes the int, grows it then just sets one bit
1975  * as required.
1976  */
1977 static int mp_2expt (mp_int * a, int b)
1978 {
1979   int     res;
1980
1981   /* zero a as per default */
1982   mp_zero (a);
1983
1984   /* grow a to accomodate the single bit */
1985   if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
1986     return res;
1987   }
1988
1989   /* set the used count of where the bit will go */
1990   a->used = b / DIGIT_BIT + 1;
1991
1992   /* put the single bit in its place */
1993   a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT);
1994
1995   return MP_OKAY;
1996 }
1997
1998
1999 /* pre-calculate the value required for Barrett reduction
2000  * For a given modulus "b" it calulates the value required in "a"
2001  */
2002 static int mp_reduce_setup (mp_int * a, mp_int * b)
2003 {
2004   int     res;
2005
2006   if ((res = mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) {
2007     return res;
2008   }
2009   return mp_div (a, b, a, NULL);
2010 }
2011
2012
2013 /* reduces x mod m, assumes 0 < x < m**2, mu is
2014  * precomputed via mp_reduce_setup.
2015  * From HAC pp.604 Algorithm 14.42
2016  */
2017 static int mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
2018 {
2019   mp_int  q;
2020   int     res, um = m->used;
2021
2022   /* q = x */
2023   if ((res = mp_init_copy (&q, x)) != MP_OKAY) {
2024     return res;
2025   }
2026
2027   /* q1 = x / b**(k-1)  */
2028   mp_rshd (&q, um - 1);
2029
2030   /* according to HAC this optimization is ok */
2031   if (((unsigned long) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
2032     if ((res = mp_mul (&q, mu, &q)) != MP_OKAY) {
2033       goto CLEANUP;
2034     }
2035   } else {
2036 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
2037     if ((res = s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
2038       goto CLEANUP;
2039     }
2040 #elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C)
2041     if ((res = fast_s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
2042       goto CLEANUP;
2043     }
2044 #else
2045     {
2046 #error mp_reduce would always fail
2047       res = MP_VAL;
2048       goto CLEANUP;
2049     }
2050 #endif
2051   }
2052
2053   /* q3 = q2 / b**(k+1) */
2054   mp_rshd (&q, um + 1);
2055
2056   /* x = x mod b**(k+1), quick (no division) */
2057   if ((res = mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
2058     goto CLEANUP;
2059   }
2060
2061   /* q = q * m mod b**(k+1), quick (no division) */
2062   if ((res = s_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) {
2063     goto CLEANUP;
2064   }
2065
2066   /* x = x - q */
2067   if ((res = mp_sub (x, &q, x)) != MP_OKAY) {
2068     goto CLEANUP;
2069   }
2070
2071   /* If x < 0, add b**(k+1) to it */
2072   if (mp_cmp_d (x, 0) == MP_LT) {
2073     mp_set (&q, 1);
2074     if ((res = mp_lshd (&q, um + 1)) != MP_OKAY) {
2075       goto CLEANUP;
2076     }
2077     if ((res = mp_add (x, &q, x)) != MP_OKAY) {
2078       goto CLEANUP;
2079     }
2080   }
2081
2082   /* Back off if it's too big */
2083   while (mp_cmp (x, m) != MP_LT) {
2084     if ((res = s_mp_sub (x, m, x)) != MP_OKAY) {
2085       goto CLEANUP;
2086     }
2087   }
2088
2089 CLEANUP:
2090   mp_clear (&q);
2091
2092   return res;
2093 }
2094
2095
2096 /* multiplies |a| * |b| and only computes upto digs digits of result
2097  * HAC pp. 595, Algorithm 14.12  Modified so you can control how
2098  * many digits of output are created.
2099  */
2100 static int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
2101 {
2102   mp_int  t;
2103   int     res, pa, pb, ix, iy;
2104   mp_digit u;
2105   mp_word r;
2106   mp_digit tmpx, *tmpt, *tmpy;
2107
2108   /* can we use the fast multiplier? */
2109   if (((digs) < MP_WARRAY) &&
2110       MIN (a->used, b->used) <
2111           (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
2112     return fast_s_mp_mul_digs (a, b, c, digs);
2113   }
2114
2115   if ((res = mp_init_size (&t, digs)) != MP_OKAY) {
2116     return res;
2117   }
2118   t.used = digs;
2119
2120   /* compute the digits of the product directly */
2121   pa = a->used;
2122   for (ix = 0; ix < pa; ix++) {
2123     /* set the carry to zero */
2124     u = 0;
2125
2126     /* limit ourselves to making digs digits of output */
2127     pb = MIN (b->used, digs - ix);
2128
2129     /* setup some aliases */
2130     /* copy of the digit from a used within the nested loop */
2131     tmpx = a->dp[ix];
2132
2133     /* an alias for the destination shifted ix places */
2134     tmpt = t.dp + ix;
2135
2136     /* an alias for the digits of b */
2137     tmpy = b->dp;
2138
2139     /* compute the columns of the output and propagate the carry */
2140     for (iy = 0; iy < pb; iy++) {
2141       /* compute the column as a mp_word */
2142       r       = ((mp_word)*tmpt) +
2143                 ((mp_word)tmpx) * ((mp_word)*tmpy++) +
2144                 ((mp_word) u);
2145
2146       /* the new column is the lower part of the result */
2147       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
2148
2149       /* get the carry word from the result */
2150       u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
2151     }
2152     /* set carry if it is placed below digs */
2153     if (ix + iy < digs) {
2154       *tmpt = u;
2155     }
2156   }
2157
2158   mp_clamp (&t);
2159   mp_exch (&t, c);
2160
2161   mp_clear (&t);
2162   return MP_OKAY;
2163 }
2164
2165
2166 /* Fast (comba) multiplier
2167  *
2168  * This is the fast column-array [comba] multiplier.  It is
2169  * designed to compute the columns of the product first
2170  * then handle the carries afterwards.  This has the effect
2171  * of making the nested loops that compute the columns very
2172  * simple and schedulable on super-scalar processors.
2173  *
2174  * This has been modified to produce a variable number of
2175  * digits of output so if say only a half-product is required
2176  * you don't have to compute the upper half (a feature
2177  * required for fast Barrett reduction).
2178  *
2179  * Based on Algorithm 14.12 on pp.595 of HAC.
2180  *
2181  */
2182 static int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
2183 {
2184   int     olduse, res, pa, ix, iz;
2185   mp_digit W[MP_WARRAY];
2186   register mp_word  _W;
2187
2188   /* grow the destination as required */
2189   if (c->alloc < digs) {
2190     if ((res = mp_grow (c, digs)) != MP_OKAY) {
2191       return res;
2192     }
2193   }
2194
2195   /* number of output digits to produce */
2196   pa = MIN(digs, a->used + b->used);
2197
2198   /* clear the carry */
2199   _W = 0;
2200   for (ix = 0; ix < pa; ix++) {
2201       int      tx, ty;
2202       int      iy;
2203       mp_digit *tmpx, *tmpy;
2204
2205       /* get offsets into the two bignums */
2206       ty = MIN(b->used-1, ix);
2207       tx = ix - ty;
2208
2209       /* setup temp aliases */
2210       tmpx = a->dp + tx;
2211       tmpy = b->dp + ty;
2212
2213       /* this is the number of times the loop will iterrate, essentially
2214          while (tx++ < a->used && ty-- >= 0) { ... }
2215        */
2216       iy = MIN(a->used-tx, ty+1);
2217
2218       /* execute loop */
2219       for (iz = 0; iz < iy; ++iz) {
2220          _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
2221
2222       }
2223
2224       /* store term */
2225       W[ix] = ((mp_digit)_W) & MP_MASK;
2226
2227       /* make next carry */
2228       _W = _W >> ((mp_word)DIGIT_BIT);
2229  }
2230
2231   /* setup dest */
2232   olduse  = c->used;
2233   c->used = pa;
2234
2235   {
2236     register mp_digit *tmpc;
2237     tmpc = c->dp;
2238     for (ix = 0; ix < pa+1; ix++) {
2239       /* now extract the previous digit [below the carry] */
2240       *tmpc++ = W[ix];
2241     }
2242
2243     /* clear unused digits [that existed in the old copy of c] */
2244     for (; ix < olduse; ix++) {
2245       *tmpc++ = 0;
2246     }
2247   }
2248   mp_clamp (c);
2249   return MP_OKAY;
2250 }
2251
2252
2253 /* init an mp_init for a given size */
2254 static int mp_init_size (mp_int * a, int size)
2255 {
2256   int x;
2257
2258   /* pad size so there are always extra digits */
2259   size += (MP_PREC * 2) - (size % MP_PREC);
2260
2261   /* alloc mem */
2262   a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * size);
2263   if (a->dp == NULL) {
2264     return MP_MEM;
2265   }
2266
2267   /* set the members */
2268   a->used  = 0;
2269   a->alloc = size;
2270   a->sign  = MP_ZPOS;
2271
2272   /* zero the digits */
2273   for (x = 0; x < size; x++) {
2274       a->dp[x] = 0;
2275   }
2276
2277   return MP_OKAY;
2278 }
2279
2280
2281 /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
2282 static int s_mp_sqr (mp_int * a, mp_int * b)
2283 {
2284   mp_int  t;
2285   int     res, ix, iy, pa;
2286   mp_word r;
2287   mp_digit u, tmpx, *tmpt;
2288
2289   pa = a->used;
2290   if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
2291     return res;
2292   }
2293
2294   /* default used is maximum possible size */
2295   t.used = 2*pa + 1;
2296
2297   for (ix = 0; ix < pa; ix++) {
2298     /* first calculate the digit at 2*ix */
2299     /* calculate double precision result */
2300     r = ((mp_word) t.dp[2*ix]) +
2301         ((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);
2302
2303     /* store lower part in result */
2304     t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));
2305
2306     /* get the carry */
2307     u           = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
2308
2309     /* left hand side of A[ix] * A[iy] */
2310     tmpx        = a->dp[ix];
2311
2312     /* alias for where to store the results */
2313     tmpt        = t.dp + (2*ix + 1);
2314
2315     for (iy = ix + 1; iy < pa; iy++) {
2316       /* first calculate the product */
2317       r       = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
2318
2319       /* now calculate the double precision result, note we use
2320        * addition instead of *2 since it's easier to optimize
2321        */
2322       r       = ((mp_word) *tmpt) + r + r + ((mp_word) u);
2323
2324       /* store lower part */
2325       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
2326
2327       /* get carry */
2328       u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
2329     }
2330     /* propagate upwards */
2331     while (u != ((mp_digit) 0)) {
2332       r       = ((mp_word) *tmpt) + ((mp_word) u);
2333       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
2334       u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
2335     }
2336   }
2337
2338   mp_clamp (&t);
2339   mp_exch (&t, b);
2340   mp_clear (&t);
2341   return MP_OKAY;
2342 }
2343
2344
2345 /* multiplies |a| * |b| and does not compute the lower digs digits
2346  * [meant to get the higher part of the product]
2347  */
2348 static int s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
2349 {
2350   mp_int  t;
2351   int     res, pa, pb, ix, iy;
2352   mp_digit u;
2353   mp_word r;
2354   mp_digit tmpx, *tmpt, *tmpy;
2355
2356   /* can we use the fast multiplier? */
2357 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
2358   if (((a->used + b->used + 1) < MP_WARRAY)
2359       && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
2360     return fast_s_mp_mul_high_digs (a, b, c, digs);
2361   }
2362 #endif
2363
2364   if ((res = mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) {
2365     return res;
2366   }
2367   t.used = a->used + b->used + 1;
2368
2369   pa = a->used;
2370   pb = b->used;
2371   for (ix = 0; ix < pa; ix++) {
2372     /* clear the carry */
2373     u = 0;
2374
2375     /* left hand side of A[ix] * B[iy] */
2376     tmpx = a->dp[ix];
2377
2378     /* alias to the address of where the digits will be stored */
2379     tmpt = &(t.dp[digs]);
2380
2381     /* alias for where to read the right hand side from */
2382     tmpy = b->dp + (digs - ix);
2383
2384     for (iy = digs - ix; iy < pb; iy++) {
2385       /* calculate the double precision result */
2386       r       = ((mp_word)*tmpt) +
2387                 ((mp_word)tmpx) * ((mp_word)*tmpy++) +
2388                 ((mp_word) u);
2389
2390       /* get the lower part */
2391       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
2392
2393       /* carry the carry */
2394       u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
2395     }
2396     *tmpt = u;
2397   }
2398   mp_clamp (&t);
2399   mp_exch (&t, c);
2400   mp_clear (&t);
2401   return MP_OKAY;
2402 }
2403
2404
2405 #ifdef BN_MP_MONTGOMERY_SETUP_C
2406 /* setups the montgomery reduction stuff */
2407 static int
2408 mp_montgomery_setup (mp_int * n, mp_digit * rho)
2409 {
2410   mp_digit x, b;
2411
2412 /* fast inversion mod 2**k
2413  *
2414  * Based on the fact that
2415  *
2416  * XA = 1 (mod 2**n)  =>  (X(2-XA)) A = 1 (mod 2**2n)
2417  *                    =>  2*X*A - X*X*A*A = 1
2418  *                    =>  2*(1) - (1)     = 1
2419  */
2420   b = n->dp[0];
2421
2422   if ((b & 1) == 0) {
2423     return MP_VAL;
2424   }
2425
2426   x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
2427   x *= 2 - b * x;               /* here x*a==1 mod 2**8 */
2428 #if !defined(MP_8BIT)
2429   x *= 2 - b * x;               /* here x*a==1 mod 2**16 */
2430 #endif
2431 #if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT))
2432   x *= 2 - b * x;               /* here x*a==1 mod 2**32 */
2433 #endif
2434 #ifdef MP_64BIT
2435   x *= 2 - b * x;               /* here x*a==1 mod 2**64 */
2436 #endif
2437
2438   /* rho = -1/m mod b */
2439   *rho = (unsigned long)(((mp_word)1 << ((mp_word) DIGIT_BIT)) - x) & MP_MASK;
2440
2441   return MP_OKAY;
2442 }
2443 #endif
2444
2445
2446 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
2447 /* computes xR**-1 == x (mod N) via Montgomery Reduction
2448  *
2449  * This is an optimized implementation of montgomery_reduce
2450  * which uses the comba method to quickly calculate the columns of the
2451  * reduction.
2452  *
2453  * Based on Algorithm 14.32 on pp.601 of HAC.
2454 */
2455 int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
2456 {
2457   int     ix, res, olduse;
2458   mp_word W[MP_WARRAY];
2459
2460   /* get old used count */
2461   olduse = x->used;
2462
2463   /* grow a as required */
2464   if (x->alloc < n->used + 1) {
2465     if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) {
2466       return res;
2467     }
2468   }
2469
2470   /* first we have to get the digits of the input into
2471    * an array of double precision words W[...]
2472    */
2473   {
2474     register mp_word *_W;
2475     register mp_digit *tmpx;
2476
2477     /* alias for the W[] array */
2478     _W   = W;
2479
2480     /* alias for the digits of  x*/
2481     tmpx = x->dp;
2482
2483     /* copy the digits of a into W[0..a->used-1] */
2484     for (ix = 0; ix < x->used; ix++) {
2485       *_W++ = *tmpx++;
2486     }
2487
2488     /* zero the high words of W[a->used..m->used*2] */
2489     for (; ix < n->used * 2 + 1; ix++) {
2490       *_W++ = 0;
2491     }
2492   }
2493
2494   /* now we proceed to zero successive digits
2495    * from the least significant upwards
2496    */
2497   for (ix = 0; ix < n->used; ix++) {
2498     /* mu = ai * m' mod b
2499      *
2500      * We avoid a double precision multiplication (which isn't required)
2501      * by casting the value down to a mp_digit.  Note this requires
2502      * that W[ix-1] have  the carry cleared (see after the inner loop)
2503      */
2504     register mp_digit mu;
2505     mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
2506
2507     /* a = a + mu * m * b**i
2508      *
2509      * This is computed in place and on the fly.  The multiplication
2510      * by b**i is handled by offseting which columns the results
2511      * are added to.
2512      *
2513      * Note the comba method normally doesn't handle carries in the
2514      * inner loop In this case we fix the carry from the previous
2515      * column since the Montgomery reduction requires digits of the
2516      * result (so far) [see above] to work.  This is
2517      * handled by fixing up one carry after the inner loop.  The
2518      * carry fixups are done in order so after these loops the
2519      * first m->used words of W[] have the carries fixed
2520      */
2521     {
2522       register int iy;
2523       register mp_digit *tmpn;
2524       register mp_word *_W;
2525
2526       /* alias for the digits of the modulus */
2527       tmpn = n->dp;
2528
2529       /* Alias for the columns set by an offset of ix */
2530       _W = W + ix;
2531
2532       /* inner loop */
2533       for (iy = 0; iy < n->used; iy++) {
2534           *_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
2535       }
2536     }
2537
2538     /* now fix carry for next digit, W[ix+1] */
2539     W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
2540   }
2541
2542   /* now we have to propagate the carries and
2543    * shift the words downward [all those least
2544    * significant digits we zeroed].
2545    */
2546   {
2547     register mp_digit *tmpx;
2548     register mp_word *_W, *_W1;
2549
2550     /* nox fix rest of carries */
2551
2552     /* alias for current word */
2553     _W1 = W + ix;
2554
2555     /* alias for next word, where the carry goes */
2556     _W = W + ++ix;
2557
2558     for (; ix <= n->used * 2 + 1; ix++) {
2559       *_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
2560     }
2561
2562     /* copy out, A = A/b**n
2563      *
2564      * The result is A/b**n but instead of converting from an
2565      * array of mp_word to mp_digit than calling mp_rshd
2566      * we just copy them in the right order
2567      */
2568
2569     /* alias for destination word */
2570     tmpx = x->dp;
2571
2572     /* alias for shifted double precision result */
2573     _W = W + n->used;
2574
2575     for (ix = 0; ix < n->used + 1; ix++) {
2576       *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
2577     }
2578
2579     /* zero oldused digits, if the input a was larger than
2580      * m->used+1 we'll have to clear the digits
2581      */
2582     for (; ix < olduse; ix++) {
2583       *tmpx++ = 0;
2584     }
2585   }
2586
2587   /* set the max used and clamp */
2588   x->used = n->used + 1;
2589   mp_clamp (x);
2590
2591   /* if A >= m then A = A - m */
2592   if (mp_cmp_mag (x, n) != MP_LT) {
2593     return s_mp_sub (x, n, x);
2594   }
2595   return MP_OKAY;
2596 }
2597 #endif
2598
2599
2600 #ifdef BN_MP_MUL_2_C
2601 /* b = a*2 */
2602 static int mp_mul_2(mp_int * a, mp_int * b)
2603 {
2604   int     x, res, oldused;
2605
2606   /* grow to accomodate result */
2607   if (b->alloc < a->used + 1) {
2608     if ((res = mp_grow (b, a->used + 1)) != MP_OKAY) {
2609       return res;
2610     }
2611   }
2612
2613   oldused = b->used;
2614   b->used = a->used;
2615
2616   {
2617     register mp_digit r, rr, *tmpa, *tmpb;
2618
2619     /* alias for source */
2620     tmpa = a->dp;
2621
2622     /* alias for dest */
2623     tmpb = b->dp;
2624
2625     /* carry */
2626     r = 0;
2627     for (x = 0; x < a->used; x++) {
2628
2629       /* get what will be the *next* carry bit from the
2630        * MSB of the current digit
2631        */
2632       rr = *tmpa >> ((mp_digit)(DIGIT_BIT - 1));
2633
2634       /* now shift up this digit, add in the carry [from the previous] */
2635       *tmpb++ = ((*tmpa++ << ((mp_digit)1)) | r) & MP_MASK;
2636
2637       /* copy the carry that would be from the source
2638        * digit into the next iteration
2639        */
2640       r = rr;
2641     }
2642
2643     /* new leading digit? */
2644     if (r != 0) {
2645       /* add a MSB which is always 1 at this point */
2646       *tmpb = 1;
2647       ++(b->used);
2648     }
2649
2650     /* now zero any excess digits on the destination
2651      * that we didn't write to
2652      */
2653     tmpb = b->dp + b->used;
2654     for (x = b->used; x < oldused; x++) {
2655       *tmpb++ = 0;
2656     }
2657   }
2658   b->sign = a->sign;
2659   return MP_OKAY;
2660 }
2661 #endif
2662
2663
2664 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
2665 /*
2666  * shifts with subtractions when the result is greater than b.
2667  *
2668  * The method is slightly modified to shift B unconditionally upto just under
2669  * the leading bit of b.  This saves alot of multiple precision shifting.
2670  */
2671 static int mp_montgomery_calc_normalization (mp_int * a, mp_int * b)
2672 {
2673   int     x, bits, res;
2674
2675   /* how many bits of last digit does b use */
2676   bits = mp_count_bits (b) % DIGIT_BIT;
2677
2678   if (b->used > 1) {
2679      if ((res = mp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) != MP_OKAY) {
2680         return res;
2681      }
2682   } else {
2683      mp_set(a, 1);
2684      bits = 1;
2685   }
2686
2687
2688   /* now compute C = A * B mod b */
2689   for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
2690     if ((res = mp_mul_2 (a, a)) != MP_OKAY) {
2691       return res;
2692     }
2693     if (mp_cmp_mag (a, b) != MP_LT) {
2694       if ((res = s_mp_sub (a, b, a)) != MP_OKAY) {
2695         return res;
2696       }
2697     }
2698   }
2699
2700   return MP_OKAY;
2701 }
2702 #endif
2703
2704
2705 #ifdef BN_MP_EXPTMOD_FAST_C
2706 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
2707  *
2708  * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
2709  * The value of k changes based on the size of the exponent.
2710  *
2711  * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
2712  */
2713
2714 static int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
2715 {
2716   mp_int  M[TAB_SIZE], res;
2717   mp_digit buf, mp;
2718   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
2719
2720   /* use a pointer to the reduction algorithm.  This allows us to use
2721    * one of many reduction algorithms without modding the guts of
2722    * the code with if statements everywhere.
2723    */
2724   int     (*redux)(mp_int*,mp_int*,mp_digit);
2725
2726   /* find window size */
2727   x = mp_count_bits (X);
2728   if (x <= 7) {
2729     winsize = 2;
2730   } else if (x <= 36) {
2731     winsize = 3;
2732   } else if (x <= 140) {
2733     winsize = 4;
2734   } else if (x <= 450) {
2735     winsize = 5;
2736   } else if (x <= 1303) {
2737     winsize = 6;
2738   } else if (x <= 3529) {
2739     winsize = 7;
2740   } else {
2741     winsize = 8;
2742   }
2743
2744 #ifdef MP_LOW_MEM
2745   if (winsize > 5) {
2746      winsize = 5;
2747   }
2748 #endif
2749
2750   /* init M array */
2751   /* init first cell */
2752   if ((err = mp_init(&M[1])) != MP_OKAY) {
2753      return err;
2754   }
2755
2756   /* now init the second half of the array */
2757   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2758     if ((err = mp_init(&M[x])) != MP_OKAY) {
2759       for (y = 1<<(winsize-1); y < x; y++) {
2760         mp_clear (&M[y]);
2761       }
2762       mp_clear(&M[1]);
2763       return err;
2764     }
2765   }
2766
2767   /* determine and setup reduction code */
2768   if (redmode == 0) {
2769 #ifdef BN_MP_MONTGOMERY_SETUP_C
2770      /* now setup montgomery  */
2771      if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
2772         goto LBL_M;
2773      }
2774 #else
2775      err = MP_VAL;
2776      goto LBL_M;
2777 #endif
2778
2779      /* automatically pick the comba one if available (saves quite a few calls/ifs) */
2780 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
2781      if (((P->used * 2 + 1) < MP_WARRAY) &&
2782           P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
2783         redux = fast_mp_montgomery_reduce;
2784      } else
2785 #endif
2786      {
2787 #ifdef BN_MP_MONTGOMERY_REDUCE_C
2788         /* use slower baseline Montgomery method */
2789         redux = mp_montgomery_reduce;
2790 #else
2791         err = MP_VAL;
2792         goto LBL_M;
2793 #endif
2794      }
2795   } else if (redmode == 1) {
2796 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
2797      /* setup DR reduction for moduli of the form B**k - b */
2798      mp_dr_setup(P, &mp);
2799      redux = mp_dr_reduce;
2800 #else
2801      err = MP_VAL;
2802      goto LBL_M;
2803 #endif
2804   } else {
2805 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
2806      /* setup DR reduction for moduli of the form 2**k - b */
2807      if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
2808         goto LBL_M;
2809      }
2810      redux = mp_reduce_2k;
2811 #else
2812      err = MP_VAL;
2813      goto LBL_M;
2814 #endif
2815   }
2816
2817   /* setup result */
2818   if ((err = mp_init (&res)) != MP_OKAY) {
2819     goto LBL_M;
2820   }
2821
2822   /* create M table
2823    *
2824
2825    *
2826    * The first half of the table is not computed though accept for M[0] and M[1]
2827    */
2828
2829   if (redmode == 0) {
2830 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
2831      /* now we need R mod m */
2832      if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
2833        goto LBL_RES;
2834      }
2835 #else
2836      err = MP_VAL;
2837      goto LBL_RES;
2838 #endif
2839
2840      /* now set M[1] to G * R mod m */
2841      if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
2842        goto LBL_RES;
2843      }
2844   } else {
2845      mp_set(&res, 1);
2846      if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
2847         goto LBL_RES;
2848      }
2849   }
2850
2851   /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
2852   if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
2853     goto LBL_RES;
2854   }
2855
2856   for (x = 0; x < (winsize - 1); x++) {
2857     if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
2858       goto LBL_RES;
2859     }
2860     if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
2861       goto LBL_RES;
2862     }
2863   }
2864
2865   /* create upper table */
2866   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
2867     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
2868       goto LBL_RES;
2869     }
2870     if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
2871       goto LBL_RES;
2872     }
2873   }
2874
2875   /* set initial mode and bit cnt */
2876   mode   = 0;
2877   bitcnt = 1;
2878   buf    = 0;
2879   digidx = X->used - 1;
2880   bitcpy = 0;
2881   bitbuf = 0;
2882
2883   for (;;) {
2884     /* grab next digit as required */
2885     if (--bitcnt == 0) {
2886       /* if digidx == -1 we are out of digits so break */
2887       if (digidx == -1) {
2888         break;
2889       }
2890       /* read next digit and reset bitcnt */
2891       buf    = X->dp[digidx--];
2892       bitcnt = (int)DIGIT_BIT;
2893     }
2894
2895     /* grab the next msb from the exponent */
2896     y     = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
2897     buf <<= (mp_digit)1;
2898
2899     /* if the bit is zero and mode == 0 then we ignore it
2900      * These represent the leading zero bits before the first 1 bit
2901      * in the exponent.  Technically this opt is not required but it
2902      * does lower the # of trivial squaring/reductions used
2903      */
2904     if (mode == 0 && y == 0) {
2905       continue;
2906     }
2907
2908     /* if the bit is zero and mode == 1 then we square */
2909     if (mode == 1 && y == 0) {
2910       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
2911         goto LBL_RES;
2912       }
2913       if ((err = redux (&res, P, mp)) != MP_OKAY) {
2914         goto LBL_RES;
2915       }
2916       continue;
2917     }
2918
2919     /* else we add it to the window */
2920     bitbuf |= (y << (winsize - ++bitcpy));
2921     mode    = 2;
2922
2923     if (bitcpy == winsize) {
2924       /* ok window is filled so square as required and multiply  */
2925       /* square first */
2926       for (x = 0; x < winsize; x++) {
2927         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
2928           goto LBL_RES;
2929         }
2930         if ((err = redux (&res, P, mp)) != MP_OKAY) {
2931           goto LBL_RES;
2932         }
2933       }
2934
2935       /* then multiply */
2936       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
2937         goto LBL_RES;
2938       }
2939       if ((err = redux (&res, P, mp)) != MP_OKAY) {
2940         goto LBL_RES;
2941       }
2942
2943       /* empty window and reset */
2944       bitcpy = 0;
2945       bitbuf = 0;
2946       mode   = 1;
2947     }
2948   }
2949
2950   /* if bits remain then square/multiply */
2951   if (mode == 2 && bitcpy > 0) {
2952     /* square then multiply if the bit is set */
2953     for (x = 0; x < bitcpy; x++) {
2954       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
2955         goto LBL_RES;
2956       }
2957       if ((err = redux (&res, P, mp)) != MP_OKAY) {
2958         goto LBL_RES;
2959       }
2960
2961       /* get next bit of the window */
2962       bitbuf <<= 1;
2963       if ((bitbuf & (1 << winsize)) != 0) {
2964         /* then multiply */
2965         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
2966           goto LBL_RES;
2967         }
2968         if ((err = redux (&res, P, mp)) != MP_OKAY) {
2969           goto LBL_RES;
2970         }
2971       }
2972     }
2973   }
2974
2975   if (redmode == 0) {
2976      /* fixup result if Montgomery reduction is used
2977       * recall that any value in a Montgomery system is
2978       * actually multiplied by R mod n.  So we have
2979       * to reduce one more time to cancel out the factor
2980       * of R.
2981       */
2982      if ((err = redux(&res, P, mp)) != MP_OKAY) {
2983        goto LBL_RES;
2984      }
2985   }
2986
2987   /* swap res with Y */
2988   mp_exch (&res, Y);
2989   err = MP_OKAY;
2990 LBL_RES:mp_clear (&res);
2991 LBL_M:
2992   mp_clear(&M[1]);
2993   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2994     mp_clear (&M[x]);
2995   }
2996   return err;
2997 }
2998 #endif
2999
3000
3001 #ifdef BN_FAST_S_MP_SQR_C
3002 /* the jist of squaring...
3003  * you do like mult except the offset of the tmpx [one that
3004  * starts closer to zero] can't equal the offset of tmpy.
3005  * So basically you set up iy like before then you min it with
3006  * (ty-tx) so that it never happens.  You double all those
3007  * you add in the inner loop
3008
3009 After that loop you do the squares and add them in.
3010 */
3011
3012 static int fast_s_mp_sqr (mp_int * a, mp_int * b)
3013 {
3014   int       olduse, res, pa, ix, iz;
3015   mp_digit   W[MP_WARRAY], *tmpx;
3016   mp_word   W1;
3017
3018   /* grow the destination as required */
3019   pa = a->used + a->used;
3020   if (b->alloc < pa) {
3021     if ((res = mp_grow (b, pa)) != MP_OKAY) {
3022       return res;
3023     }
3024   }
3025
3026   /* number of output digits to produce */
3027   W1 = 0;
3028   for (ix = 0; ix < pa; ix++) {
3029       int      tx, ty, iy;
3030       mp_word  _W;
3031       mp_digit *tmpy;
3032
3033       /* clear counter */
3034       _W = 0;
3035
3036       /* get offsets into the two bignums */
3037       ty = MIN(a->used-1, ix);
3038       tx = ix - ty;
3039
3040       /* setup temp aliases */
3041       tmpx = a->dp + tx;
3042       tmpy = a->dp + ty;
3043
3044       /* this is the number of times the loop will iterrate, essentially
3045          while (tx++ < a->used && ty-- >= 0) { ... }
3046        */
3047       iy = MIN(a->used-tx, ty+1);
3048
3049       /* now for squaring tx can never equal ty
3050        * we halve the distance since they approach at a rate of 2x
3051        * and we have to round because odd cases need to be executed
3052        */
3053       iy = MIN(iy, (ty-tx+1)>>1);
3054
3055       /* execute loop */
3056       for (iz = 0; iz < iy; iz++) {
3057          _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
3058       }
3059
3060       /* double the inner product and add carry */
3061       _W = _W + _W + W1;
3062
3063       /* even columns have the square term in them */
3064       if ((ix&1) == 0) {
3065          _W += ((mp_word)a->dp[ix>>1])*((mp_word)a->dp[ix>>1]);
3066       }
3067
3068       /* store it */
3069       W[ix] = (mp_digit)(_W & MP_MASK);
3070
3071       /* make next carry */
3072       W1 = _W >> ((mp_word)DIGIT_BIT);
3073   }
3074
3075   /* setup dest */
3076   olduse  = b->used;
3077   b->used = a->used+a->used;
3078
3079   {
3080     mp_digit *tmpb;
3081     tmpb = b->dp;
3082     for (ix = 0; ix < pa; ix++) {
3083       *tmpb++ = W[ix] & MP_MASK;
3084     }
3085
3086     /* clear unused digits [that existed in the old copy of c] */
3087     for (; ix < olduse; ix++) {
3088       *tmpb++ = 0;
3089     }
3090   }
3091   mp_clamp (b);
3092   return MP_OKAY;
3093 }
3094 #endif