networking/tls_pstm.c

   1 /*
   2  * Copyright (C) 2017 Denys Vlasenko
   3  *
   4  * Licensed under GPLv2, see file LICENSE in this source tree.
   5  */
   6 #include "tls.h"
   7
   8 /* The file is taken almost verbatim from matrixssl-3-7-2b-open/crypto/math/.
   9  * Changes are flagged with //bbox
  10  */
  11
  12 /**
  13  *      @file    pstm.c
  14  *      @version 33ef80f (HEAD, tag: MATRIXSSL-3-7-2-OPEN, tag: MATRIXSSL-3-7-2-COMM, origin/master, origin/HEAD, master)
  15  *
  16  *      Multiprecision number implementation.
  17  */
  18 /*
  19  *      Copyright (c) 2013-2015 INSIDE Secure Corporation
  20  *      Copyright (c) PeerSec Networks, 2002-2011
  21  *      All Rights Reserved
  22  *
  23  *      The latest version of this code is available at http://www.matrixssl.org
  24  *
  25  *      This software is open source; you can redistribute it and/or modify
  26  *      it under the terms of the GNU General Public License as published by
  27  *      the Free Software Foundation; either version 2 of the License, or
  28  *      (at your option) any later version.
  29  *
  30  *      This General Public License does NOT permit incorporating this software
  31  *      into proprietary programs.  If you are unable to comply with the GPL, a
  32  *      commercial license for this software may be purchased from INSIDE at
  33  *      http://www.insidesecure.com/eng/Company/Locations
  34  *
  35  *      This program is distributed in WITHOUT ANY WARRANTY; without even the
  36  *      implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  37  *      See the GNU General Public License for more details.
  38  *
  39  *      You should have received a copy of the GNU General Public License
  40  *      along with this program; if not, write to the Free Software
  41  *      Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  42  *      http://www.gnu.org/copyleft/gpl.html
  43  */
  44 /******************************************************************************/
  45
  46 //bbox
  47 //#include "../cryptoApi.h"
  48 #ifndef DISABLE_PSTM
  49
  50 #undef pstm_mul_2d
  51 static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c); //bbox: was int16 b
  52 #define pstm_mul_2d(a, b, c) (pstm_mul_2d(a, b, c), PSTM_OKAY)
  53
  54 /******************************************************************************/
  55 /*
  56         init an pstm_int for a given size
  57  */
  58 #undef pstm_init_size
  59 #define pstm_init_size(pool, a, size) \
  60         pstm_init_size(      a, size)
  61 int32 FAST_FUNC pstm_init_size(psPool_t *pool, pstm_int * a, uint32 size)
  62 {
  63 //bbox
  64 //      uint16          x;
  65
  66 /*
  67         alloc mem
  68  */
  69         a->dp = xzalloc(sizeof (pstm_digit) * size);//bbox
  70 //bbox  a->pool = pool;
  71         a->used  = 0;
  72         a->alloc = size;
  73         a->sign  = PSTM_ZPOS;
  74 /*
  75         zero the digits
  76  */
  77 //bbox
  78 //      for (x = 0; x < size; x++) {
  79 //              a->dp[x] = 0;
  80 //      }
  81         return PSTM_OKAY;
  82 }
  83 #undef pstm_init_size
  84 #define pstm_init_size(pool, a, size) (pstm_init_size(a, size), PSTM_OKAY)
  85
  86 /******************************************************************************/
  87 /*
  88         Init a new pstm_int.
  89 */
  90 #undef pstm_init
  91 #define pstm_init(pool, a) \
  92         pstm_init(      a)
  93 static int32 pstm_init(psPool_t *pool, pstm_int * a)
  94 {
  95 //bbox
  96 //      int32           i;
  97 /*
  98         allocate memory required and clear it
  99  */
 100         a->dp = xzalloc(sizeof (pstm_digit) * PSTM_DEFAULT_INIT);//bbox
 101 /*
 102         set the digits to zero
 103  */
 104 //bbox
 105 //      for (i = 0; i < PSTM_DEFAULT_INIT; i++) {
 106 //              a->dp[i] = 0;
 107 //      }
 108 /*
 109         set the used to zero, allocated digits to the default precision and sign
 110         to positive
 111  */
 112 //bbox  a->pool = pool;
 113         a->used  = 0;
 114         a->alloc = PSTM_DEFAULT_INIT;
 115         a->sign  = PSTM_ZPOS;
 116
 117         return PSTM_OKAY;
 118 }
 119 #undef pstm_init
 120 #define pstm_init(pool, a) (pstm_init(a), PSTM_OKAY)
 121
 122 /******************************************************************************/
 123 /*
 124         Grow as required
 125  */
 126 #undef pstm_grow
 127 int32 FAST_FUNC pstm_grow(pstm_int * a, int size)
 128 {
 129         int                     i; //bbox: was int16
 130         pstm_digit              *tmp;
 131
 132 /*
 133         If the alloc size is smaller alloc more ram.
 134  */
 135         if (a->alloc < size) {
 136 /*
 137                 Reallocate the array a->dp
 138
 139                 We store the return in a temporary variable in case the operation
 140                 failed we don't want to overwrite the dp member of a.
 141 */
 142                 tmp = xrealloc(a->dp, sizeof (pstm_digit) * size);//bbox
 143 /*
 144                 reallocation succeeded so set a->dp
 145  */
 146                 a->dp = tmp;
 147 /*
 148                 zero excess digits
 149  */
 150                 i                       = a->alloc;
 151                 a->alloc        = size;
 152                 for (; i < a->alloc; i++) {
 153                         a->dp[i] = 0;
 154                 }
 155         }
 156         return PSTM_OKAY;
 157 }
 158 #define pstm_grow(a, size) (pstm_grow(a, size), PSTM_OKAY)
 159
 160 /******************************************************************************/
 161 /*
 162         copy, b = a (b must be pre-allocated)
 163  */
 164 #undef pstm_copy
 165 int32 pstm_copy(pstm_int * a, pstm_int * b)
 166 {
 167         int32   res, n;
 168
 169 /*
 170         If dst == src do nothing
 171  */
 172         if (a == b) {
 173                 return PSTM_OKAY;
 174         }
 175 /*
 176         Grow dest
 177  */
 178         if (b->alloc < a->used) {
 179                 if ((res = pstm_grow (b, a->used)) != PSTM_OKAY) {
 180                         return res;
 181                 }
 182         }
 183 /*
 184         Zero b and copy the parameters over
 185  */
 186         {
 187                 register pstm_digit *tmpa, *tmpb;
 188
 189                 /* pointer aliases */
 190                 /* source */
 191                 tmpa = a->dp;
 192
 193                 /* destination */
 194                 tmpb = b->dp;
 195
 196                 /* copy all the digits */
 197                 for (n = 0; n < a->used; n++) {
 198                         *tmpb++ = *tmpa++;
 199                 }
 200
 201                 /* clear high digits */
 202                 for (; n < b->used; n++) {
 203                         *tmpb++ = 0;
 204                 }
 205         }
 206 /*
 207         copy used count and sign
 208  */
 209         b->used = a->used;
 210         b->sign = a->sign;
 211         return PSTM_OKAY;
 212 }
 213 #define pstm_copy(a, b) (pstm_copy(a, b), PSTM_OKAY)
 214
 215 /******************************************************************************/
 216 /*
 217         Trim unused digits
 218
 219         This is used to ensure that leading zero digits are trimed and the
 220         leading "used" digit will be non-zero. Typically very fast.  Also fixes
 221         the sign if there are no more leading digits
 222 */
 223 void FAST_FUNC pstm_clamp(pstm_int * a)
 224 {
 225 /*      decrease used while the most significant digit is zero. */
 226         while (a->used > 0 && a->dp[a->used - 1] == 0) {
 227                 --(a->used);
 228         }
 229 /*      reset the sign flag if used == 0 */
 230         if (a->used == 0) {
 231                 a->sign = PSTM_ZPOS;
 232         }
 233 }
 234
 235 /******************************************************************************/
 236 /*
 237         clear one (frees).
 238  */
 239 void FAST_FUNC pstm_clear(pstm_int * a)
 240 {
 241         int32           i;
 242 /*
 243         only do anything if a hasn't been freed previously
 244  */
 245         if (a != NULL && a->dp != NULL) {
 246 /*
 247                 first zero the digits
 248  */
 249                 for (i = 0; i < a->used; i++) {
 250                         a->dp[i] = 0;
 251                 }
 252
 253                 psFree (a->dp, a->pool);
 254 /*
 255                 reset members to make debugging easier
 256  */
 257                 a->dp           = NULL;
 258                 a->alloc        = a->used = 0;
 259                 a->sign         = PSTM_ZPOS;
 260         }
 261 }
 262
 263 /******************************************************************************/
 264 /*
 265         clear many (frees).
 266  */
 267 #if 0 //UNUSED
 268 void pstm_clear_multi(pstm_int *mp0, pstm_int *mp1, pstm_int *mp2,
 269                                         pstm_int *mp3, pstm_int *mp4, pstm_int *mp5,
 270                                         pstm_int *mp6, pstm_int *mp7)
 271 {
 272         int32           n;              /* Number of ok inits */
 273
 274         pstm_int        *tempArray[9];
 275
 276         tempArray[0] = mp0;
 277         tempArray[1] = mp1;
 278         tempArray[2] = mp2;
 279         tempArray[3] = mp3;
 280         tempArray[4] = mp4;
 281         tempArray[5] = mp5;
 282         tempArray[6] = mp6;
 283         tempArray[7] = mp7;
 284         tempArray[8] = NULL;
 285
 286         for (n = 0; tempArray[n] != NULL; n++) {
 287                 if ((tempArray[n] != NULL) && (tempArray[n]->dp != NULL)) {
 288                         pstm_clear(tempArray[n]);
 289                 }
 290         }
 291 }
 292 #endif
 293
 294 /******************************************************************************/
 295 /*
 296         Set to zero.
 297  */
 298 static void pstm_zero(pstm_int * a)
 299 {
 300         int32           n;
 301         pstm_digit      *tmp;
 302
 303         a->sign = PSTM_ZPOS;
 304         a->used = 0;
 305
 306         tmp = a->dp;
 307         for (n = 0; n < a->alloc; n++) {
 308                 *tmp++ = 0;
 309         }
 310 }
 311
 312
 313 /******************************************************************************/
 314 /*
 315         Compare maginitude of two ints (unsigned).
 316  */
 317 int32 FAST_FUNC pstm_cmp_mag(pstm_int * a, pstm_int * b)
 318 {
 319         int             n; //bbox: was int16
 320         pstm_digit      *tmpa, *tmpb;
 321
 322 /*
 323         compare based on # of non-zero digits
 324  */
 325         if (a->used > b->used) {
 326                 return PSTM_GT;
 327         }
 328
 329         if (a->used < b->used) {
 330                 return PSTM_LT;
 331         }
 332
 333         /* alias for a */
 334         tmpa = a->dp + (a->used - 1);
 335
 336         /* alias for b */
 337         tmpb = b->dp + (a->used - 1);
 338
 339 /*
 340         compare based on digits
 341  */
 342         for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
 343                 if (*tmpa > *tmpb) {
 344                         return PSTM_GT;
 345                 }
 346                 if (*tmpa < *tmpb) {
 347                         return PSTM_LT;
 348                 }
 349         }
 350         return PSTM_EQ;
 351 }
 352
 353 /******************************************************************************/
 354 /*
 355         Compare two ints (signed)
 356  */
 357 int32 FAST_FUNC pstm_cmp(pstm_int * a, pstm_int * b)
 358 {
 359 /*
 360         compare based on sign
 361  */
 362         if (a->sign != b->sign) {
 363                 if (a->sign == PSTM_NEG) {
 364                         return PSTM_LT;
 365                 } else {
 366                         return PSTM_GT;
 367                 }
 368         }
 369 /*
 370         compare digits
 371  */
 372         if (a->sign == PSTM_NEG) {
 373                 /* if negative compare opposite direction */
 374                 return pstm_cmp_mag(b, a);
 375         } else {
 376                 return pstm_cmp_mag(a, b);
 377         }
 378 }
 379
 380 /******************************************************************************/
 381 /*
 382         pstm_ints can be initialized more precisely when they will populated
 383         using pstm_read_unsigned_bin since the length of the byte stream is known
 384 */
 385 int32 FAST_FUNC pstm_init_for_read_unsigned_bin(psPool_t *pool, pstm_int *a, uint32 len)
 386 {
 387         int32 size;
 388 /*
 389         Need to set this based on how many words max it will take to store the bin.
 390         The magic + 2:
 391                 1 to round up for the remainder of this integer math
 392                 1 for the initial carry of '1' bits that fall between DIGIT_BIT and 8
 393 */
 394         size = (((len / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
 395                 / DIGIT_BIT) + 2;
 396         return pstm_init_size(pool, a, size);
 397 }
 398
 399
 400 /******************************************************************************/
 401 /*
 402         Reads a unsigned char array into pstm_int format.  User should have
 403         called pstm_init_for_read_unsigned_bin first.  There is some grow logic
 404         here if the default pstm_init was used but we don't really want to hit it.
 405 */
 406 int32 FAST_FUNC pstm_read_unsigned_bin(pstm_int *a, unsigned char *b, int32 c)
 407 {
 408         /* zero the int */
 409         pstm_zero (a);
 410
 411 /*
 412         If we know the endianness of this architecture, and we're using
 413         32-bit pstm_digits, we can optimize this
 414 */
 415 #if (defined(ENDIAN_LITTLE) || defined(ENDIAN_BIG)) && !defined(PSTM_64BIT)
 416   /* But not for both simultaneously */
 417 #if defined(ENDIAN_LITTLE) && defined(ENDIAN_BIG)
 418 #error Both ENDIAN_LITTLE and ENDIAN_BIG defined.
 419 #endif
 420         {
 421                 unsigned char *pd;
 422                 if ((unsigned)c > (PSTM_MAX_SIZE * sizeof(pstm_digit))) {
 423                         uint32 excess = c - (PSTM_MAX_SIZE * sizeof(pstm_digit));
 424                         c -= excess;
 425                         b += excess;
 426                 }
 427                 a->used = ((c + sizeof(pstm_digit) - 1)/sizeof(pstm_digit));
 428                 if (a->alloc < a->used) {
 429                         if (pstm_grow(a, a->used) != PSTM_OKAY) {
 430                                 return PSTM_MEM;
 431                         }
 432                 }
 433                 pd = (unsigned char *)a->dp;
 434                 /* read the bytes in */
 435 #ifdef ENDIAN_BIG
 436                 {
 437                         /* Use Duff's device to unroll the loop. */
 438                         int32 idx = (c - 1) & ~3;
 439                         switch (c % 4) {
 440                                 case 0: do { pd[idx+0] = *b++;
 441                                         case 3:      pd[idx+1] = *b++;
 442                                         case 2:      pd[idx+2] = *b++;
 443                                         case 1:      pd[idx+3] = *b++;
 444                                         idx -= 4;
 445                                 } while ((c -= 4) > 0);
 446                         }
 447                 }
 448 #else
 449                 for (c -= 1; c >= 0; c -= 1) {
 450                         pd[c] = *b++;
 451                 }
 452 #endif
 453         }
 454 #else
 455         /* Big enough based on the len? */
 456         a->used = (((c / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
 457                 / DIGIT_BIT) + 2;
 458
 459         if (a->alloc < a->used) {
 460                 if (pstm_grow(a, a->used) != PSTM_OKAY) {
 461                         return PSTM_MEM;
 462                 }
 463         }
 464         /* read the bytes in */
 465         for (; c > 0; c--) {
 466                 if (pstm_mul_2d (a, 8, a) != PSTM_OKAY) {
 467                         return PS_MEM_FAIL;
 468                 }
 469                 a->dp[0] |= *b++;
 470                 a->used += 1;
 471         }
 472 #endif
 473
 474         pstm_clamp (a);
 475         return PS_SUCCESS;
 476 }
 477
 478 /******************************************************************************/
 479 /*
 480 */
 481 static int pstm_count_bits(pstm_int * a)
 482 {
 483         int     r; //bbox: was int16
 484         pstm_digit q;
 485
 486         if (a->used == 0) {
 487                 return 0;
 488         }
 489
 490         /* get number of digits and add that */
 491         r = (a->used - 1) * DIGIT_BIT;
 492
 493         /* take the last digit and count the bits in it */
 494         q = a->dp[a->used - 1];
 495         while (q > ((pstm_digit) 0)) {
 496                 ++r;
 497                 q >>= ((pstm_digit) 1);
 498         }
 499         return r;
 500 }
 501
 502 /******************************************************************************/
 503 int32 FAST_FUNC pstm_unsigned_bin_size(pstm_int *a)
 504 {
 505         int32     size = pstm_count_bits (a);
 506         return (size / 8 + ((size & 7) != 0 ? 1 : 0));
 507 }
 508
 509 /******************************************************************************/
 510 static void pstm_set(pstm_int *a, pstm_digit b)
 511 {
 512    pstm_zero(a);
 513    a->dp[0] = b;
 514    a->used  = a->dp[0] ? 1 : 0;
 515 }
 516
 517 /******************************************************************************/
 518 /*
 519         Right shift
 520 */
 521 static void pstm_rshd(pstm_int *a, int x)
 522 {
 523         int y; //bbox: was int16
 524
 525         /* too many digits just zero and return */
 526         if (x >= a->used) {
 527                 pstm_zero(a);
 528                 return;
 529         }
 530
 531         /* shift */
 532         for (y = 0; y < a->used - x; y++) {
 533                 a->dp[y] = a->dp[y+x];
 534         }
 535
 536         /* zero rest */
 537         for (; y < a->used; y++) {
 538                 a->dp[y] = 0;
 539         }
 540
 541         /* decrement count */
 542         a->used -= x;
 543         pstm_clamp(a);
 544 }
 545
 546 /******************************************************************************/
 547 /*
 548         Shift left a certain amount of digits.
 549  */
 550 #undef pstm_lshd
 551 static int32 pstm_lshd(pstm_int * a, int b)
 552 {
 553         int     x; //bbox: was int16
 554         int32   res;
 555
 556 /*
 557         If its less than zero return.
 558  */
 559         if (b <= 0) {
 560                 return PSTM_OKAY;
 561         }
 562 /*
 563         Grow to fit the new digits.
 564  */
 565         if (a->alloc < a->used + b) {
 566                 if ((res = pstm_grow (a, a->used + b)) != PSTM_OKAY) {
 567                         return res;
 568                 }
 569         }
 570
 571         {
 572                 register pstm_digit *top, *bottom;
 573 /*
 574                 Increment the used by the shift amount then copy upwards.
 575  */
 576                 a->used += b;
 577
 578                 /* top */
 579                 top = a->dp + a->used - 1;
 580
 581                 /* base */
 582                 bottom = a->dp + a->used - 1 - b;
 583 /*
 584                 This is implemented using a sliding window except the window goes the
 585                 other way around.  Copying from the bottom to the top.
 586  */
 587                 for (x = a->used - 1; x >= b; x--) {
 588                         *top-- = *bottom--;
 589                 }
 590
 591                 /* zero the lower digits */
 592                 top = a->dp;
 593                 for (x = 0; x < b; x++) {
 594                         *top++ = 0;
 595                 }
 596         }
 597         return PSTM_OKAY;
 598 }
 599 #define pstm_lshd(a, b) (pstm_lshd(a, b), PSTM_OKAY)
 600
 601 /******************************************************************************/
 602 /*
 603         computes a = 2**b
 604 */
 605 static int32 pstm_2expt(pstm_int *a, int b)
 606 {
 607         int     z; //bbox: was int16
 608
 609    /* zero a as per default */
 610         pstm_zero (a);
 611
 612         if (b < 0) {
 613                 return PSTM_OKAY;
 614         }
 615
 616         z = b / DIGIT_BIT;
 617         if (z >= PSTM_MAX_SIZE) {
 618                 return PS_LIMIT_FAIL;
 619         }
 620
 621         /* set the used count of where the bit will go */
 622         a->used = z + 1;
 623
 624         if (a->used > a->alloc) {
 625                 if (pstm_grow(a, a->used) != PSTM_OKAY) {
 626                         return PS_MEM_FAIL;
 627                 }
 628         }
 629
 630         /* put the single bit in its place */
 631         a->dp[z] = ((pstm_digit)1) << (b % DIGIT_BIT);
 632         return PSTM_OKAY;
 633 }
 634
 635 /******************************************************************************/
 636 /*
 637
 638 */
 639 int32 FAST_FUNC pstm_mul_2(pstm_int * a, pstm_int * b)
 640 {
 641         int32   res;
 642         int     x, oldused; //bbox: was int16
 643
 644 /*
 645         grow to accomodate result
 646  */
 647         if (b->alloc < a->used + 1) {
 648                 if ((res = pstm_grow (b, a->used + 1)) != PSTM_OKAY) {
 649                         return res;
 650                 }
 651         }
 652         oldused = b->used;
 653         b->used = a->used;
 654
 655         {
 656                 register pstm_digit r, rr, *tmpa, *tmpb;
 657
 658                 /* alias for source */
 659                 tmpa = a->dp;
 660
 661                 /* alias for dest */
 662                 tmpb = b->dp;
 663
 664                 /* carry */
 665                 r = 0;
 666                 for (x = 0; x < a->used; x++) {
 667 /*
 668                         get what will be the *next* carry bit from the
 669                         MSB of the current digit
 670 */
 671                         rr = *tmpa >> ((pstm_digit)(DIGIT_BIT - 1));
 672 /*
 673                         now shift up this digit, add in the carry [from the previous]
 674 */
 675                         *tmpb++ = ((*tmpa++ << ((pstm_digit)1)) | r);
 676 /*
 677                         copy the carry that would be from the source
 678                         digit into the next iteration
 679 */
 680                         r = rr;
 681                 }
 682
 683                 /* new leading digit? */
 684                 if (r != 0 && b->used != (PSTM_MAX_SIZE-1)) {
 685                         /* add a MSB which is always 1 at this point */
 686                         *tmpb = 1;
 687                         ++(b->used);
 688                 }
 689 /*
 690                 now zero any excess digits on the destination that we didn't write to
 691 */
 692                 tmpb = b->dp + b->used;
 693                 for (x = b->used; x < oldused; x++) {
 694                         *tmpb++ = 0;
 695                 }
 696         }
 697         b->sign = a->sign;
 698         return PSTM_OKAY;
 699 }
 700
 701 /******************************************************************************/
 702 /*
 703         unsigned subtraction ||a|| >= ||b|| ALWAYS!
 704 */
 705 int32 FAST_FUNC s_pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
 706 {
 707         int             oldbused, oldused; //bbox: was int16
 708         int32           x;
 709         pstm_word       t;
 710
 711         if (b->used > a->used) {
 712                 return PS_LIMIT_FAIL;
 713         }
 714         if (c->alloc < a->used) {
 715                 if ((x = pstm_grow (c, a->used)) != PSTM_OKAY) {
 716                         return x;
 717                 }
 718         }
 719         oldused  = c->used;
 720         oldbused = b->used;
 721         c->used  = a->used;
 722         t = 0;
 723
 724         for (x = 0; x < oldbused; x++) {
 725                 t = ((pstm_word)a->dp[x]) - (((pstm_word)b->dp[x]) + t);
 726                 c->dp[x] = (pstm_digit)t;
 727                 t = (t >> DIGIT_BIT)&1;
 728         }
 729         for (; x < a->used; x++) {
 730                 t = ((pstm_word)a->dp[x]) - t;
 731                 c->dp[x] = (pstm_digit)t;
 732                 t = (t >> DIGIT_BIT);
 733         }
 734         for (; x < oldused; x++) {
 735                 c->dp[x] = 0;
 736         }
 737         pstm_clamp(c);
 738         return PSTM_OKAY;
 739 }
 740
 741 /******************************************************************************/
 742 /*
 743         unsigned addition
 744 */
 745 static int32 s_pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
 746 {
 747         int                             x, y, oldused; //bbox: was int16
 748         register pstm_word      t, adp, bdp;
 749
 750         y = a->used;
 751         if (b->used > y) {
 752                 y = b->used;
 753         }
 754         oldused = c->used;
 755         c->used = y;
 756
 757         if (c->used > c->alloc) {
 758                 if (pstm_grow(c, c->used) != PSTM_OKAY) {
 759                         return PS_MEM_FAIL;
 760                 }
 761         }
 762
 763         t = 0;
 764         for (x = 0; x < y; x++) {
 765                 if (a->used < x) {
 766                         adp = 0;
 767                 } else {
 768                         adp = (pstm_word)a->dp[x];
 769                 }
 770                 if (b->used < x) {
 771                         bdp = 0;
 772                 } else {
 773                         bdp = (pstm_word)b->dp[x];
 774                 }
 775                 t         += (adp) + (bdp);
 776                 c->dp[x]   = (pstm_digit)t;
 777                 t        >>= DIGIT_BIT;
 778         }
 779         if (t != 0 && x < PSTM_MAX_SIZE) {
 780                 if (c->used == c->alloc) {
 781                         if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
 782                                 return PS_MEM_FAIL;
 783                         }
 784                 }
 785                 c->dp[c->used++] = (pstm_digit)t;
 786                 ++x;
 787         }
 788
 789         c->used = x;
 790         for (; x < oldused; x++) {
 791                 c->dp[x] = 0;
 792         }
 793         pstm_clamp(c);
 794         return PSTM_OKAY;
 795 }
 796
 797
 798 /******************************************************************************/
 799 /*
 800
 801 */
 802 int32 FAST_FUNC pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
 803 {
 804         int32   res;
 805         int     sa, sb; //bbox: was int16
 806
 807         sa = a->sign;
 808         sb = b->sign;
 809
 810         if (sa != sb) {
 811 /*
 812                 subtract a negative from a positive, OR a positive from a negative.
 813                 For both, ADD their magnitudes, and use the sign of the first number.
 814  */
 815                 c->sign = sa;
 816                 if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
 817                         return res;
 818                 }
 819         } else {
 820 /*
 821                 subtract a positive from a positive, OR a negative from a negative.
 822                 First, take the difference between their magnitudes, then...
 823  */
 824                 if (pstm_cmp_mag (a, b) != PSTM_LT) {
 825                         /* Copy the sign from the first */
 826                         c->sign = sa;
 827                         /* The first has a larger or equal magnitude */
 828                         if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
 829                                 return res;
 830                         }
 831                 } else {
 832                         /* The result has the _opposite_ sign from the first number. */
 833                         c->sign = (sa == PSTM_ZPOS) ? PSTM_NEG : PSTM_ZPOS;
 834                         /* The second has a larger magnitude */
 835                         if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
 836                                 return res;
 837                         }
 838                 }
 839         }
 840         return PS_SUCCESS;
 841 }
 842
 843 /******************************************************************************/
 844 /*
 845         c = a - b
 846 */
 847 #if 0 //UNUSED
 848 int32 pstm_sub_d(psPool_t *pool, pstm_int *a, pstm_digit b, pstm_int *c)
 849 {
 850         pstm_int        tmp;
 851         int32           res;
 852
 853         if (pstm_init_size(pool, &tmp, sizeof(pstm_digit)) != PSTM_OKAY) {
 854                 return PS_MEM_FAIL;
 855         }
 856         pstm_set(&tmp, b);
 857         res = pstm_sub(a, &tmp, c);
 858         pstm_clear(&tmp);
 859         return res;
 860 }
 861 #endif
 862
 863 /******************************************************************************/
 864 /*
 865         setups the montgomery reduction
 866 */
 867 static int32 pstm_montgomery_setup(pstm_int *a, pstm_digit *rho)
 868 {
 869         pstm_digit x, b;
 870
 871 /*
 872         fast inversion mod 2**k
 873         Based on the fact that
 874         XA = 1 (mod 2**n)       =>  (X(2-XA)) A = 1 (mod 2**2n)
 875                                                 =>  2*X*A - X*X*A*A = 1
 876                                                 =>  2*(1) - (1)     = 1
 877  */
 878         b = a->dp[0];
 879
 880         if ((b & 1) == 0) {
 881                 psTraceCrypto("pstm_montogomery_setup failure\n");
 882                 return PS_ARG_FAIL;
 883         }
 884
 885         x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
 886         x *= 2 - b * x;               /* here x*a==1 mod 2**8 */
 887         x *= 2 - b * x;               /* here x*a==1 mod 2**16 */
 888         x *= 2 - b * x;               /* here x*a==1 mod 2**32 */
 889 #ifdef PSTM_64BIT
 890         x *= 2 - b * x;               /* here x*a==1 mod 2**64 */
 891 #endif
 892         /* rho = -1/m mod b */
 893         *rho = (pstm_digit)(((pstm_word) 1 << ((pstm_word) DIGIT_BIT)) -
 894                 ((pstm_word)x));
 895         return PSTM_OKAY;
 896 }
 897
 898 /******************************************************************************/
 899 /*
 900  *      computes a = B**n mod b without division or multiplication useful for
 901  *      normalizing numbers in a Montgomery system.
 902  */
 903 static int32 pstm_montgomery_calc_normalization(pstm_int *a, pstm_int *b)
 904 {
 905         int32     x;
 906         int     bits; //bbox: was int16
 907
 908         /* how many bits of last digit does b use */
 909         bits = pstm_count_bits (b) % DIGIT_BIT;
 910         if (!bits) bits = DIGIT_BIT;
 911
 912         /* compute A = B^(n-1) * 2^(bits-1) */
 913         if (b->used > 1) {
 914                 if ((x = pstm_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) !=
 915                                 PSTM_OKAY) {
 916                         return x;
 917                 }
 918         } else {
 919                 pstm_set(a, 1);
 920                 bits = 1;
 921         }
 922
 923         /* now compute C = A * B mod b */
 924         for (x = bits - 1; x < (int32)DIGIT_BIT; x++) {
 925                 if (pstm_mul_2 (a, a) != PSTM_OKAY) {
 926                         return PS_MEM_FAIL;
 927                 }
 928                 if (pstm_cmp_mag (a, b) != PSTM_LT) {
 929                         if (s_pstm_sub (a, b, a) != PSTM_OKAY) {
 930                                 return PS_MEM_FAIL;
 931                         }
 932                 }
 933         }
 934         return PSTM_OKAY;
 935 }
 936
 937 /******************************************************************************/
 938 /*
 939         c = a * 2**d
 940 */
 941 #undef pstm_mul_2d
 942 static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c)
 943 {
 944         pstm_digit      carry, carrytmp, shift;
 945         int             x; //bbox: was int16
 946
 947         /* copy it */
 948         if (pstm_copy(a, c) != PSTM_OKAY) {
 949                 return PS_MEM_FAIL;
 950         }
 951
 952         /* handle whole digits */
 953         if (b >= DIGIT_BIT) {
 954                 if (pstm_lshd(c, b/DIGIT_BIT) != PSTM_OKAY) {
 955                         return PS_MEM_FAIL;
 956                 }
 957         }
 958         b %= DIGIT_BIT;
 959
 960         /* shift the digits */
 961         if (b != 0) {
 962                 carry = 0;
 963                 shift = DIGIT_BIT - b;
 964                 for (x = 0; x < c->used; x++) {
 965                         carrytmp = c->dp[x] >> shift;
 966                         c->dp[x] = (c->dp[x] << b) + carry;
 967                         carry = carrytmp;
 968                 }
 969                 /* store last carry if room */
 970                 if (carry && x < PSTM_MAX_SIZE) {
 971                         if (c->used == c->alloc) {
 972                                 if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
 973                                         return PS_MEM_FAIL;
 974                                 }
 975                         }
 976                         c->dp[c->used++] = carry;
 977                 }
 978         }
 979         pstm_clamp(c);
 980         return PSTM_OKAY;
 981 }
 982 #define pstm_mul_2d(a, b, c) (pstm_mul_2d(a, b, c), PSTM_OKAY)
 983
 984 /******************************************************************************/
 985 /*
 986         c = a mod 2**d
 987 */
 988 #undef pstm_mod_2d
 989 static int32 pstm_mod_2d(pstm_int *a, int b, pstm_int *c) //bbox: was int16 b
 990 {
 991         int     x; //bbox: was int16
 992
 993         /* zero if count less than or equal to zero */
 994         if (b <= 0) {
 995                 pstm_zero(c);
 996                 return PSTM_OKAY;
 997         }
 998
 999         /* get copy of input */
1000         if (pstm_copy(a, c) != PSTM_OKAY) {
1001                 return PS_MEM_FAIL;
1002         }
1003
1004         /* if 2**d is larger than we just return */
1005         if (b >= (DIGIT_BIT * a->used)) {
1006                 return PSTM_OKAY;
1007         }
1008
1009         /* zero digits above the last digit of the modulus */
1010         for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++)
1011         {
1012                 c->dp[x] = 0;
1013         }
1014         /* clear the digit that is not completely outside/inside the modulus */
1015         c->dp[b / DIGIT_BIT] &= ~((pstm_digit)0) >> (DIGIT_BIT - b);
1016         pstm_clamp (c);
1017         return PSTM_OKAY;
1018 }
1019 #define pstm_mod_2d(a, b, c) (pstm_mod_2d(a, b, c), PSTM_OKAY)
1020
1021
1022 /******************************************************************************/
1023 /*
1024         c = a * b
1025 */
1026 #undef pstm_mul_d
1027 static int32 pstm_mul_d(pstm_int *a, pstm_digit b, pstm_int *c)
1028 {
1029         pstm_word       w;
1030         int32           res;
1031         int             x, oldused; //bbox: was int16
1032
1033         if (c->alloc < a->used + 1) {
1034                 if ((res = pstm_grow (c, a->used + 1)) != PSTM_OKAY) {
1035                         return res;
1036                 }
1037         }
1038         oldused = c->used;
1039         c->used = a->used;
1040         c->sign = a->sign;
1041         w       = 0;
1042         for (x = 0; x < a->used; x++) {
1043                 w         = ((pstm_word)a->dp[x]) * ((pstm_word)b) + w;
1044                 c->dp[x]  = (pstm_digit)w;
1045                 w         = w >> DIGIT_BIT;
1046         }
1047         if (w != 0 && (a->used != PSTM_MAX_SIZE)) {
1048                 c->dp[c->used++] = (pstm_digit)w;
1049                 ++x;
1050         }
1051         for (; x < oldused; x++) {
1052                 c->dp[x] = 0;
1053         }
1054         pstm_clamp(c);
1055         return PSTM_OKAY;
1056 }
1057 #define pstm_mul_d(a, b, c) (pstm_mul_d(a, b, c), PSTM_OKAY)
1058
1059 /******************************************************************************/
1060 /*
1061         c = a / 2**b
1062 */
1063 #undef pstm_div_2d
1064 #define pstm_div_2d(pool, a, b, c, d) \
1065         pstm_div_2d(      a, b, c, d)
1066 static int32 pstm_div_2d(psPool_t *pool, pstm_int *a, int b, pstm_int *c,
1067                                         pstm_int *d)
1068 {
1069         pstm_digit      D, r, rr;
1070         int32           res;
1071         int             x; //bbox: was int16
1072         pstm_int        t;
1073
1074         /* if the shift count is <= 0 then we do no work */
1075         if (b <= 0) {
1076                 if (pstm_copy (a, c) != PSTM_OKAY) {
1077                         return PS_MEM_FAIL;
1078                 }
1079                 if (d != NULL) {
1080                         pstm_zero (d);
1081                 }
1082                 return PSTM_OKAY;
1083         }
1084
1085         /* get the remainder */
1086         if (d != NULL) {
1087                 if (pstm_init(pool, &t) != PSTM_OKAY) {
1088                         return PS_MEM_FAIL;
1089                 }
1090                 if (pstm_mod_2d (a, b, &t) != PSTM_OKAY) {
1091                         res = PS_MEM_FAIL;
1092                         goto LBL_DONE;
1093                 }
1094         }
1095
1096         /* copy */
1097         if (pstm_copy(a, c) != PSTM_OKAY) {
1098                 res = PS_MEM_FAIL;
1099                 goto LBL_DONE;
1100         }
1101
1102         /* shift by as many digits in the bit count */
1103         if (b >= (int32)DIGIT_BIT) {
1104                 pstm_rshd (c, b / DIGIT_BIT);
1105         }
1106
1107         /* shift any bit count < DIGIT_BIT */
1108         D = (pstm_digit) (b % DIGIT_BIT);
1109         if (D != 0) {
1110                 register pstm_digit *tmpc, mask, shift;
1111
1112                 /* mask */
1113                 mask = (((pstm_digit)1) << D) - 1;
1114
1115                 /* shift for lsb */
1116                 shift = DIGIT_BIT - D;
1117
1118                 /* alias */
1119                 tmpc = c->dp + (c->used - 1);
1120
1121                 /* carry */
1122                 r = 0;
1123                 for (x = c->used - 1; x >= 0; x--) {
1124                         /* get the lower  bits of this word in a temp */
1125                         rr = *tmpc & mask;
1126
1127                         /* shift the current word and mix in the carry bits from previous */
1128                         *tmpc = (*tmpc >> D) | (r << shift);
1129                         --tmpc;
1130
1131                         /* set the carry to the carry bits of the current word above */
1132                         r = rr;
1133                 }
1134         }
1135         pstm_clamp (c);
1136
1137         res = PSTM_OKAY;
1138 LBL_DONE:
1139         if (d != NULL) {
1140                 if (pstm_copy(&t, d) != PSTM_OKAY) {
1141                         res = PS_MEM_FAIL;
1142                 }
1143                 pstm_clear(&t);
1144         }
1145         return res;
1146 }
1147 #undef pstm_div_2d
1148 #define pstm_div_2d(pool, a, b, c, d) (pstm_div_2d(a, b, c, d), PSTM_OKAY)
1149
1150 /******************************************************************************/
1151 /*
1152         b = a/2
1153 */
1154 #if 0 //UNUSED
1155 int32 pstm_div_2(pstm_int * a, pstm_int * b)
1156 {
1157         int     x, oldused; //bbox: was int16
1158
1159         if (b->alloc < a->used) {
1160                 if (pstm_grow(b, a->used) != PSTM_OKAY) {
1161                         return PS_MEM_FAIL;
1162                 }
1163         }
1164         oldused = b->used;
1165         b->used = a->used;
1166         {
1167                 register pstm_digit r, rr, *tmpa, *tmpb;
1168
1169                 /* source alias */
1170                 tmpa = a->dp + b->used - 1;
1171
1172                 /* dest alias */
1173                 tmpb = b->dp + b->used - 1;
1174
1175                 /* carry */
1176                 r = 0;
1177                 for (x = b->used - 1; x >= 0; x--) {
1178                         /* get the carry for the next iteration */
1179                         rr = *tmpa & 1;
1180
1181                         /* shift the current digit, add in carry and store */
1182                         *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1183
1184                         /* forward carry to next iteration */
1185                         r = rr;
1186                 }
1187
1188                 /* zero excess digits */
1189                 tmpb = b->dp + b->used;
1190                 for (x = b->used; x < oldused; x++) {
1191                         *tmpb++ = 0;
1192                 }
1193         }
1194         b->sign = a->sign;
1195         pstm_clamp (b);
1196         return PSTM_OKAY;
1197 }
1198 #endif
1199
1200 /******************************************************************************/
1201 /*
1202         Creates "a" then copies b into it
1203  */
1204 #undef pstm_init_copy
1205 #define pstm_init_copy(pool, a, b, toSqr) \
1206         pstm_init_copy(      a, b, toSqr)
1207 static int32 pstm_init_copy(psPool_t *pool, pstm_int * a, pstm_int * b, int toSqr)
1208 {
1209         int     x; //bbox: was int16
1210         int32   res;
1211
1212         if (a == b) {
1213                 return PSTM_OKAY;
1214         }
1215         x = b->alloc;
1216
1217         if (toSqr) {
1218 /*
1219                 Smart-size:  Increasing size of a if b->used is roughly half
1220                 of b->alloc because usage has shown that a lot of these copies
1221                 go on to be squared and need these extra digits
1222 */
1223                 if ((b->used * 2) + 2 >= x) {
1224                         x = (b->used * 2) + 3;
1225                 }
1226         }
1227         if ((res = pstm_init_size(pool, a, x)) != PSTM_OKAY) {
1228                 return res;
1229         }
1230         return pstm_copy(b, a);
1231 }
1232 #undef pstm_init_copy
1233 #define pstm_init_copy(pool, a, b, toSqr) (pstm_init_copy(a, b, toSqr), PSTM_OKAY)
1234
1235 /******************************************************************************/
1236 /*
1237         With some compilers, we have seen issues linking with the builtin
1238         64 bit division routine. The issues with either manifest in a failure
1239         to find 'udivdi3' at link time, or a runtime invalid instruction fault
1240         during an RSA operation.
1241         The routine below divides a 64 bit unsigned int by a 32 bit unsigned int
1242         explicitly, rather than using the division operation
1243                 The 64 bit result is placed in the 'numerator' parameter
1244                 The 32 bit mod (remainder) of the division is the return parameter
1245         Based on implementations by:
1246                 Copyright (C) 2003 Bernardo Innocenti <bernie@develer.com>
1247                 Copyright (C) 1999 Hewlett-Packard Co
1248                 Copyright (C) 1999 David Mosberger-Tang <davidm@hpl.hp.com>
1249 */
1250 #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
1251 static uint32 psDiv64(uint64 *numerator, uint32 denominator)
1252 {
1253         uint64  rem = *numerator;
1254         uint64  b = denominator;
1255         uint64  res = 0;
1256         uint64  d = 1;
1257         uint32  high = rem >> 32;
1258
1259         if (high >= denominator) {
1260                 high /= denominator;
1261                 res = (uint64) high << 32;
1262                 rem -= (uint64) (high * denominator) << 32;
1263         }
1264         while ((int64)b > 0 && b < rem) {
1265                 b = b+b;
1266                 d = d+d;
1267         }
1268         do {
1269                 if (rem >= b) {
1270                         rem -= b;
1271                         res += d;
1272                 }
1273                 b >>= 1;
1274                 d >>= 1;
1275         } while (d);
1276         *numerator = res;
1277         return rem;
1278 }
1279 #endif /* USE_MATRIX_DIV64 */
1280
1281 #if defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1282 typedef unsigned long   uint128 __attribute__ ((mode(TI)));
1283 static uint64 psDiv128(uint128 *numerator, uint64 denominator)
1284 {
1285         uint128 rem = *numerator;
1286         uint128 b = denominator;
1287         uint128 res = 0;
1288         uint128 d = 1;
1289         uint64  high = rem >> 64;
1290
1291         if (high >= denominator) {
1292                 high /= denominator;
1293                 res = (uint128) high << 64;
1294                 rem -= (uint128) (high * denominator) << 64;
1295         }
1296         while ((uint128)b > 0 && b < rem) {
1297                 b = b+b;
1298                 d = d+d;
1299         }
1300         do {
1301                 if (rem >= b) {
1302                         rem -= b;
1303                         res += d;
1304                 }
1305                 b >>= 1;
1306                 d >>= 1;
1307         } while (d);
1308         *numerator = res;
1309         return rem;
1310 }
1311 #endif /* USE_MATRIX_DIV128 */
1312
1313 /******************************************************************************/
1314 /*
1315         a/b => cb + d == a
1316 */
1317 static int32 pstm_div(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
1318                                 pstm_int *d)
1319 {
1320         pstm_int        q, x, y, t1, t2;
1321         int32           res;
1322         int             n, t, i, norm, neg; //bbox: was int16
1323
1324         /* is divisor zero ? */
1325         if (pstm_iszero (b) == 1) {
1326                 return PS_LIMIT_FAIL;
1327         }
1328
1329         /* if a < b then q=0, r = a */
1330         if (pstm_cmp_mag (a, b) == PSTM_LT) {
1331                 if (d != NULL) {
1332                         if (pstm_copy(a, d) != PSTM_OKAY) {
1333                                 return PS_MEM_FAIL;
1334                         }
1335                 }
1336                 if (c != NULL) {
1337                         pstm_zero (c);
1338                 }
1339                 return PSTM_OKAY;
1340         }
1341 /*
1342         Smart-size inits
1343 */
1344         if ((res = pstm_init_size(pool, &t1, a->alloc)) != PSTM_OKAY) {
1345                 return res;
1346         }
1347         if ((res = pstm_init_size(pool, &t2, 3)) != PSTM_OKAY) {
1348                 goto LBL_T1;
1349         }
1350         if ((res = pstm_init_copy(pool, &x, a, 0)) != PSTM_OKAY) {
1351                 goto LBL_T2;
1352         }
1353 /*
1354         Used to be an init_copy on b but pstm_grow was always hit with triple size
1355 */
1356         if ((res = pstm_init_size(pool, &y, b->used * 3)) != PSTM_OKAY) {
1357                 goto LBL_X;
1358         }
1359         if ((res = pstm_copy(b, &y)) != PSTM_OKAY) {
1360                 goto LBL_Y;
1361         }
1362
1363         /* fix the sign */
1364         neg = (a->sign == b->sign) ? PSTM_ZPOS : PSTM_NEG;
1365         x.sign = y.sign = PSTM_ZPOS;
1366
1367         /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
1368         norm = pstm_count_bits(&y) % DIGIT_BIT;
1369         if (norm < (int32)(DIGIT_BIT-1)) {
1370                 norm = (DIGIT_BIT-1) - norm;
1371                 if ((res = pstm_mul_2d(&x, norm, &x)) != PSTM_OKAY) {
1372                         goto LBL_Y;
1373                 }
1374                 if ((res = pstm_mul_2d(&y, norm, &y)) != PSTM_OKAY) {
1375                         goto LBL_Y;
1376                 }
1377         } else {
1378                 norm = 0;
1379         }
1380
1381         /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
1382         n = x.used - 1;
1383         t = y.used - 1;
1384
1385         if ((res = pstm_init_size(pool, &q, n - t + 1)) != PSTM_OKAY) {
1386                 goto LBL_Y;
1387         }
1388         q.used = n - t + 1;
1389
1390         /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
1391         if ((res = pstm_lshd(&y, n - t)) != PSTM_OKAY) { /* y = y*b**{n-t} */
1392                 goto LBL_Q;
1393         }
1394
1395         while (pstm_cmp (&x, &y) != PSTM_LT) {
1396                 ++(q.dp[n - t]);
1397                 if ((res = pstm_sub(&x, &y, &x)) != PSTM_OKAY) {
1398                         goto LBL_Q;
1399                 }
1400         }
1401
1402         /* reset y by shifting it back down */
1403         pstm_rshd (&y, n - t);
1404
1405         /* step 3. for i from n down to (t + 1) */
1406         for (i = n; i >= (t + 1); i--) {
1407                 if (i > x.used) {
1408                         continue;
1409                 }
1410
1411                 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
1412                  * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
1413                 if (x.dp[i] == y.dp[t]) {
1414                         q.dp[i - t - 1] = (pstm_digit)((((pstm_word)1) << DIGIT_BIT) - 1);
1415                 } else {
1416                         pstm_word tmp;
1417                         tmp = ((pstm_word) x.dp[i]) << ((pstm_word) DIGIT_BIT);
1418                         tmp |= ((pstm_word) x.dp[i - 1]);
1419 #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
1420                         psDiv64(&tmp, y.dp[t]);
1421 #elif defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1422                         psDiv128(&tmp, y.dp[t]);
1423 #else
1424                         tmp /= ((pstm_word) y.dp[t]);
1425 #endif /* USE_MATRIX_DIV64 */
1426                         q.dp[i - t - 1] = (pstm_digit) (tmp);
1427                 }
1428
1429                 /* while (q{i-t-1} * (yt * b + y{t-1})) >
1430                          xi * b**2 + xi-1 * b + xi-2
1431
1432                         do q{i-t-1} -= 1;
1433                 */
1434                 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
1435                 do {
1436                         q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
1437
1438                         /* find left hand */
1439                         pstm_zero (&t1);
1440                         t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
1441                         t1.dp[1] = y.dp[t];
1442                         t1.used = 2;
1443                         if ((res = pstm_mul_d (&t1, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
1444                                 goto LBL_Q;
1445                         }
1446
1447                         /* find right hand */
1448                         t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
1449                         t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
1450                         t2.dp[2] = x.dp[i];
1451                         t2.used = 3;
1452                 } while (pstm_cmp_mag(&t1, &t2) == PSTM_GT);
1453
1454                 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
1455                 if ((res = pstm_mul_d(&y, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
1456                         goto LBL_Q;
1457                 }
1458
1459                 if ((res = pstm_lshd(&t1, i - t - 1)) != PSTM_OKAY) {
1460                         goto LBL_Q;
1461                 }
1462
1463                 if ((res = pstm_sub(&x, &t1, &x)) != PSTM_OKAY) {
1464                         goto LBL_Q;
1465                 }
1466
1467                 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
1468                 if (x.sign == PSTM_NEG) {
1469                         if ((res = pstm_copy(&y, &t1)) != PSTM_OKAY) {
1470                                 goto LBL_Q;
1471                         }
1472                         if ((res = pstm_lshd (&t1, i - t - 1)) != PSTM_OKAY) {
1473                                 goto LBL_Q;
1474                         }
1475                         if ((res = pstm_add (&x, &t1, &x)) != PSTM_OKAY) {
1476                                 goto LBL_Q;
1477                         }
1478                         q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
1479                 }
1480         }
1481 /*
1482         now q is the quotient and x is the remainder (which we have to normalize)
1483 */
1484         /* get sign before writing to c */
1485         x.sign = x.used == 0 ? PSTM_ZPOS : a->sign;
1486
1487         if (c != NULL) {
1488                 pstm_clamp (&q);
1489                 if (pstm_copy (&q, c) != PSTM_OKAY) {
1490                         res = PS_MEM_FAIL;
1491                         goto LBL_Q;
1492                 }
1493                 c->sign = neg;
1494         }
1495
1496         if (d != NULL) {
1497                 if ((res = pstm_div_2d (pool, &x, norm, &x, NULL)) != PSTM_OKAY) {
1498                         goto LBL_Q;
1499                 }
1500 /*
1501                 the following is a kludge, essentially we were seeing the right
1502                 remainder but with excess digits that should have been zero
1503  */
1504                 for (i = b->used; i < x.used; i++) {
1505                         x.dp[i] = 0;
1506                 }
1507                 pstm_clamp(&x);
1508                 if (pstm_copy (&x, d) != PSTM_OKAY) {
1509                         res = PS_MEM_FAIL;
1510                         goto LBL_Q;
1511                 }
1512         }
1513
1514         res = PSTM_OKAY;
1515
1516 LBL_Q:pstm_clear (&q);
1517 LBL_Y:pstm_clear (&y);
1518 LBL_X:pstm_clear (&x);
1519 LBL_T2:pstm_clear (&t2);
1520 LBL_T1:pstm_clear (&t1);
1521
1522         return res;
1523 }
1524
1525 /******************************************************************************/
1526 /*
1527         Swap the elements of two integers, for cases where you can't simply swap
1528         the pstm_int pointers around
1529 */
1530 static void pstm_exch(pstm_int * a, pstm_int * b)
1531 {
1532         pstm_int                t;
1533
1534         t       = *a;
1535         *a      = *b;
1536         *b      = t;
1537 }
1538
1539 /******************************************************************************/
1540 /*
1541         c = a mod b, 0 <= c < b
1542 */
1543 static int32 pstm_mod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
1544 {
1545         pstm_int        t;
1546         int32           err;
1547 /*
1548         Smart-size
1549 */
1550         if ((err = pstm_init_size(pool, &t, b->alloc)) != PSTM_OKAY) {
1551                 return err;
1552         }
1553         if ((err = pstm_div(pool, a, b, NULL, &t)) != PSTM_OKAY) {
1554                 pstm_clear (&t);
1555                 return err;
1556         }
1557         if (t.sign != b->sign) {
1558                 err = pstm_add(&t, b, c);
1559         } else {
1560                 pstm_exch (&t, c);
1561         }
1562         pstm_clear (&t);
1563         return err;
1564 }
1565
1566 /******************************************************************************/
1567 /*
1568         d = a * b (mod c)
1569 */
1570 int32 FAST_FUNC pstm_mulmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
1571                         pstm_int *d)
1572 {
1573         int32           res;
1574         int             size; //bbox: was int16
1575         pstm_int        tmp;
1576
1577 /*
1578         Smart-size pstm_inits.  d is an output that is influenced by this local 't'
1579         so don't shrink 'd' if it wants to becuase this will lead to an pstm_grow
1580         in RSA operations
1581 */
1582         size = a->used + b->used + 1;
1583         if ((a == d) && (size < a->alloc)) {
1584                 size = a->alloc;
1585         }
1586         if ((res = pstm_init_size(pool, &tmp, size)) != PSTM_OKAY) {
1587                 return res;
1588         }
1589         if ((res = pstm_mul_comba(pool, a, b, &tmp, NULL, 0)) != PSTM_OKAY) {
1590                 pstm_clear(&tmp);
1591                 return res;
1592         }
1593         res = pstm_mod(pool, &tmp, c, d);
1594         pstm_clear(&tmp);
1595         return res;
1596 }
1597
1598 /******************************************************************************/
1599 /*
1600  *      y = g**x (mod b)
1601  *      Some restrictions... x must be positive and < b
1602  */
1603 int32 FAST_FUNC pstm_exptmod(psPool_t *pool, pstm_int *G, pstm_int *X, pstm_int *P,
1604                         pstm_int *Y)
1605 {
1606         pstm_int        M[32], res; /* Keep this winsize based: (1 << max_winsize) */
1607         pstm_digit      buf, mp;
1608         pstm_digit      *paD;
1609         int32           err, bitbuf;
1610         int             bitcpy, bitcnt, mode, digidx, x, y, winsize; //bbox: was int16
1611         uint32          paDlen;
1612
1613         /* set window size from what user set as optimization */
1614         x = pstm_count_bits(X);
1615         if (x < 50) {
1616                 winsize = 2;
1617         } else {
1618                 winsize = PS_EXPTMOD_WINSIZE;
1619         }
1620
1621         /* now setup montgomery  */
1622         if ((err = pstm_montgomery_setup (P, &mp)) != PSTM_OKAY) {
1623                 return err;
1624         }
1625
1626         /* setup result */
1627         if ((err = pstm_init_size(pool, &res, (P->used * 2) + 1)) != PSTM_OKAY) {
1628                 return err;
1629         }
1630 /*
1631         create M table
1632         The M table contains powers of the input base, e.g. M[x] = G^x mod P
1633         The first half of the table is not computed though except for M[0] and M[1]
1634  */
1635         /* now we need R mod m */
1636         if ((err = pstm_montgomery_calc_normalization (&res, P)) != PSTM_OKAY) {
1637                 goto LBL_RES;
1638         }
1639 /*
1640         init M array
1641         init first cell
1642  */
1643         if ((err = pstm_init_size(pool, &M[1], res.used)) != PSTM_OKAY) {
1644                 goto LBL_RES;
1645         }
1646
1647         /* now set M[1] to G * R mod m */
1648         if (pstm_cmp_mag(P, G) != PSTM_GT) {
1649                 /* G > P so we reduce it first */
1650                 if ((err = pstm_mod(pool, G, P, &M[1])) != PSTM_OKAY) {
1651                         goto LBL_M;
1652                 }
1653         } else {
1654                 if ((err = pstm_copy(G, &M[1])) != PSTM_OKAY) {
1655                         goto LBL_M;
1656                 }
1657         }
1658         if ((err = pstm_mulmod (pool, &M[1], &res, P, &M[1])) != PSTM_OKAY) {
1659                 goto LBL_M;
1660         }
1661 /*
1662         Pre-allocated digit.  Used for mul, sqr, AND reduce
1663 */
1664         paDlen = ((M[1].used + 3) * 2) * sizeof(pstm_digit);
1665         paD = xzalloc(paDlen);//bbox
1666 /*
1667         compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times
1668  */
1669         if (pstm_init_copy(pool, &M[1 << (winsize - 1)], &M[1], 1) != PSTM_OKAY) {
1670                 err = PS_MEM_FAIL;
1671                 goto LBL_PAD;
1672         }
1673         for (x = 0; x < (winsize - 1); x++) {
1674                 if ((err = pstm_sqr_comba (pool, &M[1 << (winsize - 1)],
1675                                 &M[1 << (winsize - 1)], paD, paDlen)) != PSTM_OKAY) {
1676                         goto LBL_PAD;
1677                 }
1678                 if ((err = pstm_montgomery_reduce(pool, &M[1 << (winsize - 1)], P, mp,
1679                                 paD, paDlen)) != PSTM_OKAY) {
1680                         goto LBL_PAD;
1681                 }
1682         }
1683 /*
1684         now init the second half of the array
1685 */
1686         for (x = (1<<(winsize-1)) + 1; x < (1 << winsize); x++) {
1687                 if ((err = pstm_init_size(pool, &M[x], M[1<<(winsize-1)].alloc + 1))
1688                                 != PSTM_OKAY) {
1689                         for (y = 1<<(winsize-1); y < x; y++) {
1690                                 pstm_clear(&M[y]);
1691                         }
1692                         goto LBL_PAD;
1693                 }
1694         }
1695
1696         /* create upper table */
1697         for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
1698                 if ((err = pstm_mul_comba(pool, &M[x - 1], &M[1], &M[x], paD, paDlen))
1699                                 != PSTM_OKAY) {
1700                         goto LBL_MARRAY;
1701                 }
1702                 if ((err = pstm_montgomery_reduce(pool, &M[x], P, mp, paD, paDlen)) !=
1703                                 PSTM_OKAY) {
1704                         goto LBL_MARRAY;
1705                 }
1706         }
1707
1708         /* set initial mode and bit cnt */
1709         mode   = 0;
1710         bitcnt = 1;
1711         buf    = 0;
1712         digidx = X->used - 1;
1713         bitcpy = 0;
1714         bitbuf = 0;
1715
1716         for (;;) {
1717                 /* grab next digit as required */
1718                 if (--bitcnt == 0) {
1719                         /* if digidx == -1 we are out of digits so break */
1720                         if (digidx == -1) {
1721                                 break;
1722                         }
1723                         /* read next digit and reset bitcnt */
1724                         buf    = X->dp[digidx--];
1725                         bitcnt = (int32)DIGIT_BIT;
1726                 }
1727
1728                 /* grab the next msb from the exponent */
1729                 y     = (pstm_digit)(buf >> (DIGIT_BIT - 1)) & 1;
1730                 buf <<= (pstm_digit)1;
1731 /*
1732                  If the bit is zero and mode == 0 then we ignore it.
1733                  These represent the leading zero bits before the first 1 bit
1734                  in the exponent.  Technically this opt is not required but it
1735                  does lower the # of trivial squaring/reductions used
1736 */
1737                 if (mode == 0 && y == 0) {
1738                         continue;
1739                 }
1740
1741                 /* if the bit is zero and mode == 1 then we square */
1742                 if (mode == 1 && y == 0) {
1743                         if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1744                                         PSTM_OKAY) {
1745                                 goto LBL_MARRAY;
1746                         }
1747                         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1748                                         != PSTM_OKAY) {
1749                                 goto LBL_MARRAY;
1750                         }
1751                         continue;
1752                 }
1753
1754                 /* else we add it to the window */
1755                 bitbuf |= (y << (winsize - ++bitcpy));
1756                 mode    = 2;
1757
1758                 if (bitcpy == winsize) {
1759                         /* ok window is filled so square as required and mul square first */
1760                         for (x = 0; x < winsize; x++) {
1761                                 if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1762                                                 PSTM_OKAY) {
1763                                         goto LBL_MARRAY;
1764                                 }
1765                                 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
1766                                                 paDlen)) != PSTM_OKAY) {
1767                                         goto LBL_MARRAY;
1768                                 }
1769                         }
1770
1771                         /* then multiply */
1772                         if ((err = pstm_mul_comba(pool, &res, &M[bitbuf], &res, paD,
1773                                         paDlen)) != PSTM_OKAY) {
1774                                 goto LBL_MARRAY;
1775                         }
1776                         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1777                                         != PSTM_OKAY) {
1778                                 goto LBL_MARRAY;
1779                         }
1780
1781                         /* empty window and reset */
1782                         bitcpy = 0;
1783                         bitbuf = 0;
1784                         mode   = 1;
1785                 }
1786         }
1787
1788         /* if bits remain then square/multiply */
1789         if (mode == 2 && bitcpy > 0) {
1790                 /* square then multiply if the bit is set */
1791                 for (x = 0; x < bitcpy; x++) {
1792                         if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1793                                         PSTM_OKAY) {
1794                                 goto LBL_MARRAY;
1795                         }
1796                         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1797                                         != PSTM_OKAY) {
1798                                 goto LBL_MARRAY;
1799                         }
1800
1801                         /* get next bit of the window */
1802                         bitbuf <<= 1;
1803                         if ((bitbuf & (1 << winsize)) != 0) {
1804                         /* then multiply */
1805                                 if ((err = pstm_mul_comba(pool, &res, &M[1], &res, paD, paDlen))
1806                                                 != PSTM_OKAY) {
1807                                         goto LBL_MARRAY;
1808                                 }
1809                                 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
1810                                                 paDlen)) != PSTM_OKAY) {
1811                                         goto LBL_MARRAY;
1812                                 }
1813                         }
1814                 }
1815         }
1816 /*
1817         Fix up result if Montgomery reduction is used recall that any value in a
1818         Montgomery system is actually multiplied by R mod n.  So we have to reduce
1819         one more time to cancel out the factor of R.
1820 */
1821         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen)) !=
1822                         PSTM_OKAY) {
1823                 goto LBL_MARRAY;
1824         }
1825         /* swap res with Y */
1826         if ((err = pstm_copy (&res, Y)) != PSTM_OKAY) {
1827                 goto LBL_MARRAY;
1828         }
1829         err = PSTM_OKAY;
1830 LBL_MARRAY:
1831         for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
1832                 pstm_clear(&M[x]);
1833         }
1834 LBL_PAD:psFree(paD, pool);
1835 LBL_M: pstm_clear(&M[1]);
1836 LBL_RES:pstm_clear(&res);
1837         return err;
1838 }
1839
1840 /******************************************************************************/
1841 /*
1842
1843 */
1844 int32 FAST_FUNC pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
1845 {
1846         int32   res;
1847         int     sa, sb; //bbox: was int16
1848
1849         /* get sign of both inputs */
1850         sa = a->sign;
1851         sb = b->sign;
1852
1853         /* handle two cases, not four */
1854         if (sa == sb) {
1855                 /* both positive or both negative, add their mags, copy the sign */
1856                 c->sign = sa;
1857                 if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
1858                         return res;
1859                 }
1860         } else {
1861 /*
1862                 one positive, the other negative
1863                 subtract the one with the greater magnitude from the one of the lesser
1864                 magnitude. The result gets the sign of the one with the greater mag.
1865  */
1866                 if (pstm_cmp_mag (a, b) == PSTM_LT) {
1867                         c->sign = sb;
1868                         if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
1869                                 return res;
1870                         }
1871                 } else {
1872                         c->sign = sa;
1873                         if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
1874                                 return res;
1875                         }
1876                 }
1877         }
1878         return PS_SUCCESS;
1879 }
1880
1881 /******************************************************************************/
1882 /*
1883         reverse an array, used for radix code
1884 */
1885 static void pstm_reverse (unsigned char *s, int len) //bbox: was int16 len
1886 {
1887         int32     ix, iy;
1888         unsigned char t;
1889
1890         ix = 0;
1891         iy = len - 1;
1892         while (ix < iy) {
1893                 t     = s[ix];
1894                 s[ix] = s[iy];
1895                 s[iy] = t;
1896                 ++ix;
1897                 --iy;
1898         }
1899 }
1900 /******************************************************************************/
1901 /*
1902         No reverse.  Useful in some of the EIP-154 PKA stuff where special byte
1903         order seems to come into play more often
1904 */
1905 #if 0 //UNUSED
1906 int32 pstm_to_unsigned_bin_nr(psPool_t *pool, pstm_int *a, unsigned char *b)
1907 {
1908         int32     res;
1909         int     x; //bbox: was int16
1910         pstm_int  t = { 0 };
1911
1912         if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
1913                 return res;
1914         }
1915
1916         x = 0;
1917         while (pstm_iszero (&t) == 0) {
1918                 b[x++] = (unsigned char) (t.dp[0] & 255);
1919                 if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
1920                         pstm_clear(&t);
1921                         return res;
1922                 }
1923         }
1924         pstm_clear(&t);
1925         return PS_SUCCESS;
1926 }
1927 #endif
1928 /******************************************************************************/
1929 /*
1930
1931 */
1932 int32 FAST_FUNC pstm_to_unsigned_bin(psPool_t *pool, pstm_int *a, unsigned char *b)
1933 {
1934         int32     res;
1935         int     x; //bbox: was int16
1936         pstm_int  t = { 0 };
1937
1938         if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
1939                 return res;
1940         }
1941
1942         x = 0;
1943         while (pstm_iszero (&t) == 0) {
1944                 b[x++] = (unsigned char) (t.dp[0] & 255);
1945                 if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
1946                         pstm_clear(&t);
1947                         return res;
1948                 }
1949         }
1950         pstm_reverse (b, x);
1951         pstm_clear(&t);
1952         return PS_SUCCESS;
1953 }
1954
1955 #if 0 //UNUSED
1956 /******************************************************************************/
1957 /*
1958         compare against a single digit
1959 */
1960 static int32 pstm_cmp_d(pstm_int *a, pstm_digit b)
1961 {
1962         /* compare based on sign */
1963         if ((b && a->used == 0) || a->sign == PSTM_NEG) {
1964                 return PSTM_LT;
1965         }
1966
1967         /* compare based on magnitude */
1968         if (a->used > 1) {
1969                 return PSTM_GT;
1970         }
1971
1972         /* compare the only digit of a to b */
1973         if (a->dp[0] > b) {
1974                 return PSTM_GT;
1975         } else if (a->dp[0] < b) {
1976                 return PSTM_LT;
1977         } else {
1978                 return PSTM_EQ;
1979         }
1980 }
1981
1982 /*
1983         Need invmod for ECC and also private key loading for hardware crypto
1984         in cases where dQ > dP.  The values must be switched and a new qP must be
1985         calculated using this function
1986 */
1987 //bbox: pool unused
1988 #define pstm_invmod_slow(pool, a, b, c) \
1989         pstm_invmod_slow(      a, b, c)
1990 static int32 pstm_invmod_slow(psPool_t *pool, pstm_int * a, pstm_int * b,
1991                                 pstm_int * c)
1992 {
1993         pstm_int  x, y, u, v, A, B, C, D;
1994         int32     res;
1995
1996         /* b cannot be negative */
1997         if (b->sign == PSTM_NEG || pstm_iszero(b) == 1) {
1998                 return PS_LIMIT_FAIL;
1999         }
2000
2001         /* init temps */
2002         if (pstm_init_size(pool, &x, b->used) != PSTM_OKAY) {
2003                 return PS_MEM_FAIL;
2004         }
2005
2006         /* x = a, y = b */
2007         if ((res = pstm_mod(pool, a, b, &x)) != PSTM_OKAY) {
2008                 goto LBL_X;
2009         }
2010
2011         if (pstm_init_copy(pool, &y, b, 0) != PSTM_OKAY) {
2012                 goto LBL_X;
2013         }
2014
2015         /* 2. [modified] if x,y are both even then return an error! */
2016         if (pstm_iseven (&x) == 1 && pstm_iseven (&y) == 1) {
2017                 res = PS_FAILURE;
2018                 goto LBL_Y;
2019         }
2020
2021         /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2022         if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
2023                 goto LBL_Y;
2024         }
2025         if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
2026                 goto LBL_U;
2027         }
2028
2029         if ((res = pstm_init_size(pool, &A, sizeof(pstm_digit))) != PSTM_OKAY) {
2030                 goto LBL_V;
2031         }
2032
2033         if ((res = pstm_init_size(pool, &D, sizeof(pstm_digit))) != PSTM_OKAY) {
2034                 goto LBL_A;
2035         }
2036         pstm_set (&A, 1);
2037         pstm_set (&D, 1);
2038
2039         if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
2040                 goto LBL_D;
2041         }
2042         if ((res = pstm_init(pool, &C)) != PSTM_OKAY) {
2043                 goto LBL_B;
2044         }
2045
2046 top:
2047         /* 4.  while u is even do */
2048         while (pstm_iseven (&u) == 1) {
2049                 /* 4.1 u = u/2 */
2050                 if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
2051                         goto LBL_C;
2052                 }
2053
2054                 /* 4.2 if A or B is odd then */
2055                 if (pstm_isodd (&A) == 1 || pstm_isodd (&B) == 1) {
2056                         /* A = (A+y)/2, B = (B-x)/2 */
2057                         if ((res = pstm_add (&A, &y, &A)) != PSTM_OKAY) {
2058                                 goto LBL_C;
2059                         }
2060                         if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
2061                                 goto LBL_C;
2062                         }
2063                 }
2064                 /* A = A/2, B = B/2 */
2065                 if ((res = pstm_div_2 (&A, &A)) != PSTM_OKAY) {
2066                         goto LBL_C;
2067                 }
2068                 if ((res = pstm_div_2 (&B, &B)) != PSTM_OKAY) {
2069                         goto LBL_C;
2070                 }
2071         }
2072
2073         /* 5.  while v is even do */
2074         while (pstm_iseven (&v) == 1) {
2075                 /* 5.1 v = v/2 */
2076                 if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
2077                         goto LBL_C;
2078                 }
2079
2080                 /* 5.2 if C or D is odd then */
2081                 if (pstm_isodd (&C) == 1 || pstm_isodd (&D) == 1) {
2082                         /* C = (C+y)/2, D = (D-x)/2 */
2083                         if ((res = pstm_add (&C, &y, &C)) != PSTM_OKAY) {
2084                                 goto LBL_C;
2085                         }
2086                         if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
2087                                 goto LBL_C;
2088                         }
2089                 }
2090                 /* C = C/2, D = D/2 */
2091                 if ((res = pstm_div_2 (&C, &C)) != PSTM_OKAY) {
2092                         goto LBL_C;
2093                 }
2094                 if ((res = pstm_div_2 (&D, &D)) != PSTM_OKAY) {
2095                         goto LBL_C;
2096                 }
2097         }
2098
2099         /* 6.  if u >= v then */
2100         if (pstm_cmp (&u, &v) != PSTM_LT) {
2101                 /* u = u - v, A = A - C, B = B - D */
2102                 if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
2103                                 goto LBL_C;
2104                 }
2105                 if ((res = pstm_sub (&A, &C, &A)) != PSTM_OKAY) {
2106                                 goto LBL_C;
2107                 }
2108                 if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
2109                                 goto LBL_C;
2110                 }
2111         } else {
2112                 /* v - v - u, C = C - A, D = D - B */
2113                 if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
2114                                 goto LBL_C;
2115                 }
2116                 if ((res = pstm_sub (&C, &A, &C)) != PSTM_OKAY) {
2117                                 goto LBL_C;
2118                 }
2119                 if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
2120                                 goto LBL_C;
2121                 }
2122         }
2123
2124         /* if not zero goto step 4 */
2125         if (pstm_iszero (&u) == 0)
2126                 goto top;
2127
2128         /* now a = C, b = D, gcd == g*v */
2129
2130         /* if v != 1 then there is no inverse */
2131         if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
2132                 res = PS_FAILURE;
2133                 goto LBL_C;
2134         }
2135
2136         /* if its too low */
2137         while (pstm_cmp_d(&C, 0) == PSTM_LT) {
2138                 if ((res = pstm_add(&C, b, &C)) != PSTM_OKAY) {
2139                         goto LBL_C;
2140                 }
2141         }
2142
2143         /* too big */
2144         while (pstm_cmp_mag(&C, b) != PSTM_LT) {
2145                 if ((res = pstm_sub(&C, b, &C)) != PSTM_OKAY) {
2146                         goto LBL_C;
2147                 }
2148         }
2149
2150         /* C is now the inverse */
2151         if ((res = pstm_copy(&C, c)) != PSTM_OKAY) {
2152                 goto LBL_C;
2153         }
2154         res = PSTM_OKAY;
2155
2156 LBL_C: pstm_clear(&C);
2157 LBL_D: pstm_clear(&D);
2158 LBL_B: pstm_clear(&B);
2159 LBL_A: pstm_clear(&A);
2160 LBL_V: pstm_clear(&v);
2161 LBL_U: pstm_clear(&u);
2162 LBL_Y: pstm_clear(&y);
2163 LBL_X: pstm_clear(&x);
2164
2165         return res;
2166 }
2167
2168 /* c = 1/a (mod b) for odd b only */
2169 int32 pstm_invmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
2170 {
2171         pstm_int        x, y, u, v, B, D;
2172         int32           res;
2173         int             neg, sanity; //bbox: was uint16
2174
2175         /* 2. [modified] b must be odd   */
2176         if (pstm_iseven (b) == 1) {
2177                 return pstm_invmod_slow(pool, a,b,c);
2178         }
2179
2180         /* x == modulus, y == value to invert */
2181         if ((res = pstm_init_copy(pool, &x, b, 0)) != PSTM_OKAY) {
2182                 return res;
2183         }
2184
2185         if ((res = pstm_init_size(pool, &y, a->alloc)) != PSTM_OKAY) {
2186                 goto LBL_X;
2187         }
2188
2189         /* we need y = |a| */
2190         pstm_abs(a, &y);
2191
2192         /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2193         if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
2194                 goto LBL_Y;
2195         }
2196         if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
2197                 goto LBL_U;
2198         }
2199         if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
2200                 goto LBL_V;
2201         }
2202         if ((res = pstm_init(pool, &D)) != PSTM_OKAY) {
2203                 goto LBL_B;
2204         }
2205
2206         pstm_set (&D, 1);
2207
2208         sanity = 0;
2209 top:
2210         /* 4.  while u is even do */
2211         while (pstm_iseven (&u) == 1) {
2212                 /* 4.1 u = u/2 */
2213                 if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
2214                         goto LBL_D;
2215                 }
2216
2217                 /* 4.2 if B is odd then */
2218                 if (pstm_isodd (&B) == 1) {
2219                         if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
2220                                 goto LBL_D;
2221                         }
2222                 }
2223                 /* B = B/2 */
2224                 if ((res = pstm_div_2 (&B, &B)) !=  PSTM_OKAY) {
2225                         goto LBL_D;
2226                 }
2227         }
2228
2229         /* 5.  while v is even do */
2230         while (pstm_iseven (&v) == 1) {
2231                 /* 5.1 v = v/2 */
2232                 if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
2233                         goto LBL_D;
2234                 }
2235                 /* 5.2 if D is odd then */
2236                 if (pstm_isodd (&D) == 1) {
2237                         /* D = (D-x)/2 */
2238                         if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
2239                                 goto LBL_D;
2240                         }
2241                 }
2242                 /* D = D/2 */
2243                 if ((res = pstm_div_2 (&D, &D)) !=  PSTM_OKAY) {
2244                         goto LBL_D;
2245                 }
2246         }
2247
2248         /* 6.  if u >= v then */
2249         if (pstm_cmp (&u, &v) != PSTM_LT) {
2250                 /* u = u - v, B = B - D */
2251                 if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
2252                                 goto LBL_D;
2253                 }
2254                 if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
2255                                 goto LBL_D;
2256                 }
2257         } else {
2258                 /* v - v - u, D = D - B */
2259                 if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
2260                                 goto LBL_D;
2261                 }
2262                 if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
2263                                 goto LBL_D;
2264                 }
2265         }
2266
2267         /* if not zero goto step 4 */
2268         if (sanity++ > 1000) {
2269                 res = PS_LIMIT_FAIL;
2270                 goto LBL_D;
2271         }
2272         if (pstm_iszero (&u) == 0) {
2273                 goto top;
2274         }
2275
2276         /* now a = C, b = D, gcd == g*v */
2277
2278         /* if v != 1 then there is no inverse */
2279         if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
2280                 res = PS_FAILURE;
2281                 goto LBL_D;
2282         }
2283
2284         /* b is now the inverse */
2285         neg = a->sign;
2286         while (D.sign == PSTM_NEG) {
2287                 if ((res = pstm_add (&D, b, &D)) != PSTM_OKAY) {
2288                         goto LBL_D;
2289                 }
2290         }
2291         if ((res = pstm_copy (&D, c)) != PSTM_OKAY) {
2292                 goto LBL_D;
2293         }
2294         c->sign = neg;
2295         res = PSTM_OKAY;
2296
2297 LBL_D: pstm_clear(&D);
2298 LBL_B: pstm_clear(&B);
2299 LBL_V: pstm_clear(&v);
2300 LBL_U: pstm_clear(&u);
2301 LBL_Y: pstm_clear(&y);
2302 LBL_X: pstm_clear(&x);
2303         return res;
2304 }
2305 #endif //UNUSED
2306
2307 #endif /* !DISABLE_PSTM */
2308 /******************************************************************************/