luna - wolfssl/wolfcrypt/src/integer.c

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