library/bignum.c

   1 /*
   2  *      Multi-precision integer library
   3  *
   4  *      Based on XySSL: Copyright (C) 2006-2008  Christophe Devine
   5  *
   6  *      Copyright (C) 2009      Paul Bakker <polarssl_maintainer at polarssl dot org>
   7  *      Copyright (C) 2010 StackFoundry LLC <yann@stackfoundry.com>
   8  *      All rights reserved.
   9  *
  10  *      Redistribution and use in source and binary forms, with or without
  11  *      modification, are permitted provided that the following conditions
  12  *      are met:
  13  *
  14  *        * Redistributions of source code must retain the above copyright
  15  *              notice, this list of conditions and the following disclaimer.
  16  *        * Redistributions in binary form must reproduce the above copyright
  17  *              notice, this list of conditions and the following disclaimer in the
  18  *              documentation and/or other materials provided with the distribution.
  19  *        * Neither the names of Tropic SSL,
  20  *              PolarSSL or XySSL nor the names of its contributors
  21  *              may be used to endorse or promote products derived from this software
  22  *              without specific prior written permission.
  23  *
  24  *      THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  25  *      "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  26  *      LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  27  *      FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  28  *      OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  29  *      SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
  30  *      TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  31  *      PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  32  *      LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  33  *      NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  34  *      SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  35  */
  36 /*
  37  *      This MPI implementation is based on:
  38  *
  39  *      http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
  40  *      http://www.stillhq.com/extracted/gnupg-api/mpi/
  41  *      http://math.libtomcrypt.com/files/tommath.pdf
  42  */
  43
  44 #include "tropicssl/config.h"
  45
  46 #if defined(TROPICSSL_BIGNUM_C)
  47
  48 #include "tropicssl/bignum.h"
  49 #include "tropicssl/bn_mul.h"
  50
  51 #include <string.h>
  52 #include <stdlib.h>
  53 #include <stdarg.h>
  54
  55 #define ciL        ((int) sizeof(t_int))        /* chars in limb  */
  56 #define biL        (ciL << 3)   /* bits  in limb  */
  57 #define biH        (ciL << 2)   /* half limb size */
  58
  59 /*
  60  * Convert between bits/chars and number of limbs
  61  */
  62 #define BITS_TO_LIMBS(i)  (((i) + biL - 1) / biL)
  63 #define CHARS_TO_LIMBS(i) (((i) + ciL - 1) / ciL)
  64
  65 /*
  66  * Initialize one MPI
  67  */
  68 void mpi_init(mpi * X)
  69 {
  70         if (X == NULL) {
  71                 return;
  72         }
  73         X->s = 1;
  74         X->n = 0;
  75         X->p = NULL;
  76 }
  77
  78 /*
  79  * Unallocate one MPI
  80  */
  81 void mpi_free(mpi * X)
  82 {
  83         if (X == NULL) {
  84                 return;
  85         }
  86
  87         if (X->p != NULL) {
  88                 memset(X->p, 0, X->n * ciL);
  89                 free(X->p);
  90         }
  91
  92         X->s = 1;
  93         X->n = 0;
  94         X->p = NULL;
  95 }
  96
  97 /*
  98  * Enlarge to the specified number of limbs
  99  */
 100 int mpi_grow(mpi * X, int nblimbs)
 101 {
 102         t_int *p;
 103
 104         if (X->n < nblimbs) {
 105                 if ((p = (t_int *) malloc(nblimbs * ciL)) == NULL)
 106                         return (1);
 107
 108                 memset(p, 0, nblimbs * ciL);
 109
 110                 if (X->p != NULL) {
 111                         memcpy(p, X->p, X->n * ciL);
 112                         memset(X->p, 0, X->n * ciL);
 113                         free(X->p);
 114                 }
 115
 116                 X->n = nblimbs;
 117                 X->p = p;
 118         }
 119
 120         return (0);
 121 }
 122
 123 /*
 124  * Copy the contents of Y into X
 125  */
 126 int mpi_copy(mpi * X, const mpi * Y)
 127 {
 128         int ret, i;
 129
 130         if (X == Y)
 131                 return (0);
 132
 133         for (i = Y->n - 1; i > 0; i--)
 134                 if (Y->p[i] != 0)
 135                         break;
 136         i++;
 137
 138         X->s = Y->s;
 139
 140         MPI_CHK(mpi_grow(X, i));
 141
 142         memset(X->p, 0, X->n * ciL);
 143         memcpy(X->p, Y->p, i * ciL);
 144
 145 cleanup:
 146
 147         return (ret);
 148 }
 149
 150 /*
 151  * Swap the contents of X and Y
 152  */
 153 void mpi_swap(mpi * X, mpi * Y)
 154 {
 155         mpi T;
 156
 157         memcpy(&T, X, sizeof(mpi));
 158         memcpy(X, Y, sizeof(mpi));
 159         memcpy(Y, &T, sizeof(mpi));
 160 }
 161
 162 /*
 163  * Set value from integer
 164  */
 165 int mpi_lset(mpi * X, int z)
 166 {
 167         int ret;
 168
 169         MPI_CHK(mpi_grow(X, 1));
 170         memset(X->p, 0, X->n * ciL);
 171
 172         X->p[0] = (z < 0) ? -z : z;
 173         X->s = (z < 0) ? -1 : 1;
 174
 175 cleanup:
 176
 177         return (ret);
 178 }
 179
 180 /*
 181  * Return the number of least significant bits
 182  */
 183 int mpi_lsb(const mpi * X)
 184 {
 185         int i, j, count = 0;
 186
 187         for (i = 0; i < X->n; i++)
 188                 for (j = 0; j < (int)biL; j++, count++)
 189                         if (((X->p[i] >> j) & 1) != 0)
 190                                 return (count);
 191
 192         return (0);
 193 }
 194
 195 /*
 196  * Return the number of most significant bits
 197  */
 198 int mpi_msb(const mpi * X)
 199 {
 200         int i, j;
 201
 202         for (i = X->n - 1; i > 0; i--)
 203                 if (X->p[i] != 0)
 204                         break;
 205
 206         for (j = biL - 1; j >= 0; j--)
 207                 if (((X->p[i] >> j) & 1) != 0)
 208                         break;
 209
 210         return ((i * biL) + j + 1);
 211 }
 212
 213 /*
 214  * Return the total size in bytes
 215  */
 216 int mpi_size(const mpi * X)
 217 {
 218         return ((mpi_msb(X) + 7) >> 3);
 219 }
 220
 221 /*
 222  * Convert an ASCII character to digit value
 223  */
 224 static int mpi_get_digit(t_int * d, int radix, char c)
 225 {
 226         *d = 255;
 227
 228         if (c >= 0x30 && c <= 0x39)
 229                 *d = c - 0x30;
 230         if (c >= 0x41 && c <= 0x46)
 231                 *d = c - 0x37;
 232         if (c >= 0x61 && c <= 0x66)
 233                 *d = c - 0x57;
 234
 235         if (*d >= (t_int) radix)
 236                 return (TROPICSSL_ERR_MPI_INVALID_CHARACTER);
 237
 238         return (0);
 239 }
 240
 241 /*
 242  * Import from an ASCII string
 243  */
 244 int mpi_read_string(mpi * X, int radix, const char *s)
 245 {
 246         int ret, i, j, n;
 247         t_int d;
 248         mpi T;
 249
 250         if (radix < 2 || radix > 16)
 251                 return (TROPICSSL_ERR_MPI_BAD_INPUT_DATA);
 252
 253         mpi_init(&T);
 254
 255         if (radix == 16) {
 256                 n = BITS_TO_LIMBS(strlen(s) << 2);
 257
 258                 MPI_CHK(mpi_grow(X, n));
 259                 MPI_CHK(mpi_lset(X, 0));
 260
 261                 for (i = strlen(s) - 1, j = 0; i >= 0; i--, j++) {
 262                         if (i == 0 && s[i] == '-') {
 263                                 X->s = -1;
 264                                 break;
 265                         }
 266
 267                         MPI_CHK(mpi_get_digit(&d, radix, s[i]));
 268                         X->p[j / (2 * ciL)] |= d << ((j % (2 * ciL)) << 2);
 269                 }
 270         } else {
 271                 MPI_CHK(mpi_lset(X, 0));
 272
 273                 for (i = 0; i < (int)strlen(s); i++) {
 274                         if (i == 0 && s[i] == '-') {
 275                                 X->s = -1;
 276                                 continue;
 277                         }
 278
 279                         MPI_CHK(mpi_get_digit(&d, radix, s[i]));
 280                         MPI_CHK(mpi_mul_int(&T, X, radix));
 281                         MPI_CHK(mpi_add_int(X, &T, d));
 282                 }
 283         }
 284
 285 cleanup:
 286
 287         mpi_free(&T);
 288
 289         return (ret);
 290 }
 291
 292 /*
 293  * Helper to write the digits high-order first
 294  */
 295 static int mpi_write_hlp(mpi * X, int radix, char **p)
 296 {
 297         int ret;
 298         t_int r;
 299
 300         if (radix < 2 || radix > 16)
 301                 return (TROPICSSL_ERR_MPI_BAD_INPUT_DATA);
 302
 303         MPI_CHK(mpi_mod_int(&r, X, radix));
 304         MPI_CHK(mpi_div_int(X, NULL, X, radix));
 305
 306         if (mpi_cmp_int(X, 0) != 0)
 307                 MPI_CHK(mpi_write_hlp(X, radix, p));
 308
 309         if (r < 10)
 310                 *(*p)++ = (char)(r + 0x30);
 311         else
 312                 *(*p)++ = (char)(r + 0x37);
 313
 314 cleanup:
 315
 316         return (ret);
 317 }
 318
 319 /*
 320  * Export into an ASCII string
 321  */
 322 int mpi_write_string(const mpi * X, int radix, char *s, int *slen)
 323 {
 324         int ret = 0, n;
 325         char *p;
 326         mpi T;
 327
 328         if (radix < 2 || radix > 16)
 329                 return (TROPICSSL_ERR_MPI_BAD_INPUT_DATA);
 330
 331         n = mpi_msb(X);
 332         if (radix >= 4)
 333                 n >>= 1;
 334         if (radix >= 16)
 335                 n >>= 1;
 336         n += 3;
 337
 338         if (*slen < n) {
 339                 *slen = n;
 340                 return (TROPICSSL_ERR_MPI_BUFFER_TOO_SMALL);
 341         }
 342
 343         p = s;
 344         mpi_init(&T);
 345
 346         if (X->s == -1)
 347                 *p++ = '-';
 348
 349         if (radix == 16) {
 350                 int c, i, j, k;
 351
 352                 for (i = X->n - 1, k = 0; i >= 0; i--) {
 353                         for (j = ciL - 1; j >= 0; j--) {
 354                                 c = (X->p[i] >> (j << 3)) & 0xFF;
 355
 356                                 if (c == 0 && k == 0 && (i + j) != 0)
 357                                         continue;
 358
 359                                 p += sprintf(p, "%02X", c);
 360                                 k = 1;
 361                         }
 362                 }
 363         } else {
 364                 MPI_CHK(mpi_copy(&T, X));
 365                 MPI_CHK(mpi_write_hlp(&T, radix, &p));
 366         }
 367
 368         *p++ = '\0';
 369         *slen = p - s;
 370
 371 cleanup:
 372
 373         mpi_free(&T);
 374
 375         return (ret);
 376 }
 377
 378 /*
 379  * Read X from an opened file
 380  */
 381 int mpi_read_file(mpi * X, int radix, FILE * fin)
 382 {
 383         t_int d;
 384         int slen;
 385         char *p;
 386         char s[1024];
 387
 388         memset(s, 0, sizeof(s));
 389         if (fgets(s, sizeof(s) - 1, fin) == NULL)
 390                 return (TROPICSSL_ERR_MPI_FILE_IO_ERROR);
 391
 392         slen = strlen(s);
 393         if (s[slen - 1] == '\n') {
 394                 slen--;
 395                 s[slen] = '\0';
 396         }
 397         if (s[slen - 1] == '\r') {
 398                 slen--;
 399                 s[slen] = '\0';
 400         }
 401
 402         p = s + slen;
 403         while (--p >= s)
 404                 if (mpi_get_digit(&d, radix, *p) != 0)
 405                         break;
 406
 407         return (mpi_read_string(X, radix, p + 1));
 408 }
 409
 410 /*
 411  * Write X into an opened file (or stdout if fout == NULL)
 412  */
 413 int mpi_write_file(const char *p, const mpi * X, int radix, FILE * fout)
 414 {
 415         int n, ret;
 416         size_t slen;
 417         size_t plen;
 418         char s[1024];
 419
 420         n = sizeof(s);
 421         memset(s, 0, n);
 422         n -= 2;
 423
 424         MPI_CHK(mpi_write_string(X, radix, s, (int *)&n));
 425
 426         if (p == NULL)
 427                 p = "";
 428
 429         plen = strlen(p);
 430         slen = strlen(s);
 431         s[slen++] = '\r';
 432         s[slen++] = '\n';
 433
 434         if (fout != NULL) {
 435                 if (fwrite(p, 1, plen, fout) != plen ||
 436                     fwrite(s, 1, slen, fout) != slen)
 437                         return (TROPICSSL_ERR_MPI_FILE_IO_ERROR);
 438         } else
 439                 printf("%s%s", p, s);
 440
 441 cleanup:
 442
 443         return (ret);
 444 }
 445
 446 /*
 447  * Import X from unsigned binary data, big endian
 448  */
 449 int mpi_read_binary(mpi * X, const unsigned char *buf, int buflen)
 450 {
 451         int ret, i, j, n;
 452
 453         for (n = 0; n < buflen; n++)
 454                 if (buf[n] != 0)
 455                         break;
 456
 457         MPI_CHK(mpi_grow(X, CHARS_TO_LIMBS(buflen - n)));
 458         MPI_CHK(mpi_lset(X, 0));
 459
 460         for (i = buflen - 1, j = 0; i >= n; i--, j++)
 461                 X->p[j / ciL] |= ((t_int) buf[i]) << ((j % ciL) << 3);
 462
 463 cleanup:
 464
 465         return (ret);
 466 }
 467
 468 /*
 469  * Export X into unsigned binary data, big endian
 470  */
 471 int mpi_write_binary(const mpi * X, unsigned char *buf, int buflen)
 472 {
 473         int i, j, n;
 474
 475         n = mpi_size(X);
 476
 477         if (buflen < n)
 478                 return (TROPICSSL_ERR_MPI_BUFFER_TOO_SMALL);
 479
 480         memset(buf, 0, buflen);
 481
 482         for (i = buflen - 1, j = 0; n > 0; i--, j++, n--)
 483                 buf[i] = (unsigned char)(X->p[j / ciL] >> ((j % ciL) << 3));
 484
 485         return (0);
 486 }
 487
 488 /*
 489  * Left-shift: X <<= count
 490  */
 491 int mpi_shift_l(mpi * X, int count)
 492 {
 493         int ret, i, v0, t1;
 494         t_int r0 = 0, r1;
 495
 496         v0 = count / (biL);
 497         t1 = count & (biL - 1);
 498
 499         i = mpi_msb(X) + count;
 500
 501         if (X->n * (int)biL < i)
 502                 MPI_CHK(mpi_grow(X, BITS_TO_LIMBS(i)));
 503
 504         ret = 0;
 505
 506         /*
 507          * shift by count / limb_size
 508          */
 509         if (v0 > 0) {
 510                 for (i = X->n - 1; i >= v0; i--)
 511                         X->p[i] = X->p[i - v0];
 512
 513                 for (; i >= 0; i--)
 514                         X->p[i] = 0;
 515         }
 516
 517         /*
 518          * shift by count % limb_size
 519          */
 520         if (t1 > 0) {
 521                 for (i = v0; i < X->n; i++) {
 522                         r1 = X->p[i] >> (biL - t1);
 523                         X->p[i] <<= t1;
 524                         X->p[i] |= r0;
 525                         r0 = r1;
 526                 }
 527         }
 528
 529 cleanup:
 530
 531         return (ret);
 532 }
 533
 534 /*
 535  * Right-shift: X >>= count
 536  */
 537 int mpi_shift_r(mpi * X, int count)
 538 {
 539         int i, v0, v1;
 540         t_int r0 = 0, r1;
 541
 542         v0 = count / biL;
 543         v1 = count & (biL - 1);
 544
 545         /*
 546          * shift by count / limb_size
 547          */
 548         if (v0 > 0) {
 549                 for (i = 0; i < X->n - v0; i++)
 550                         X->p[i] = X->p[i + v0];
 551
 552                 for (; i < X->n; i++)
 553                         X->p[i] = 0;
 554         }
 555
 556         /*
 557          * shift by count % limb_size
 558          */
 559         if (v1 > 0) {
 560                 for (i = X->n - 1; i >= 0; i--) {
 561                         r1 = X->p[i] << (biL - v1);
 562                         X->p[i] >>= v1;
 563                         X->p[i] |= r0;
 564                         r0 = r1;
 565                 }
 566         }
 567
 568         return (0);
 569 }
 570
 571 /*
 572  * Compare unsigned values
 573  */
 574 int mpi_cmp_abs(const mpi * X, const mpi * Y)
 575 {
 576         int i, j;
 577
 578         for (i = X->n - 1; i >= 0; i--)
 579                 if (X->p[i] != 0)
 580                         break;
 581
 582         for (j = Y->n - 1; j >= 0; j--)
 583                 if (Y->p[j] != 0)
 584                         break;
 585
 586         if (i < 0 && j < 0)
 587                 return (0);
 588
 589         if (i > j)
 590                 return (1);
 591         if (j > i)
 592                 return (-1);
 593
 594         for (; i >= 0; i--) {
 595                 if (X->p[i] > Y->p[i])
 596                         return (1);
 597                 if (X->p[i] < Y->p[i])
 598                         return (-1);
 599         }
 600
 601         return (0);
 602 }
 603
 604 /*
 605  * Compare signed values
 606  */
 607 int mpi_cmp_mpi(const mpi * X, const mpi * Y)
 608 {
 609         int i, j;
 610
 611         for (i = X->n - 1; i >= 0; i--)
 612                 if (X->p[i] != 0)
 613                         break;
 614
 615         for (j = Y->n - 1; j >= 0; j--)
 616                 if (Y->p[j] != 0)
 617                         break;
 618
 619         if (i < 0 && j < 0)
 620                 return (0);
 621
 622         if (i > j)
 623                 return (X->s);
 624         if (j > i)
 625                 return (-X->s);
 626
 627         if (X->s > 0 && Y->s < 0)
 628                 return (1);
 629         if (Y->s > 0 && X->s < 0)
 630                 return (-1);
 631
 632         for (; i >= 0; i--) {
 633                 if (X->p[i] > Y->p[i])
 634                         return (X->s);
 635                 if (X->p[i] < Y->p[i])
 636                         return (-X->s);
 637         }
 638
 639         return (0);
 640 }
 641
 642 /*
 643  * Compare signed values
 644  */
 645 int mpi_cmp_int(const mpi * X, int z)
 646 {
 647         mpi Y;
 648         t_int p[1];
 649
 650         *p = (z < 0) ? -z : z;
 651         Y.s = (z < 0) ? -1 : 1;
 652         Y.n = 1;
 653         Y.p = p;
 654
 655         return (mpi_cmp_mpi(X, &Y));
 656 }
 657
 658 /*
 659  * Unsigned addition: X = |A| + |B|      (HAC 14.7)
 660  */
 661 int mpi_add_abs(mpi * X, const mpi * A, const mpi * B)
 662 {
 663         int ret, i, j;
 664         t_int *o, *p, c;
 665
 666         if (X == B) {
 667                 const mpi *T = A;
 668                 A = X;
 669                 B = T;
 670         }
 671
 672         if (X != A)
 673                 MPI_CHK(mpi_copy(X, A));
 674
 675         for (j = B->n - 1; j >= 0; j--)
 676                 if (B->p[j] != 0)
 677                         break;
 678
 679         MPI_CHK(mpi_grow(X, j + 1));
 680
 681         o = B->p;
 682         p = X->p;
 683         c = 0;
 684
 685         for (i = 0; i <= j; i++, o++, p++) {
 686                 *p += c;
 687                 c = (*p < c);
 688                 *p += *o;
 689                 c += (*p < *o);
 690         }
 691
 692         while (c != 0) {
 693                 if (i >= X->n) {
 694                         MPI_CHK(mpi_grow(X, i + 1));
 695                         p = X->p + i;
 696                 }
 697
 698                 *p += c;
 699                 c = (*p < c);
 700                 i++;
 701         }
 702
 703 cleanup:
 704
 705         return (ret);
 706 }
 707
 708 /*
 709  * Helper for mpi substraction
 710  */
 711 static void mpi_sub_hlp(int n, t_int * s, t_int * d)
 712 {
 713         int i;
 714         t_int c, z;
 715
 716         for (i = c = 0; i < n; i++, s++, d++) {
 717                 z = (*d < c);
 718                 *d -= c;
 719                 c = (*d < *s) + z;
 720                 *d -= *s;
 721         }
 722
 723         while (c != 0) {
 724                 z = (*d < c);
 725                 *d -= c;
 726                 c = z;
 727                 i++;
 728                 d++;
 729         }
 730 }
 731
 732 /*
 733  * Unsigned substraction: X = |A| - |B|  (HAC 14.9)
 734  */
 735 int mpi_sub_abs(mpi * X, const mpi * A, const mpi * B)
 736 {
 737         mpi TB;
 738         int ret, n;
 739
 740         if (mpi_cmp_abs(A, B) < 0)
 741                 return (TROPICSSL_ERR_MPI_NEGATIVE_VALUE);
 742
 743         mpi_init(&TB);
 744
 745         if (X == B) {
 746                 MPI_CHK(mpi_copy(&TB, B));
 747                 B = &TB;
 748         }
 749
 750         if (X != A)
 751                 MPI_CHK(mpi_copy(X, A));
 752
 753         ret = 0;
 754
 755         for (n = B->n - 1; n >= 0; n--)
 756                 if (B->p[n] != 0)
 757                         break;
 758
 759         mpi_sub_hlp(n + 1, B->p, X->p);
 760
 761 cleanup:
 762
 763         mpi_free(&TB);
 764
 765         return (ret);
 766 }
 767
 768 /*
 769  * Signed addition: X = A + B
 770  */
 771 int mpi_add_mpi(mpi * X, const mpi * A, const mpi * B)
 772 {
 773         int ret, s = A->s;
 774
 775         if (A->s * B->s < 0) {
 776                 if (mpi_cmp_abs(A, B) >= 0) {
 777                         MPI_CHK(mpi_sub_abs(X, A, B));
 778                         X->s = s;
 779                 } else {
 780                         MPI_CHK(mpi_sub_abs(X, B, A));
 781                         X->s = -s;
 782                 }
 783         } else {
 784                 MPI_CHK(mpi_add_abs(X, A, B));
 785                 X->s = s;
 786         }
 787
 788 cleanup:
 789
 790         return (ret);
 791 }
 792
 793 /*
 794  * Signed substraction: X = A - B
 795  */
 796 int mpi_sub_mpi(mpi * X, const mpi * A, const mpi * B)
 797 {
 798         int ret, s = A->s;
 799
 800         if (A->s * B->s > 0) {
 801                 if (mpi_cmp_abs(A, B) >= 0) {
 802                         MPI_CHK(mpi_sub_abs(X, A, B));
 803                         X->s = s;
 804                 } else {
 805                         MPI_CHK(mpi_sub_abs(X, B, A));
 806                         X->s = -s;
 807                 }
 808         } else {
 809                 MPI_CHK(mpi_add_abs(X, A, B));
 810                 X->s = s;
 811         }
 812
 813 cleanup:
 814
 815         return (ret);
 816 }
 817
 818 /*
 819  * Signed addition: X = A + b
 820  */
 821 int mpi_add_int(mpi * X, const mpi * A, int b)
 822 {
 823         mpi _B;
 824         t_int p[1];
 825
 826         p[0] = (b < 0) ? -b : b;
 827         _B.s = (b < 0) ? -1 : 1;
 828         _B.n = 1;
 829         _B.p = p;
 830
 831         return (mpi_add_mpi(X, A, &_B));
 832 }
 833
 834 /*
 835  * Signed substraction: X = A - b
 836  */
 837 int mpi_sub_int(mpi * X, const mpi * A, int b)
 838 {
 839         mpi _B;
 840         t_int p[1];
 841
 842         p[0] = (b < 0) ? -b : b;
 843         _B.s = (b < 0) ? -1 : 1;
 844         _B.n = 1;
 845         _B.p = p;
 846
 847         return (mpi_sub_mpi(X, A, &_B));
 848 }
 849
 850 /*
 851  * Helper for mpi multiplication
 852  */
 853 static void mpi_mul_hlp(int i, t_int * s, t_int * d, t_int b)
 854 {
 855         t_int c = 0, t = 0;
 856
 857 #if defined(MULADDC_HUIT)
 858         for (; i >= 8; i -= 8) {
 859         MULADDC_INIT MULADDC_HUIT MULADDC_STOP} for (; i > 0; i--) {
 860         MULADDC_INIT MULADDC_CORE MULADDC_STOP}
 861 #else
 862         for (; i >= 16; i -= 16) {
 863         MULADDC_INIT
 864                     MULADDC_CORE MULADDC_CORE
 865                     MULADDC_CORE MULADDC_CORE
 866                     MULADDC_CORE MULADDC_CORE
 867                     MULADDC_CORE MULADDC_CORE
 868                     MULADDC_CORE MULADDC_CORE
 869                     MULADDC_CORE MULADDC_CORE
 870                     MULADDC_CORE MULADDC_CORE
 871                     MULADDC_CORE MULADDC_CORE MULADDC_STOP}
 872
 873         for (; i >= 8; i -= 8) {
 874         MULADDC_INIT
 875                     MULADDC_CORE MULADDC_CORE
 876                     MULADDC_CORE MULADDC_CORE
 877                     MULADDC_CORE MULADDC_CORE
 878                     MULADDC_CORE MULADDC_CORE MULADDC_STOP}
 879
 880         for (; i > 0; i--) {
 881         MULADDC_INIT MULADDC_CORE MULADDC_STOP}
 882 #endif
 883         t++;
 884
 885         do {
 886                 *d += c;
 887                 c = (*d < c);
 888                 d++;
 889         } while (c != 0);
 890 }
 891
 892 /*
 893  * Baseline multiplication: X = A * B  (HAC 14.12)
 894  */
 895 int mpi_mul_mpi(mpi * X, const mpi * A, const mpi * B)
 896 {
 897         int ret, i, j;
 898         mpi TA, TB;
 899
 900         mpi_init(&TA); mpi_init(&TB);
 901
 902         if (X == A) {
 903                 MPI_CHK(mpi_copy(&TA, A));
 904                 A = &TA;
 905         }
 906         if (X == B) {
 907                 MPI_CHK(mpi_copy(&TB, B));
 908                 B = &TB;
 909         }
 910
 911         for (i = A->n - 1; i >= 0; i--)
 912                 if (A->p[i] != 0)
 913                         break;
 914
 915         for (j = B->n - 1; j >= 0; j--)
 916                 if (B->p[j] != 0)
 917                         break;
 918
 919         MPI_CHK(mpi_grow(X, i + j + 2));
 920         MPI_CHK(mpi_lset(X, 0));
 921
 922         for (i++; j >= 0; j--)
 923                 mpi_mul_hlp(i, A->p, X->p + j, B->p[j]);
 924
 925         X->s = A->s * B->s;
 926
 927 cleanup:
 928
 929         mpi_free(&TB); mpi_free(&TA);
 930
 931         return (ret);
 932 }
 933
 934 /*
 935  * Baseline multiplication: X = A * b
 936  */
 937 int mpi_mul_int(mpi * X, const mpi * A, t_int b)
 938 {
 939         mpi _B;
 940         t_int p[1];
 941
 942         _B.s = 1;
 943         _B.n = 1;
 944         _B.p = p;
 945         p[0] = b;
 946
 947         return (mpi_mul_mpi(X, A, &_B));
 948 }
 949
 950 /*
 951  * Division by mpi: A = Q * B + R  (HAC 14.20)
 952  */
 953 int mpi_div_mpi(mpi * Q, mpi * R, const mpi * A, const mpi * B)
 954 {
 955         int ret, i, n, t, k;
 956         mpi X, Y, Z, T1, T2;
 957
 958         if (mpi_cmp_int(B, 0) == 0)
 959                 return (TROPICSSL_ERR_MPI_DIVISION_BY_ZERO);
 960
 961         mpi_init(&X); mpi_init(&Y); mpi_init(&Z);
 962         mpi_init(&T1); mpi_init(&T2);
 963
 964         if (mpi_cmp_abs(A, B) < 0) {
 965                 if (Q != NULL)
 966                         MPI_CHK(mpi_lset(Q, 0));
 967                 if (R != NULL)
 968                         MPI_CHK(mpi_copy(R, A));
 969                 return (0);
 970         }
 971
 972         MPI_CHK(mpi_copy(&X, A));
 973         MPI_CHK(mpi_copy(&Y, B));
 974         X.s = Y.s = 1;
 975
 976         MPI_CHK(mpi_grow(&Z, A->n + 2));
 977         MPI_CHK(mpi_lset(&Z, 0));
 978         MPI_CHK(mpi_grow(&T1, 2));
 979         MPI_CHK(mpi_grow(&T2, 3));
 980
 981         k = mpi_msb(&Y) % biL;
 982         if (k < (int)biL - 1) {
 983                 k = biL - 1 - k;
 984                 MPI_CHK(mpi_shift_l(&X, k));
 985                 MPI_CHK(mpi_shift_l(&Y, k));
 986         } else
 987                 k = 0;
 988
 989         n = X.n - 1;
 990         t = Y.n - 1;
 991         mpi_shift_l(&Y, biL * (n - t));
 992
 993         while (mpi_cmp_mpi(&X, &Y) >= 0) {
 994                 Z.p[n - t]++;
 995                 mpi_sub_mpi(&X, &X, &Y);
 996         }
 997         mpi_shift_r(&Y, biL * (n - t));
 998
 999         for (i = n; i > t; i--) {
1000                 if (X.p[i] >= Y.p[t])
1001                         Z.p[i - t - 1] = ~0;
1002                 else {
1003 #if defined(TROPICSSL_HAVE_LONGLONG)
1004                         t_dbl r;
1005
1006                         r = (t_dbl) X.p[i] << biL;
1007                         r |= (t_dbl) X.p[i - 1];
1008                         r /= Y.p[t];
1009                         if (r > ((t_dbl) 1 << biL) - 1)
1010                                 r = ((t_dbl) 1 << biL) - 1;
1011
1012                         Z.p[i - t - 1] = (t_int) r;
1013 #else
1014                         /*
1015                          * __udiv_qrnnd_c, from gmp/longlong.h
1016                          */
1017                         t_int q0, q1, r0, r1;
1018                         t_int d0, d1, d, m;
1019
1020                         d = Y.p[t];
1021                         d0 = (d << biH) >> biH;
1022                         d1 = (d >> biH);
1023
1024                         q1 = X.p[i] / d1;
1025                         r1 = X.p[i] - d1 * q1;
1026                         r1 <<= biH;
1027                         r1 |= (X.p[i - 1] >> biH);
1028
1029                         m = q1 * d0;
1030                         if (r1 < m) {
1031                                 q1--, r1 += d;
1032                                 while (r1 >= d && r1 < m)
1033                                         q1--, r1 += d;
1034                         }
1035                         r1 -= m;
1036
1037                         q0 = r1 / d1;
1038                         r0 = r1 - d1 * q0;
1039                         r0 <<= biH;
1040                         r0 |= (X.p[i - 1] << biH) >> biH;
1041
1042                         m = q0 * d0;
1043                         if (r0 < m) {
1044                                 q0--, r0 += d;
1045                                 while (r0 >= d && r0 < m)
1046                                         q0--, r0 += d;
1047                         }
1048                         r0 -= m;
1049
1050                         Z.p[i - t - 1] = (q1 << biH) | q0;
1051 #endif
1052                 }
1053
1054                 Z.p[i - t - 1]++;
1055                 do {
1056                         Z.p[i - t - 1]--;
1057
1058                         MPI_CHK(mpi_lset(&T1, 0));
1059                         T1.p[0] = (t < 1) ? 0 : Y.p[t - 1];
1060                         T1.p[1] = Y.p[t];
1061                         MPI_CHK(mpi_mul_int(&T1, &T1, Z.p[i - t - 1]));
1062
1063                         MPI_CHK(mpi_lset(&T2, 0));
1064                         T2.p[0] = (i < 2) ? 0 : X.p[i - 2];
1065                         T2.p[1] = (i < 1) ? 0 : X.p[i - 1];
1066                         T2.p[2] = X.p[i];
1067                 } while (mpi_cmp_mpi(&T1, &T2) > 0);
1068
1069                 MPI_CHK(mpi_mul_int(&T1, &Y, Z.p[i - t - 1]));
1070                 MPI_CHK(mpi_shift_l(&T1, biL * (i - t - 1)));
1071                 MPI_CHK(mpi_sub_mpi(&X, &X, &T1));
1072
1073                 if (mpi_cmp_int(&X, 0) < 0) {
1074                         MPI_CHK(mpi_copy(&T1, &Y));
1075                         MPI_CHK(mpi_shift_l(&T1, biL * (i - t - 1)));
1076                         MPI_CHK(mpi_add_mpi(&X, &X, &T1));
1077                         Z.p[i - t - 1]--;
1078                 }
1079         }
1080
1081         if (Q != NULL) {
1082                 mpi_copy(Q, &Z);
1083                 Q->s = A->s * B->s;
1084         }
1085
1086         if (R != NULL) {
1087                 mpi_shift_r(&X, k);
1088                 mpi_copy(R, &X);
1089
1090                 R->s = A->s;
1091                 if (mpi_cmp_int(R, 0) == 0)
1092                         R->s = 1;
1093         }
1094
1095 cleanup:
1096
1097         mpi_free(&X); mpi_free(&Y); mpi_free(&Z);
1098         mpi_free(&T1); mpi_free(&T2);
1099
1100         return (ret);
1101 }
1102
1103 /*
1104  * Division by int: A = Q * b + R
1105  *
1106  * Returns 0 if successful
1107  *                 1 if memory allocation failed
1108  *                 TROPICSSL_ERR_MPI_DIVISION_BY_ZERO if b == 0
1109  */
1110 int mpi_div_int(mpi * Q, mpi * R, const mpi * A, int b)
1111 {
1112         mpi _B;
1113         t_int p[1];
1114
1115         p[0] = (b < 0) ? -b : b;
1116         _B.s = (b < 0) ? -1 : 1;
1117         _B.n = 1;
1118         _B.p = p;
1119
1120         return (mpi_div_mpi(Q, R, A, &_B));
1121 }
1122
1123 /*
1124  * Modulo: R = A mod B
1125  */
1126 int mpi_mod_mpi(mpi * R, const mpi * A, const mpi * B)
1127 {
1128         int ret;
1129
1130         MPI_CHK(mpi_div_mpi(NULL, R, A, B));
1131
1132         while (mpi_cmp_int(R, 0) < 0)
1133                 MPI_CHK(mpi_add_mpi(R, R, B));
1134
1135         while (mpi_cmp_mpi(R, B) >= 0)
1136                 MPI_CHK(mpi_sub_mpi(R, R, B));
1137
1138 cleanup:
1139
1140         return (ret);
1141 }
1142
1143 /*
1144  * Modulo: r = A mod b
1145  */
1146 int mpi_mod_int(t_int * r, const mpi * A, int b)
1147 {
1148         int i;
1149         t_int x, y, z;
1150
1151         if (b == 0)
1152                 return (TROPICSSL_ERR_MPI_DIVISION_BY_ZERO);
1153
1154         if (b < 0)
1155                 b = -b;
1156
1157         /*
1158          * handle trivial cases
1159          */
1160         if (b == 1) {
1161                 *r = 0;
1162                 return (0);
1163         }
1164
1165         if (b == 2) {
1166                 *r = A->p[0] & 1;
1167                 return (0);
1168         }
1169
1170         /*
1171          * general case
1172          */
1173         for (i = A->n - 1, y = 0; i >= 0; i--) {
1174                 x = A->p[i];
1175                 y = (y << biH) | (x >> biH);
1176                 z = y / b;
1177                 y -= z * b;
1178
1179                 x <<= biH;
1180                 y = (y << biH) | (x >> biH);
1181                 z = y / b;
1182                 y -= z * b;
1183         }
1184
1185         *r = y;
1186
1187         return (0);
1188 }
1189
1190 /*
1191  * Fast Montgomery initialization (thanks to Tom St Denis)
1192  */
1193 static void mpi_montg_init(t_int * mm, const mpi * N)
1194 {
1195         t_int x, m0 = N->p[0];
1196
1197         x = m0;
1198         x += ((m0 + 2) & 4) << 1;
1199         x *= (2 - (m0 * x));
1200
1201         if (biL >= 16)
1202                 x *= (2 - (m0 * x));
1203         if (biL >= 32)
1204                 x *= (2 - (m0 * x));
1205         if (biL >= 64)
1206                 x *= (2 - (m0 * x));
1207
1208         *mm = ~x + 1;
1209 }
1210
1211 /*
1212  * Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
1213  */
1214 static void mpi_montmul(mpi * A, const mpi * B, const mpi * N, t_int mm, mpi * T)
1215 {
1216         int i, n, m;
1217         t_int u0, u1, *d;
1218
1219         memset(T->p, 0, T->n * ciL);
1220
1221         d = T->p;
1222         n = N->n;
1223         m = (B->n < n) ? B->n : n;
1224
1225         for (i = 0; i < n; i++) {
1226                 /*
1227                  * T = (T + u0*B + u1*N) / 2^biL
1228                  */
1229                 u0 = A->p[i];
1230                 u1 = (d[0] + u0 * B->p[0]) * mm;
1231
1232                 mpi_mul_hlp(m, B->p, d, u0);
1233                 mpi_mul_hlp(n, N->p, d, u1);
1234
1235                 *d++ = u0;
1236                 d[n + 1] = 0;
1237         }
1238
1239         memcpy(A->p, d, (n + 1) * ciL);
1240
1241         if (mpi_cmp_abs(A, N) >= 0)
1242                 mpi_sub_hlp(n, N->p, A->p);
1243         else
1244                 /* prevent timing attacks */
1245                 mpi_sub_hlp(n, A->p, T->p);
1246 }
1247
1248 /*
1249  * Montgomery reduction: A = A * R^-1 mod N
1250  */
1251 static void mpi_montred(mpi * A, const mpi * N, t_int mm, mpi * T)
1252 {
1253         t_int z = 1;
1254         mpi U;
1255
1256         U.n = U.s = z;
1257         U.p = &z;
1258
1259         mpi_montmul(A, &U, N, mm, T);
1260 }
1261
1262 /*
1263  * Sliding-window exponentiation: X = A^E mod N  (HAC 14.85)
1264  */
1265 int mpi_exp_mod(mpi * X, const mpi * A, const mpi * E, const mpi * N, mpi * _RR)
1266 {
1267         int ret, i, j, wsize, wbits;
1268         int bufsize, nblimbs, nbits;
1269         t_int ei, mm, state;
1270         mpi RR, T, W[64];
1271
1272         if (mpi_cmp_int(N, 0) < 0 || (N->p[0] & 1) == 0)
1273                 return (TROPICSSL_ERR_MPI_BAD_INPUT_DATA);
1274
1275         /*
1276          * Init temps and window size
1277          */
1278         mpi_montg_init(&mm, N);
1279         mpi_init(&RR); mpi_init(&T);
1280         memset(W, 0, sizeof(W));
1281
1282         i = mpi_msb(E);
1283
1284         wsize = (i > 671) ? 6 : (i > 239) ? 5 : (i > 79) ? 4 : (i > 23) ? 3 : 1;
1285
1286         j = N->n + 1;
1287         MPI_CHK(mpi_grow(X, j));
1288         MPI_CHK(mpi_grow(&W[1], j));
1289         MPI_CHK(mpi_grow(&T, j * 2));
1290
1291         /*
1292          * If 1st call, pre-compute R^2 mod N
1293          */
1294         if (_RR == NULL || _RR->p == NULL) {
1295                 MPI_CHK(mpi_lset(&RR, 1));
1296                 MPI_CHK(mpi_shift_l(&RR, N->n * 2 * biL));
1297                 MPI_CHK(mpi_mod_mpi(&RR, &RR, N));
1298
1299                 if (_RR != NULL)
1300                         memcpy(_RR, &RR, sizeof(mpi));
1301         } else
1302                 memcpy(&RR, _RR, sizeof(mpi));
1303
1304         /*
1305          * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1306          */
1307         if (mpi_cmp_mpi(A, N) >= 0)
1308                 mpi_mod_mpi(&W[1], A, N);
1309         else
1310                 mpi_copy(&W[1], A);
1311
1312         mpi_montmul(&W[1], &RR, N, mm, &T);
1313
1314         /*
1315          * X = R^2 * R^-1 mod N = R mod N
1316          */
1317         MPI_CHK(mpi_copy(X, &RR));
1318         mpi_montred(X, N, mm, &T);
1319
1320         if (wsize > 1) {
1321                 /*
1322                  * W[1 << (wsize - 1)] = W[1] ^ (wsize - 1)
1323                  */
1324                 j = 1 << (wsize - 1);
1325
1326                 MPI_CHK(mpi_grow(&W[j], N->n + 1));
1327                 MPI_CHK(mpi_copy(&W[j], &W[1]));
1328
1329                 for (i = 0; i < wsize - 1; i++)
1330                         mpi_montmul(&W[j], &W[j], N, mm, &T);
1331
1332                 /*
1333                  * W[i] = W[i - 1] * W[1]
1334                  */
1335                 for (i = j + 1; i < (1 << wsize); i++) {
1336                         MPI_CHK(mpi_grow(&W[i], N->n + 1));
1337                         MPI_CHK(mpi_copy(&W[i], &W[i - 1]));
1338
1339                         mpi_montmul(&W[i], &W[1], N, mm, &T);
1340                 }
1341         }
1342
1343         nblimbs = E->n;
1344         bufsize = 0;
1345         nbits = 0;
1346         wbits = 0;
1347         state = 0;
1348
1349         while (1) {
1350                 if (bufsize == 0) {
1351                         if (nblimbs-- == 0)
1352                                 break;
1353
1354                         bufsize = sizeof(t_int) << 3;
1355                 }
1356
1357                 bufsize--;
1358
1359                 ei = (E->p[nblimbs] >> bufsize) & 1;
1360
1361                 /*
1362                  * skip leading 0s
1363                  */
1364                 if (ei == 0 && state == 0)
1365                         continue;
1366
1367                 if (ei == 0 && state == 1) {
1368                         /*
1369                          * out of window, square X
1370                          */
1371                         mpi_montmul(X, X, N, mm, &T);
1372                         continue;
1373                 }
1374
1375                 /*
1376                  * add ei to current window
1377                  */
1378                 state = 2;
1379
1380                 nbits++;
1381                 wbits |= (ei << (wsize - nbits));
1382
1383                 if (nbits == wsize) {
1384                         /*
1385                          * X = X^wsize R^-1 mod N
1386                          */
1387                         for (i = 0; i < wsize; i++)
1388                                 mpi_montmul(X, X, N, mm, &T);
1389
1390                         /*
1391                          * X = X * W[wbits] R^-1 mod N
1392                          */
1393                         mpi_montmul(X, &W[wbits], N, mm, &T);
1394
1395                         state--;
1396                         nbits = 0;
1397                         wbits = 0;
1398                 }
1399         }
1400
1401         /*
1402          * process the remaining bits
1403          */
1404         for (i = 0; i < nbits; i++) {
1405                 mpi_montmul(X, X, N, mm, &T);
1406
1407                 wbits <<= 1;
1408
1409                 if ((wbits & (1 << wsize)) != 0)
1410                         mpi_montmul(X, &W[1], N, mm, &T);
1411         }
1412
1413         /*
1414          * X = A^E * R * R^-1 mod N = A^E mod N
1415          */
1416         mpi_montred(X, N, mm, &T);
1417
1418 cleanup:
1419
1420         for (i = (1 << (wsize - 1)); i < (1 << wsize); i++)
1421                 mpi_free(&W[i]);
1422
1423         mpi_free(&W[1]); mpi_free(&T);
1424         if (_RR == NULL) {
1425                 mpi_free(&RR);
1426         }
1427
1428         return (ret);
1429 }
1430
1431 /*
1432  * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
1433  */
1434 int mpi_gcd(mpi * G, const mpi * A, const mpi * B)
1435 {
1436         int ret, lz, lzt;
1437         mpi TG, TA, TB;
1438
1439         mpi_init(&TG); mpi_init(&TA); mpi_init(&TB);
1440
1441         MPI_CHK(mpi_copy(&TA, A));
1442         MPI_CHK(mpi_copy(&TB, B));
1443
1444         lz = mpi_lsb(&TA);
1445         lzt = mpi_lsb(&TB);
1446
1447         if (lzt < lz)
1448                 lz = lzt;
1449
1450         MPI_CHK(mpi_shift_r(&TA, lz));
1451         MPI_CHK(mpi_shift_r(&TB, lz));
1452
1453         TA.s = TB.s = 1;
1454
1455         while (mpi_cmp_int(&TA, 0) != 0) {
1456                 MPI_CHK(mpi_shift_r(&TA, mpi_lsb(&TA)));
1457                 MPI_CHK(mpi_shift_r(&TB, mpi_lsb(&TB)));
1458
1459                 if (mpi_cmp_mpi(&TA, &TB) >= 0) {
1460                         MPI_CHK(mpi_sub_abs(&TA, &TA, &TB));
1461                         MPI_CHK(mpi_shift_r(&TA, 1));
1462                 } else {
1463                         MPI_CHK(mpi_sub_abs(&TB, &TB, &TA));
1464                         MPI_CHK(mpi_shift_r(&TB, 1));
1465                 }
1466         }
1467
1468         MPI_CHK(mpi_shift_l(&TB, lz));
1469         MPI_CHK(mpi_copy(G, &TB));
1470
1471 cleanup:
1472
1473         mpi_free(&TB); mpi_free(&TA); mpi_free(&TG);
1474
1475         return (ret);
1476 }
1477
1478 #if defined(TROPICSSL_GENPRIME)
1479
1480 /*
1481  * Modular inverse: X = A^-1 mod N      (HAC 14.61 / 14.64)
1482  */
1483 int mpi_inv_mod(mpi * X, const mpi * A, const mpi * N)
1484 {
1485         int ret;
1486         mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1487
1488         if (mpi_cmp_int(N, 0) <= 0)
1489                 return (TROPICSSL_ERR_MPI_BAD_INPUT_DATA);
1490
1491         mpi_init(&TA); mpi_init(&TU); mpi_init(&U1); mpi_init(&U2);
1492         mpi_init(&G); mpi_init(&TB); mpi_init(&TV);
1493         mpi_init(&V1); mpi_init(&V2);
1494
1495         MPI_CHK(mpi_gcd(&G, A, N));
1496
1497         if (mpi_cmp_int(&G, 1) != 0) {
1498                 ret = TROPICSSL_ERR_MPI_NOT_ACCEPTABLE;
1499                 goto cleanup;
1500         }
1501
1502         MPI_CHK(mpi_mod_mpi(&TA, A, N));
1503         MPI_CHK(mpi_copy(&TU, &TA));
1504         MPI_CHK(mpi_copy(&TB, N));
1505         MPI_CHK(mpi_copy(&TV, N));
1506
1507         MPI_CHK(mpi_lset(&U1, 1));
1508         MPI_CHK(mpi_lset(&U2, 0));
1509         MPI_CHK(mpi_lset(&V1, 0));
1510         MPI_CHK(mpi_lset(&V2, 1));
1511
1512         do {
1513                 while ((TU.p[0] & 1) == 0) {
1514                         MPI_CHK(mpi_shift_r(&TU, 1));
1515
1516                         if ((U1.p[0] & 1) != 0 || (U2.p[0] & 1) != 0) {
1517                                 MPI_CHK(mpi_add_mpi(&U1, &U1, &TB));
1518                                 MPI_CHK(mpi_sub_mpi(&U2, &U2, &TA));
1519                         }
1520
1521                         MPI_CHK(mpi_shift_r(&U1, 1));
1522                         MPI_CHK(mpi_shift_r(&U2, 1));
1523                 }
1524
1525                 while ((TV.p[0] & 1) == 0) {
1526                         MPI_CHK(mpi_shift_r(&TV, 1));
1527
1528                         if ((V1.p[0] & 1) != 0 || (V2.p[0] & 1) != 0) {
1529                                 MPI_CHK(mpi_add_mpi(&V1, &V1, &TB));
1530                                 MPI_CHK(mpi_sub_mpi(&V2, &V2, &TA));
1531                         }
1532
1533                         MPI_CHK(mpi_shift_r(&V1, 1));
1534                         MPI_CHK(mpi_shift_r(&V2, 1));
1535                 }
1536
1537                 if (mpi_cmp_mpi(&TU, &TV) >= 0) {
1538                         MPI_CHK(mpi_sub_mpi(&TU, &TU, &TV));
1539                         MPI_CHK(mpi_sub_mpi(&U1, &U1, &V1));
1540                         MPI_CHK(mpi_sub_mpi(&U2, &U2, &V2));
1541                 } else {
1542                         MPI_CHK(mpi_sub_mpi(&TV, &TV, &TU));
1543                         MPI_CHK(mpi_sub_mpi(&V1, &V1, &U1));
1544                         MPI_CHK(mpi_sub_mpi(&V2, &V2, &U2));
1545                 }
1546         } while (mpi_cmp_int(&TU, 0) != 0);
1547
1548         while (mpi_cmp_int(&V1, 0) < 0)
1549                 MPI_CHK(mpi_add_mpi(&V1, &V1, N));
1550
1551         while (mpi_cmp_mpi(&V1, N) >= 0)
1552                 MPI_CHK(mpi_sub_mpi(&V1, &V1, N));
1553
1554         MPI_CHK(mpi_copy(X, &V1));
1555
1556 cleanup:
1557
1558         mpi_free(&V2); mpi_free(&V1); mpi_free(&TV); mpi_free(&TB);
1559         mpi_free(&G); mpi_free(&U2); mpi_free(&U1);
1560         mpi_free(&TU); mpi_free(&TA);
1561
1562         return (ret);
1563 }
1564
1565 static const int small_prime[] = {
1566         3, 5, 7, 11, 13, 17, 19, 23,
1567         29, 31, 37, 41, 43, 47, 53, 59,
1568         61, 67, 71, 73, 79, 83, 89, 97,
1569         101, 103, 107, 109, 113, 127, 131, 137,
1570         139, 149, 151, 157, 163, 167, 173, 179,
1571         181, 191, 193, 197, 199, 211, 223, 227,
1572         229, 233, 239, 241, 251, 257, 263, 269,
1573         271, 277, 281, 283, 293, 307, 311, 313,
1574         317, 331, 337, 347, 349, 353, 359, 367,
1575         373, 379, 383, 389, 397, 401, 409, 419,
1576         421, 431, 433, 439, 443, 449, 457, 461,
1577         463, 467, 479, 487, 491, 499, 503, 509,
1578         521, 523, 541, 547, 557, 563, 569, 571,
1579         577, 587, 593, 599, 601, 607, 613, 617,
1580         619, 631, 641, 643, 647, 653, 659, 661,
1581         673, 677, 683, 691, 701, 709, 719, 727,
1582         733, 739, 743, 751, 757, 761, 769, 773,
1583         787, 797, 809, 811, 821, 823, 827, 829,
1584         839, 853, 857, 859, 863, 877, 881, 883,
1585         887, 907, 911, 919, 929, 937, 941, 947,
1586         953, 967, 971, 977, 983, 991, 997, -103
1587 };
1588
1589 /*
1590  * Miller-Rabin primality test  (HAC 4.24)
1591  */
1592 int mpi_is_prime(mpi * X, int (*f_rng) (void *), void *p_rng)
1593 {
1594         int ret, i, j, n, s, xs;
1595         mpi W, R, T, A, RR;
1596         unsigned char *p;
1597
1598         if (mpi_cmp_int(X, 0) == 0)
1599                 return (0);
1600
1601         mpi_init(&W); mpi_init(&R); mpi_init(&T);
1602         mpi_init(&A); mpi_init(&RR);
1603
1604         xs = X->s;
1605         X->s = 1;
1606
1607         /*
1608          * test trivial factors first
1609          */
1610         if ((X->p[0] & 1) == 0)
1611                 return (TROPICSSL_ERR_MPI_NOT_ACCEPTABLE);
1612
1613         for (i = 0; small_prime[i] > 0; i++) {
1614                 t_int r;
1615
1616                 if (mpi_cmp_int(X, small_prime[i]) <= 0)
1617                         return (0);
1618
1619                 MPI_CHK(mpi_mod_int(&r, X, small_prime[i]));
1620
1621                 if (r == 0)
1622                         return (TROPICSSL_ERR_MPI_NOT_ACCEPTABLE);
1623         }
1624
1625         /*
1626          * W = |X| - 1
1627          * R = W >> lsb( W )
1628          */
1629         s = mpi_lsb(&W);
1630         MPI_CHK(mpi_sub_int(&W, X, 1));
1631         MPI_CHK(mpi_copy(&R, &W));
1632         MPI_CHK(mpi_shift_r(&R, s));
1633
1634         i = mpi_msb(X);
1635         /*
1636          * HAC, table 4.4
1637          */
1638         n = ((i >= 1300) ? 2 : (i >= 850) ? 3 :
1639              (i >= 650) ? 4 : (i >= 350) ? 8 :
1640              (i >= 250) ? 12 : (i >= 150) ? 18 : 27);
1641
1642         for (i = 0; i < n; i++) {
1643                 /*
1644                  * pick a random A, 1 < A < |X| - 1
1645                  */
1646                 MPI_CHK(mpi_grow(&A, X->n));
1647
1648                 p = (unsigned char *)A.p;
1649                 for (j = 0; j < A.n * ciL; j++)
1650                         *p++ = (unsigned char)f_rng(p_rng);
1651
1652                 j = mpi_msb(&A) - mpi_msb(&W);
1653                 MPI_CHK(mpi_shift_r(&A, j + 1));
1654                 A.p[0] |= 3;
1655
1656                 /*
1657                  * A = A^R mod |X|
1658                  */
1659                 MPI_CHK(mpi_exp_mod(&A, &A, &R, X, &RR));
1660
1661                 if (mpi_cmp_mpi(&A, &W) == 0 || mpi_cmp_int(&A, 1) == 0)
1662                         continue;
1663
1664                 j = 1;
1665                 while (j < s && mpi_cmp_mpi(&A, &W) != 0) {
1666                         /*
1667                          * A = A * A mod |X|
1668                          */
1669                         MPI_CHK(mpi_mul_mpi(&T, &A, &A));
1670                         MPI_CHK(mpi_mod_mpi(&A, &T, X));
1671
1672                         if (mpi_cmp_int(&A, 1) == 0)
1673                                 break;
1674
1675                         j++;
1676                 }
1677
1678                 /*
1679                  * not prime if A != |X| - 1 or A == 1
1680                  */
1681                 if (mpi_cmp_mpi(&A, &W) != 0 || mpi_cmp_int(&A, 1) == 0) {
1682                         ret = TROPICSSL_ERR_MPI_NOT_ACCEPTABLE;
1683                         break;
1684                 }
1685         }
1686
1687 cleanup:
1688
1689         X->s = xs;
1690
1691         mpi_free(&RR); mpi_free(&A); mpi_free(&T);
1692         mpi_free(&R); mpi_free(&W);
1693
1694         return (ret);
1695 }
1696
1697 /*
1698  * Prime number generation
1699  */
1700 int mpi_gen_prime(mpi * X, int nbits, int dh_flag,
1701                   int (*f_rng) (void *), void *p_rng)
1702 {
1703         int ret, k, n;
1704         unsigned char *p;
1705         mpi Y;
1706
1707         if (nbits < 3)
1708                 return (TROPICSSL_ERR_MPI_BAD_INPUT_DATA);
1709
1710         mpi_init(&Y);
1711
1712         n = BITS_TO_LIMBS(nbits);
1713
1714         MPI_CHK(mpi_grow(X, n));
1715         MPI_CHK(mpi_lset(X, 0));
1716
1717         p = (unsigned char *)X->p;
1718         for (k = 0; k < X->n * ciL; k++)
1719                 *p++ = (unsigned char)f_rng(p_rng);
1720
1721         k = mpi_msb(X);
1722         if (k < nbits)
1723                 MPI_CHK(mpi_shift_l(X, nbits - k));
1724         if (k > nbits)
1725                 MPI_CHK(mpi_shift_r(X, k - nbits));
1726
1727         X->p[0] |= 3;
1728
1729         if (dh_flag == 0) {
1730                 while ((ret = mpi_is_prime(X, f_rng, p_rng)) != 0) {
1731                         if (ret != TROPICSSL_ERR_MPI_NOT_ACCEPTABLE)
1732                                 goto cleanup;
1733
1734                         MPI_CHK(mpi_add_int(X, X, 2));
1735                 }
1736         } else {
1737                 MPI_CHK(mpi_sub_int(&Y, X, 1));
1738                 MPI_CHK(mpi_shift_r(&Y, 1));
1739
1740                 while (1) {
1741                         if ((ret = mpi_is_prime(X, f_rng, p_rng)) == 0) {
1742                                 if ((ret = mpi_is_prime(&Y, f_rng, p_rng)) == 0)
1743                                         break;
1744
1745                                 if (ret != TROPICSSL_ERR_MPI_NOT_ACCEPTABLE)
1746                                         goto cleanup;
1747                         }
1748
1749                         if (ret != TROPICSSL_ERR_MPI_NOT_ACCEPTABLE)
1750                                 goto cleanup;
1751
1752                         MPI_CHK(mpi_add_int(&Y, X, 1));
1753                         MPI_CHK(mpi_add_int(X, X, 2));
1754                         MPI_CHK(mpi_shift_r(&Y, 1));
1755                 }
1756         }
1757
1758 cleanup:
1759
1760         mpi_free(&Y);
1761
1762         return (ret);
1763 }
1764
1765 #endif
1766
1767 #if defined(TROPICSSL_SELF_TEST)
1768
1769 #define GCD_PAIR_COUNT  3
1770
1771 static const int gcd_pairs[GCD_PAIR_COUNT][3] = {
1772         {693, 609, 21},
1773         {1764, 868, 28},
1774         {768454923, 542167814, 1}
1775 };
1776
1777 /*
1778  * Checkup routine
1779  */
1780 int mpi_self_test(int verbose)
1781 {
1782         int ret, i;
1783         mpi A, E, N, X, Y, U, V;
1784
1785         mpi_init(&A); mpi_init(&E); mpi_init(&N); mpi_init(&X);
1786         mpi_init(&Y); mpi_init(&U); mpi_init(&V);
1787
1788         MPI_CHK(mpi_read_string(&A, 16,
1789                                 "EFE021C2645FD1DC586E69184AF4A31E"
1790                                 "D5F53E93B5F123FA41680867BA110131"
1791                                 "944FE7952E2517337780CB0DB80E61AA"
1792                                 "E7C8DDC6C5C6AADEB34EB38A2F40D5E6"));
1793
1794         MPI_CHK(mpi_read_string(&E, 16,
1795                                 "B2E7EFD37075B9F03FF989C7C5051C20"
1796                                 "34D2A323810251127E7BF8625A4F49A5"
1797                                 "F3E27F4DA8BD59C47D6DAABA4C8127BD"
1798                                 "5B5C25763222FEFCCFC38B832366C29E"));
1799
1800         MPI_CHK(mpi_read_string(&N, 16,
1801                                 "0066A198186C18C10B2F5ED9B522752A"
1802                                 "9830B69916E535C8F047518A889A43A5"
1803                                 "94B6BED27A168D31D4A52F88925AA8F5"));
1804
1805         MPI_CHK(mpi_mul_mpi(&X, &A, &N));
1806
1807         MPI_CHK(mpi_read_string(&U, 16,
1808                                 "602AB7ECA597A3D6B56FF9829A5E8B85"
1809                                 "9E857EA95A03512E2BAE7391688D264A"
1810                                 "A5663B0341DB9CCFD2C4C5F421FEC814"
1811                                 "8001B72E848A38CAE1C65F78E56ABDEF"
1812                                 "E12D3C039B8A02D6BE593F0BBBDA56F1"
1813                                 "ECF677152EF804370C1A305CAF3B5BF1"
1814                                 "30879B56C61DE584A0F53A2447A51E"));
1815
1816         if (verbose != 0)
1817                 printf("  MPI test #1 (mul_mpi): ");
1818
1819         if (mpi_cmp_mpi(&X, &U) != 0) {
1820                 if (verbose != 0)
1821                         printf("failed\n");
1822
1823                 return (1);
1824         }
1825
1826         if (verbose != 0)
1827                 printf("passed\n");
1828
1829         MPI_CHK(mpi_div_mpi(&X, &Y, &A, &N));
1830
1831         MPI_CHK(mpi_read_string(&U, 16, "256567336059E52CAE22925474705F39A94"));
1832
1833         MPI_CHK(mpi_read_string(&V, 16,
1834                                 "6613F26162223DF488E9CD48CC132C7A"
1835                                 "0AC93C701B001B092E4E5B9F73BCD27B"
1836                                 "9EE50D0657C77F374E903CDFA4C642"));
1837
1838         if (verbose != 0)
1839                 printf("  MPI test #2 (div_mpi): ");
1840
1841         if (mpi_cmp_mpi(&X, &U) != 0 || mpi_cmp_mpi(&Y, &V) != 0) {
1842                 if (verbose != 0)
1843                         printf("failed\n");
1844
1845                 return (1);
1846         }
1847
1848         if (verbose != 0)
1849                 printf("passed\n");
1850
1851         MPI_CHK(mpi_exp_mod(&X, &A, &E, &N, NULL));
1852
1853         MPI_CHK(mpi_read_string(&U, 16,
1854                                 "36E139AEA55215609D2816998ED020BB"
1855                                 "BD96C37890F65171D948E9BC7CBAA4D9"
1856                                 "325D24D6A3C12710F10A09FA08AB87"));
1857
1858         if (verbose != 0)
1859                 printf("  MPI test #3 (exp_mod): ");
1860
1861         if (mpi_cmp_mpi(&X, &U) != 0) {
1862                 if (verbose != 0)
1863                         printf("failed\n");
1864
1865                 return (1);
1866         }
1867
1868         if (verbose != 0)
1869                 printf("passed\n");
1870
1871         MPI_CHK(mpi_inv_mod(&X, &A, &N));
1872
1873         MPI_CHK(mpi_read_string(&U, 16,
1874                                 "003A0AAEDD7E784FC07D8F9EC6E3BFD5"
1875                                 "C3DBA76456363A10869622EAC2DD84EC"
1876                                 "C5B8A74DAC4D09E03B5E0BE779F2DF61"));
1877
1878         if (verbose != 0)
1879                 printf("  MPI test #4 (inv_mod): ");
1880
1881         if (mpi_cmp_mpi(&X, &U) != 0) {
1882                 if (verbose != 0)
1883                         printf("failed\n");
1884
1885                 return (1);
1886         }
1887
1888         if (verbose != 0)
1889                 printf("passed\n");
1890
1891         if (verbose != 0)
1892                 printf("  MPI test #5 (simple gcd): ");
1893
1894         for (i = 0; i < GCD_PAIR_COUNT; i++) {
1895                 MPI_CHK(mpi_lset(&X, gcd_pairs[i][0]));
1896                 MPI_CHK(mpi_lset(&Y, gcd_pairs[i][1]));
1897
1898                 MPI_CHK(mpi_gcd(&A, &X, &Y));
1899
1900                 if (mpi_cmp_int(&A, gcd_pairs[i][2]) != 0) {
1901                         if (verbose != 0)
1902                                 printf("failed at %d\n", i);
1903
1904                         return (1);
1905                 }
1906         }
1907
1908         if (verbose != 0)
1909                 printf("passed\n");
1910
1911 cleanup:
1912
1913         if (ret != 0 && verbose != 0)
1914                 printf("Unexpected error, return code = %08X\n", ret);
1915
1916         mpi_free(&V); mpi_free(&U); mpi_free(&Y); mpi_free(&X);
1917         mpi_free(&N); mpi_free(&E); mpi_free(&A);
1918
1919         if (verbose != 0)
1920                 printf("\n");
1921
1922         return (ret);
1923 }
1924
1925 #endif
1926
1927 #endif