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