polarssl/library/bignum.c

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