src/libpolarssl/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 "config.h"
  34
  35 #if defined(POLARSSL_BIGNUM_C)
  36
  37 #include "bignum.h"
  38 #include "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
 939 #if defined(__APPLE__) && defined(__arm__)
 940 /*
 941  * Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
 942  * appears to need this to prevent bad ARM code generation at -O3.
 943  */
 944 __attribute__ ((noinline))
 945 #endif
 946 void mpi_mul_hlp( size_t i, t_uint *s, t_uint *d, t_uint b )
 947 {
 948     t_uint c = 0, t = 0;
 949
 950 #if defined(MULADDC_HUIT)
 951     for( ; i >= 8; i -= 8 )
 952     {
 953         MULADDC_INIT
 954         MULADDC_HUIT
 955         MULADDC_STOP
 956     }
 957
 958     for( ; i > 0; i-- )
 959     {
 960         MULADDC_INIT
 961         MULADDC_CORE
 962         MULADDC_STOP
 963     }
 964 #else
 965     for( ; i >= 16; i -= 16 )
 966     {
 967         MULADDC_INIT
 968         MULADDC_CORE   MULADDC_CORE
 969         MULADDC_CORE   MULADDC_CORE
 970         MULADDC_CORE   MULADDC_CORE
 971         MULADDC_CORE   MULADDC_CORE
 972
 973         MULADDC_CORE   MULADDC_CORE
 974         MULADDC_CORE   MULADDC_CORE
 975         MULADDC_CORE   MULADDC_CORE
 976         MULADDC_CORE   MULADDC_CORE
 977         MULADDC_STOP
 978     }
 979
 980     for( ; i >= 8; i -= 8 )
 981     {
 982         MULADDC_INIT
 983         MULADDC_CORE   MULADDC_CORE
 984         MULADDC_CORE   MULADDC_CORE
 985
 986         MULADDC_CORE   MULADDC_CORE
 987         MULADDC_CORE   MULADDC_CORE
 988         MULADDC_STOP
 989     }
 990
 991     for( ; i > 0; i-- )
 992     {
 993         MULADDC_INIT
 994         MULADDC_CORE
 995         MULADDC_STOP
 996     }
 997 #endif
 998
 999     t++;
1000
1001     do {
1002         *d += c; c = ( *d < c ); d++;
1003     }
1004     while( c != 0 );
1005 }
1006
1007 /*
1008  * Baseline multiplication: X = A * B  (HAC 14.12)
1009  */
1010 int mpi_mul_mpi( mpi *X, const mpi *A, const mpi *B )
1011 {
1012     int ret;
1013     size_t i, j;
1014     mpi TA, TB;
1015
1016     mpi_init( &TA ); mpi_init( &TB );
1017
1018     if( X == A ) { MPI_CHK( mpi_copy( &TA, A ) ); A = &TA; }
1019     if( X == B ) { MPI_CHK( mpi_copy( &TB, B ) ); B = &TB; }
1020
1021     for( i = A->n; i > 0; i-- )
1022         if( A->p[i - 1] != 0 )
1023             break;
1024
1025     for( j = B->n; j > 0; j-- )
1026         if( B->p[j - 1] != 0 )
1027             break;
1028
1029     MPI_CHK( mpi_grow( X, i + j ) );
1030     MPI_CHK( mpi_lset( X, 0 ) );
1031
1032     for( i++; j > 0; j-- )
1033         mpi_mul_hlp( i - 1, A->p, X->p + j - 1, B->p[j - 1] );
1034
1035     X->s = A->s * B->s;
1036
1037 cleanup:
1038
1039     mpi_free( &TB ); mpi_free( &TA );
1040
1041     return( ret );
1042 }
1043
1044 /*
1045  * Baseline multiplication: X = A * b
1046  */
1047 int mpi_mul_int( mpi *X, const mpi *A, t_sint b )
1048 {
1049     mpi _B;
1050     t_uint p[1];
1051
1052     _B.s = 1;
1053     _B.n = 1;
1054     _B.p = p;
1055     p[0] = b;
1056
1057     return( mpi_mul_mpi( X, A, &_B ) );
1058 }
1059
1060 /*
1061  * Division by mpi: A = Q * B + R  (HAC 14.20)
1062  */
1063 int mpi_div_mpi( mpi *Q, mpi *R, const mpi *A, const mpi *B )
1064 {
1065     int ret;
1066     size_t i, n, t, k;
1067     mpi X, Y, Z, T1, T2;
1068
1069     if( mpi_cmp_int( B, 0 ) == 0 )
1070         return( POLARSSL_ERR_MPI_DIVISION_BY_ZERO );
1071
1072     mpi_init( &X ); mpi_init( &Y ); mpi_init( &Z );
1073     mpi_init( &T1 ); mpi_init( &T2 );
1074
1075     if( mpi_cmp_abs( A, B ) < 0 )
1076     {
1077         if( Q != NULL ) MPI_CHK( mpi_lset( Q, 0 ) );
1078         if( R != NULL ) MPI_CHK( mpi_copy( R, A ) );
1079         return( 0 );
1080     }
1081
1082     MPI_CHK( mpi_copy( &X, A ) );
1083     MPI_CHK( mpi_copy( &Y, B ) );
1084     X.s = Y.s = 1;
1085
1086     MPI_CHK( mpi_grow( &Z, A->n + 2 ) );
1087     MPI_CHK( mpi_lset( &Z,  0 ) );
1088     MPI_CHK( mpi_grow( &T1, 2 ) );
1089     MPI_CHK( mpi_grow( &T2, 3 ) );
1090
1091     k = mpi_msb( &Y ) % biL;
1092     if( k < biL - 1 )
1093     {
1094         k = biL - 1 - k;
1095         MPI_CHK( mpi_shift_l( &X, k ) );
1096         MPI_CHK( mpi_shift_l( &Y, k ) );
1097     }
1098     else k = 0;
1099
1100     n = X.n - 1;
1101     t = Y.n - 1;
1102     MPI_CHK( mpi_shift_l( &Y, biL * (n - t) ) );
1103
1104     while( mpi_cmp_mpi( &X, &Y ) >= 0 )
1105     {
1106         Z.p[n - t]++;
1107         mpi_sub_mpi( &X, &X, &Y );
1108     }
1109     mpi_shift_r( &Y, biL * (n - t) );
1110
1111     for( i = n; i > t ; i-- )
1112     {
1113         if( X.p[i] >= Y.p[t] )
1114             Z.p[i - t - 1] = ~0;
1115         else
1116         {
1117 #if defined(POLARSSL_HAVE_UDBL)
1118             t_udbl r;
1119
1120             r  = (t_udbl) X.p[i] << biL;
1121             r |= (t_udbl) X.p[i - 1];
1122             r /= Y.p[t];
1123             if( r > ((t_udbl) 1 << biL) - 1)
1124                 r = ((t_udbl) 1 << biL) - 1;
1125
1126             Z.p[i - t - 1] = (t_uint) r;
1127 #else
1128             /*
1129              * __udiv_qrnnd_c, from gmp/longlong.h
1130              */
1131             t_uint q0, q1, r0, r1;
1132             t_uint d0, d1, d, m;
1133
1134             d  = Y.p[t];
1135             d0 = ( d << biH ) >> biH;
1136             d1 = ( d >> biH );
1137
1138             q1 = X.p[i] / d1;
1139             r1 = X.p[i] - d1 * q1;
1140             r1 <<= biH;
1141             r1 |= ( X.p[i - 1] >> biH );
1142
1143             m = q1 * d0;
1144             if( r1 < m )
1145             {
1146                 q1--, r1 += d;
1147                 while( r1 >= d && r1 < m )
1148                     q1--, r1 += d;
1149             }
1150             r1 -= m;
1151
1152             q0 = r1 / d1;
1153             r0 = r1 - d1 * q0;
1154             r0 <<= biH;
1155             r0 |= ( X.p[i - 1] << biH ) >> biH;
1156
1157             m = q0 * d0;
1158             if( r0 < m )
1159             {
1160                 q0--, r0 += d;
1161                 while( r0 >= d && r0 < m )
1162                     q0--, r0 += d;
1163             }
1164             r0 -= m;
1165
1166             Z.p[i - t - 1] = ( q1 << biH ) | q0;
1167 #endif
1168         }
1169
1170         Z.p[i - t - 1]++;
1171         do
1172         {
1173             Z.p[i - t - 1]--;
1174
1175             MPI_CHK( mpi_lset( &T1, 0 ) );
1176             T1.p[0] = (t < 1) ? 0 : Y.p[t - 1];
1177             T1.p[1] = Y.p[t];
1178             MPI_CHK( mpi_mul_int( &T1, &T1, Z.p[i - t - 1] ) );
1179
1180             MPI_CHK( mpi_lset( &T2, 0 ) );
1181             T2.p[0] = (i < 2) ? 0 : X.p[i - 2];
1182             T2.p[1] = (i < 1) ? 0 : X.p[i - 1];
1183             T2.p[2] = X.p[i];
1184         }
1185         while( mpi_cmp_mpi( &T1, &T2 ) > 0 );
1186
1187         MPI_CHK( mpi_mul_int( &T1, &Y, Z.p[i - t - 1] ) );
1188         MPI_CHK( mpi_shift_l( &T1,  biL * (i - t - 1) ) );
1189         MPI_CHK( mpi_sub_mpi( &X, &X, &T1 ) );
1190
1191         if( mpi_cmp_int( &X, 0 ) < 0 )
1192         {
1193             MPI_CHK( mpi_copy( &T1, &Y ) );
1194             MPI_CHK( mpi_shift_l( &T1, biL * (i - t - 1) ) );
1195             MPI_CHK( mpi_add_mpi( &X, &X, &T1 ) );
1196             Z.p[i - t - 1]--;
1197         }
1198     }
1199
1200     if( Q != NULL )
1201     {
1202         mpi_copy( Q, &Z );
1203         Q->s = A->s * B->s;
1204     }
1205
1206     if( R != NULL )
1207     {
1208         mpi_shift_r( &X, k );
1209         X.s = A->s;
1210         mpi_copy( R, &X );
1211
1212         if( mpi_cmp_int( R, 0 ) == 0 )
1213             R->s = 1;
1214     }
1215
1216 cleanup:
1217
1218     mpi_free( &X ); mpi_free( &Y ); mpi_free( &Z );
1219     mpi_free( &T1 ); mpi_free( &T2 );
1220
1221     return( ret );
1222 }
1223
1224 /*
1225  * Division by int: A = Q * b + R
1226  */
1227 int mpi_div_int( mpi *Q, mpi *R, const mpi *A, t_sint b )
1228 {
1229     mpi _B;
1230     t_uint p[1];
1231
1232     p[0] = ( b < 0 ) ? -b : b;
1233     _B.s = ( b < 0 ) ? -1 : 1;
1234     _B.n = 1;
1235     _B.p = p;
1236
1237     return( mpi_div_mpi( Q, R, A, &_B ) );
1238 }
1239
1240 /*
1241  * Modulo: R = A mod B
1242  */
1243 int mpi_mod_mpi( mpi *R, const mpi *A, const mpi *B )
1244 {
1245     int ret;
1246
1247     if( mpi_cmp_int( B, 0 ) < 0 )
1248         return POLARSSL_ERR_MPI_NEGATIVE_VALUE;
1249
1250     MPI_CHK( mpi_div_mpi( NULL, R, A, B ) );
1251
1252     while( mpi_cmp_int( R, 0 ) < 0 )
1253       MPI_CHK( mpi_add_mpi( R, R, B ) );
1254
1255     while( mpi_cmp_mpi( R, B ) >= 0 )
1256       MPI_CHK( mpi_sub_mpi( R, R, B ) );
1257
1258 cleanup:
1259
1260     return( ret );
1261 }
1262
1263 /*
1264  * Modulo: r = A mod b
1265  */
1266 int mpi_mod_int( t_uint *r, const mpi *A, t_sint b )
1267 {
1268     size_t i;
1269     t_uint x, y, z;
1270
1271     if( b == 0 )
1272         return( POLARSSL_ERR_MPI_DIVISION_BY_ZERO );
1273
1274     if( b < 0 )
1275         return POLARSSL_ERR_MPI_NEGATIVE_VALUE;
1276
1277     /*
1278      * handle trivial cases
1279      */
1280     if( b == 1 )
1281     {
1282         *r = 0;
1283         return( 0 );
1284     }
1285
1286     if( b == 2 )
1287     {
1288         *r = A->p[0] & 1;
1289         return( 0 );
1290     }
1291
1292     /*
1293      * general case
1294      */
1295     for( i = A->n, y = 0; i > 0; i-- )
1296     {
1297         x  = A->p[i - 1];
1298         y  = ( y << biH ) | ( x >> biH );
1299         z  = y / b;
1300         y -= z * b;
1301
1302         x <<= biH;
1303         y  = ( y << biH ) | ( x >> biH );
1304         z  = y / b;
1305         y -= z * b;
1306     }
1307
1308     /*
1309      * If A is negative, then the current y represents a negative value.
1310      * Flipping it to the positive side.
1311      */
1312     if( A->s < 0 && y != 0 )
1313         y = b - y;
1314
1315     *r = y;
1316
1317     return( 0 );
1318 }
1319
1320 /*
1321  * Fast Montgomery initialization (thanks to Tom St Denis)
1322  */
1323 static void mpi_montg_init( t_uint *mm, const mpi *N )
1324 {
1325     t_uint x, m0 = N->p[0];
1326
1327     x  = m0;
1328     x += ( ( m0 + 2 ) & 4 ) << 1;
1329     x *= ( 2 - ( m0 * x ) );
1330
1331     if( biL >= 16 ) x *= ( 2 - ( m0 * x ) );
1332     if( biL >= 32 ) x *= ( 2 - ( m0 * x ) );
1333     if( biL >= 64 ) x *= ( 2 - ( m0 * x ) );
1334
1335     *mm = ~x + 1;
1336 }
1337
1338 /*
1339  * Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
1340  */
1341 static void mpi_montmul( mpi *A, const mpi *B, const mpi *N, t_uint mm, const mpi *T )
1342 {
1343     size_t i, n, m;
1344     t_uint u0, u1, *d;
1345
1346     memset( T->p, 0, T->n * ciL );
1347
1348     d = T->p;
1349     n = N->n;
1350     m = ( B->n < n ) ? B->n : n;
1351
1352     for( i = 0; i < n; i++ )
1353     {
1354         /*
1355          * T = (T + u0*B + u1*N) / 2^biL
1356          */
1357         u0 = A->p[i];
1358         u1 = ( d[0] + u0 * B->p[0] ) * mm;
1359
1360         mpi_mul_hlp( m, B->p, d, u0 );
1361         mpi_mul_hlp( n, N->p, d, u1 );
1362
1363         *d++ = u0; d[n + 1] = 0;
1364     }
1365
1366     memcpy( A->p, d, (n + 1) * ciL );
1367
1368     if( mpi_cmp_abs( A, N ) >= 0 )
1369         mpi_sub_hlp( n, N->p, A->p );
1370     else
1371         /* prevent timing attacks */
1372         mpi_sub_hlp( n, A->p, T->p );
1373 }
1374
1375 /*
1376  * Montgomery reduction: A = A * R^-1 mod N
1377  */
1378 static void mpi_montred( mpi *A, const mpi *N, t_uint mm, const mpi *T )
1379 {
1380     t_uint z = 1;
1381     mpi U;
1382
1383     U.n = U.s = (int) z;
1384     U.p = &z;
1385
1386     mpi_montmul( A, &U, N, mm, T );
1387 }
1388
1389 /*
1390  * Sliding-window exponentiation: X = A^E mod N  (HAC 14.85)
1391  */
1392 int mpi_exp_mod( mpi *X, const mpi *A, const mpi *E, const mpi *N, mpi *_RR )
1393 {
1394     int ret;
1395     size_t wbits, wsize, one = 1;
1396     size_t i, j, nblimbs;
1397     size_t bufsize, nbits;
1398     t_uint ei, mm, state;
1399     mpi RR, T, W[ 2 << POLARSSL_MPI_WINDOW_SIZE ], Apos;
1400     int neg;
1401
1402     if( mpi_cmp_int( N, 0 ) < 0 || ( N->p[0] & 1 ) == 0 )
1403         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1404
1405     if( mpi_cmp_int( E, 0 ) < 0 )
1406         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1407
1408     /*
1409      * Init temps and window size
1410      */
1411     mpi_montg_init( &mm, N );
1412     mpi_init( &RR ); mpi_init( &T );
1413     memset( W, 0, sizeof( W ) );
1414
1415     i = mpi_msb( E );
1416
1417     wsize = ( i > 671 ) ? 6 : ( i > 239 ) ? 5 :
1418             ( i >  79 ) ? 4 : ( i >  23 ) ? 3 : 1;
1419
1420     if( wsize > POLARSSL_MPI_WINDOW_SIZE )
1421         wsize = POLARSSL_MPI_WINDOW_SIZE;
1422
1423     j = N->n + 1;
1424     MPI_CHK( mpi_grow( X, j ) );
1425     MPI_CHK( mpi_grow( &W[1],  j ) );
1426     MPI_CHK( mpi_grow( &T, j * 2 ) );
1427
1428     /*
1429      * Compensate for negative A (and correct at the end)
1430      */
1431     neg = ( A->s == -1 );
1432
1433     mpi_init( &Apos );
1434     if( neg )
1435     {
1436         MPI_CHK( mpi_copy( &Apos, A ) );
1437         Apos.s = 1;
1438         A = &Apos;
1439     }
1440
1441     /*
1442      * If 1st call, pre-compute R^2 mod N
1443      */
1444     if( _RR == NULL || _RR->p == NULL )
1445     {
1446         MPI_CHK( mpi_lset( &RR, 1 ) );
1447         MPI_CHK( mpi_shift_l( &RR, N->n * 2 * biL ) );
1448         MPI_CHK( mpi_mod_mpi( &RR, &RR, N ) );
1449
1450         if( _RR != NULL )
1451             memcpy( _RR, &RR, sizeof( mpi ) );
1452     }
1453     else
1454         memcpy( &RR, _RR, sizeof( mpi ) );
1455
1456     /*
1457      * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1458      */
1459     if( mpi_cmp_mpi( A, N ) >= 0 )
1460         mpi_mod_mpi( &W[1], A, N );
1461     else   mpi_copy( &W[1], A );
1462
1463     mpi_montmul( &W[1], &RR, N, mm, &T );
1464
1465     /*
1466      * X = R^2 * R^-1 mod N = R mod N
1467      */
1468     MPI_CHK( mpi_copy( X, &RR ) );
1469     mpi_montred( X, N, mm, &T );
1470
1471     if( wsize > 1 )
1472     {
1473         /*
1474          * W[1 << (wsize - 1)] = W[1] ^ (wsize - 1)
1475          */
1476         j =  one << (wsize - 1);
1477
1478         MPI_CHK( mpi_grow( &W[j], N->n + 1 ) );
1479         MPI_CHK( mpi_copy( &W[j], &W[1]    ) );
1480
1481         for( i = 0; i < wsize - 1; i++ )
1482             mpi_montmul( &W[j], &W[j], N, mm, &T );
1483
1484         /*
1485          * W[i] = W[i - 1] * W[1]
1486          */
1487         for( i = j + 1; i < (one << wsize); i++ )
1488         {
1489             MPI_CHK( mpi_grow( &W[i], N->n + 1 ) );
1490             MPI_CHK( mpi_copy( &W[i], &W[i - 1] ) );
1491
1492             mpi_montmul( &W[i], &W[1], N, mm, &T );
1493         }
1494     }
1495
1496     nblimbs = E->n;
1497     bufsize = 0;
1498     nbits   = 0;
1499     wbits   = 0;
1500     state   = 0;
1501
1502     while( 1 )
1503     {
1504         if( bufsize == 0 )
1505         {
1506             if( nblimbs-- == 0 )
1507                 break;
1508
1509             bufsize = sizeof( t_uint ) << 3;
1510         }
1511
1512         bufsize--;
1513
1514         ei = (E->p[nblimbs] >> bufsize) & 1;
1515
1516         /*
1517          * skip leading 0s
1518          */
1519         if( ei == 0 && state == 0 )
1520             continue;
1521
1522         if( ei == 0 && state == 1 )
1523         {
1524             /*
1525              * out of window, square X
1526              */
1527             mpi_montmul( X, X, N, mm, &T );
1528             continue;
1529         }
1530
1531         /*
1532          * add ei to current window
1533          */
1534         state = 2;
1535
1536         nbits++;
1537         wbits |= (ei << (wsize - nbits));
1538
1539         if( nbits == wsize )
1540         {
1541             /*
1542              * X = X^wsize R^-1 mod N
1543              */
1544             for( i = 0; i < wsize; i++ )
1545                 mpi_montmul( X, X, N, mm, &T );
1546
1547             /*
1548              * X = X * W[wbits] R^-1 mod N
1549              */
1550             mpi_montmul( X, &W[wbits], N, mm, &T );
1551
1552             state--;
1553             nbits = 0;
1554             wbits = 0;
1555         }
1556     }
1557
1558     /*
1559      * process the remaining bits
1560      */
1561     for( i = 0; i < nbits; i++ )
1562     {
1563         mpi_montmul( X, X, N, mm, &T );
1564
1565         wbits <<= 1;
1566
1567         if( (wbits & (one << wsize)) != 0 )
1568             mpi_montmul( X, &W[1], N, mm, &T );
1569     }
1570
1571     /*
1572      * X = A^E * R * R^-1 mod N = A^E mod N
1573      */
1574     mpi_montred( X, N, mm, &T );
1575
1576     if( neg )
1577     {
1578         X->s = -1;
1579         mpi_add_mpi( X, N, X );
1580     }
1581
1582 cleanup:
1583
1584     for( i = (one << (wsize - 1)); i < (one << wsize); i++ )
1585         mpi_free( &W[i] );
1586
1587     mpi_free( &W[1] ); mpi_free( &T ); mpi_free( &Apos );
1588
1589     if( _RR == NULL )
1590         mpi_free( &RR );
1591
1592     return( ret );
1593 }
1594
1595 /*
1596  * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
1597  */
1598 int mpi_gcd( mpi *G, const mpi *A, const mpi *B )
1599 {
1600     int ret;
1601     size_t lz, lzt;
1602     mpi TG, TA, TB;
1603
1604     mpi_init( &TG ); mpi_init( &TA ); mpi_init( &TB );
1605
1606     MPI_CHK( mpi_copy( &TA, A ) );
1607     MPI_CHK( mpi_copy( &TB, B ) );
1608
1609     lz = mpi_lsb( &TA );
1610     lzt = mpi_lsb( &TB );
1611
1612     if ( lzt < lz )
1613         lz = lzt;
1614
1615     MPI_CHK( mpi_shift_r( &TA, lz ) );
1616     MPI_CHK( mpi_shift_r( &TB, lz ) );
1617
1618     TA.s = TB.s = 1;
1619
1620     while( mpi_cmp_int( &TA, 0 ) != 0 )
1621     {
1622         MPI_CHK( mpi_shift_r( &TA, mpi_lsb( &TA ) ) );
1623         MPI_CHK( mpi_shift_r( &TB, mpi_lsb( &TB ) ) );
1624
1625         if( mpi_cmp_mpi( &TA, &TB ) >= 0 )
1626         {
1627             MPI_CHK( mpi_sub_abs( &TA, &TA, &TB ) );
1628             MPI_CHK( mpi_shift_r( &TA, 1 ) );
1629         }
1630         else
1631         {
1632             MPI_CHK( mpi_sub_abs( &TB, &TB, &TA ) );
1633             MPI_CHK( mpi_shift_r( &TB, 1 ) );
1634         }
1635     }
1636
1637     MPI_CHK( mpi_shift_l( &TB, lz ) );
1638     MPI_CHK( mpi_copy( G, &TB ) );
1639
1640 cleanup:
1641
1642     mpi_free( &TG ); mpi_free( &TA ); mpi_free( &TB );
1643
1644     return( ret );
1645 }
1646
1647 int mpi_fill_random( mpi *X, size_t size,
1648                      int (*f_rng)(void *, unsigned char *, size_t),
1649                      void *p_rng )
1650 {
1651     int ret;
1652
1653     MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( size ) ) );
1654     MPI_CHK( mpi_lset( X, 0 ) );
1655
1656     MPI_CHK( f_rng( p_rng, (unsigned char *) X->p, size ) );
1657
1658 cleanup:
1659     return( ret );
1660 }
1661
1662 /*
1663  * Modular inverse: X = A^-1 mod N  (HAC 14.61 / 14.64)
1664  */
1665 int mpi_inv_mod( mpi *X, const mpi *A, const mpi *N )
1666 {
1667     int ret;
1668     mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1669
1670     if( mpi_cmp_int( N, 0 ) <= 0 )
1671         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1672
1673     mpi_init( &TA ); mpi_init( &TU ); mpi_init( &U1 ); mpi_init( &U2 );
1674     mpi_init( &G ); mpi_init( &TB ); mpi_init( &TV );
1675     mpi_init( &V1 ); mpi_init( &V2 );
1676
1677     MPI_CHK( mpi_gcd( &G, A, N ) );
1678
1679     if( mpi_cmp_int( &G, 1 ) != 0 )
1680     {
1681         ret = POLARSSL_ERR_MPI_NOT_ACCEPTABLE;
1682         goto cleanup;
1683     }
1684
1685     MPI_CHK( mpi_mod_mpi( &TA, A, N ) );
1686     MPI_CHK( mpi_copy( &TU, &TA ) );
1687     MPI_CHK( mpi_copy( &TB, N ) );
1688     MPI_CHK( mpi_copy( &TV, N ) );
1689
1690     MPI_CHK( mpi_lset( &U1, 1 ) );
1691     MPI_CHK( mpi_lset( &U2, 0 ) );
1692     MPI_CHK( mpi_lset( &V1, 0 ) );
1693     MPI_CHK( mpi_lset( &V2, 1 ) );
1694
1695     do
1696     {
1697         while( ( TU.p[0] & 1 ) == 0 )
1698         {
1699             MPI_CHK( mpi_shift_r( &TU, 1 ) );
1700
1701             if( ( U1.p[0] & 1 ) != 0 || ( U2.p[0] & 1 ) != 0 )
1702             {
1703                 MPI_CHK( mpi_add_mpi( &U1, &U1, &TB ) );
1704                 MPI_CHK( mpi_sub_mpi( &U2, &U2, &TA ) );
1705             }
1706
1707             MPI_CHK( mpi_shift_r( &U1, 1 ) );
1708             MPI_CHK( mpi_shift_r( &U2, 1 ) );
1709         }
1710
1711         while( ( TV.p[0] & 1 ) == 0 )
1712         {
1713             MPI_CHK( mpi_shift_r( &TV, 1 ) );
1714
1715             if( ( V1.p[0] & 1 ) != 0 || ( V2.p[0] & 1 ) != 0 )
1716             {
1717                 MPI_CHK( mpi_add_mpi( &V1, &V1, &TB ) );
1718                 MPI_CHK( mpi_sub_mpi( &V2, &V2, &TA ) );
1719             }
1720
1721             MPI_CHK( mpi_shift_r( &V1, 1 ) );
1722             MPI_CHK( mpi_shift_r( &V2, 1 ) );
1723         }
1724
1725         if( mpi_cmp_mpi( &TU, &TV ) >= 0 )
1726         {
1727             MPI_CHK( mpi_sub_mpi( &TU, &TU, &TV ) );
1728             MPI_CHK( mpi_sub_mpi( &U1, &U1, &V1 ) );
1729             MPI_CHK( mpi_sub_mpi( &U2, &U2, &V2 ) );
1730         }
1731         else
1732         {
1733             MPI_CHK( mpi_sub_mpi( &TV, &TV, &TU ) );
1734             MPI_CHK( mpi_sub_mpi( &V1, &V1, &U1 ) );
1735             MPI_CHK( mpi_sub_mpi( &V2, &V2, &U2 ) );
1736         }
1737     }
1738     while( mpi_cmp_int( &TU, 0 ) != 0 );
1739
1740     while( mpi_cmp_int( &V1, 0 ) < 0 )
1741         MPI_CHK( mpi_add_mpi( &V1, &V1, N ) );
1742
1743     while( mpi_cmp_mpi( &V1, N ) >= 0 )
1744         MPI_CHK( mpi_sub_mpi( &V1, &V1, N ) );
1745
1746     MPI_CHK( mpi_copy( X, &V1 ) );
1747
1748 cleanup:
1749
1750     mpi_free( &TA ); mpi_free( &TU ); mpi_free( &U1 ); mpi_free( &U2 );
1751     mpi_free( &G ); mpi_free( &TB ); mpi_free( &TV );
1752     mpi_free( &V1 ); mpi_free( &V2 );
1753
1754     return( ret );
1755 }
1756
1757 #if defined(POLARSSL_GENPRIME)
1758
1759 static const int small_prime[] =
1760 {
1761         3,    5,    7,   11,   13,   17,   19,   23,
1762        29,   31,   37,   41,   43,   47,   53,   59,
1763        61,   67,   71,   73,   79,   83,   89,   97,
1764       101,  103,  107,  109,  113,  127,  131,  137,
1765       139,  149,  151,  157,  163,  167,  173,  179,
1766       181,  191,  193,  197,  199,  211,  223,  227,
1767       229,  233,  239,  241,  251,  257,  263,  269,
1768       271,  277,  281,  283,  293,  307,  311,  313,
1769       317,  331,  337,  347,  349,  353,  359,  367,
1770       373,  379,  383,  389,  397,  401,  409,  419,
1771       421,  431,  433,  439,  443,  449,  457,  461,
1772       463,  467,  479,  487,  491,  499,  503,  509,
1773       521,  523,  541,  547,  557,  563,  569,  571,
1774       577,  587,  593,  599,  601,  607,  613,  617,
1775       619,  631,  641,  643,  647,  653,  659,  661,
1776       673,  677,  683,  691,  701,  709,  719,  727,
1777       733,  739,  743,  751,  757,  761,  769,  773,
1778       787,  797,  809,  811,  821,  823,  827,  829,
1779       839,  853,  857,  859,  863,  877,  881,  883,
1780       887,  907,  911,  919,  929,  937,  941,  947,
1781       953,  967,  971,  977,  983,  991,  997, -103
1782 };
1783
1784 /*
1785  * Miller-Rabin primality test  (HAC 4.24)
1786  */
1787 int mpi_is_prime( mpi *X,
1788                   int (*f_rng)(void *, unsigned char *, size_t),
1789                   void *p_rng )
1790 {
1791     int ret, xs;
1792     size_t i, j, n, s;
1793     mpi W, R, T, A, RR;
1794
1795     if( mpi_cmp_int( X, 0 ) == 0 ||
1796         mpi_cmp_int( X, 1 ) == 0 )
1797         return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
1798
1799     if( mpi_cmp_int( X, 2 ) == 0 )
1800         return( 0 );
1801
1802     mpi_init( &W ); mpi_init( &R ); mpi_init( &T ); mpi_init( &A );
1803     mpi_init( &RR );
1804
1805     xs = X->s; X->s = 1;
1806
1807     /*
1808      * test trivial factors first
1809      */
1810     if( ( X->p[0] & 1 ) == 0 )
1811         return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
1812
1813     for( i = 0; small_prime[i] > 0; i++ )
1814     {
1815         t_uint r;
1816
1817         if( mpi_cmp_int( X, small_prime[i] ) <= 0 )
1818             return( 0 );
1819
1820         MPI_CHK( mpi_mod_int( &r, X, small_prime[i] ) );
1821
1822         if( r == 0 )
1823             return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
1824     }
1825
1826     /*
1827      * W = |X| - 1
1828      * R = W >> lsb( W )
1829      */
1830     MPI_CHK( mpi_sub_int( &W, X, 1 ) );
1831     s = mpi_lsb( &W );
1832     MPI_CHK( mpi_copy( &R, &W ) );
1833     MPI_CHK( mpi_shift_r( &R, s ) );
1834
1835     i = mpi_msb( X );
1836     /*
1837      * HAC, table 4.4
1838      */
1839     n = ( ( i >= 1300 ) ?  2 : ( i >=  850 ) ?  3 :
1840           ( i >=  650 ) ?  4 : ( i >=  350 ) ?  8 :
1841           ( i >=  250 ) ? 12 : ( i >=  150 ) ? 18 : 27 );
1842
1843     for( i = 0; i < n; i++ )
1844     {
1845         /*
1846          * pick a random A, 1 < A < |X| - 1
1847          */
1848         MPI_CHK( mpi_fill_random( &A, X->n * ciL, f_rng, p_rng ) );
1849
1850         if( mpi_cmp_mpi( &A, &W ) >= 0 )
1851         {
1852             j = mpi_msb( &A ) - mpi_msb( &W );
1853             MPI_CHK( mpi_shift_r( &A, j + 1 ) );
1854         }
1855         A.p[0] |= 3;
1856
1857         /*
1858          * A = A^R mod |X|
1859          */
1860         MPI_CHK( mpi_exp_mod( &A, &A, &R, X, &RR ) );
1861
1862         if( mpi_cmp_mpi( &A, &W ) == 0 ||
1863             mpi_cmp_int( &A,  1 ) == 0 )
1864             continue;
1865
1866         j = 1;
1867         while( j < s && mpi_cmp_mpi( &A, &W ) != 0 )
1868         {
1869             /*
1870              * A = A * A mod |X|
1871              */
1872             MPI_CHK( mpi_mul_mpi( &T, &A, &A ) );
1873             MPI_CHK( mpi_mod_mpi( &A, &T, X  ) );
1874
1875             if( mpi_cmp_int( &A, 1 ) == 0 )
1876                 break;
1877
1878             j++;
1879         }
1880
1881         /*
1882          * not prime if A != |X| - 1 or A == 1
1883          */
1884         if( mpi_cmp_mpi( &A, &W ) != 0 ||
1885             mpi_cmp_int( &A,  1 ) == 0 )
1886         {
1887             ret = POLARSSL_ERR_MPI_NOT_ACCEPTABLE;
1888             break;
1889         }
1890     }
1891
1892 cleanup:
1893
1894     X->s = xs;
1895
1896     mpi_free( &W ); mpi_free( &R ); mpi_free( &T ); mpi_free( &A );
1897     mpi_free( &RR );
1898
1899     return( ret );
1900 }
1901
1902 /*
1903  * Prime number generation
1904  */
1905 int mpi_gen_prime( mpi *X, size_t nbits, int dh_flag,
1906                    int (*f_rng)(void *, unsigned char *, size_t),
1907                    void *p_rng )
1908 {
1909     int ret;
1910     size_t k, n;
1911     mpi Y;
1912
1913     if( nbits < 3 || nbits > POLARSSL_MPI_MAX_BITS )
1914         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1915
1916     mpi_init( &Y );
1917
1918     n = BITS_TO_LIMBS( nbits );
1919
1920     MPI_CHK( mpi_fill_random( X, n * ciL, f_rng, p_rng ) );
1921
1922     k = mpi_msb( X );
1923     if( k < nbits ) MPI_CHK( mpi_shift_l( X, nbits - k ) );
1924     if( k > nbits ) MPI_CHK( mpi_shift_r( X, k - nbits ) );
1925
1926     X->p[0] |= 3;
1927
1928     if( dh_flag == 0 )
1929     {
1930         while( ( ret = mpi_is_prime( X, f_rng, p_rng ) ) != 0 )
1931         {
1932             if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1933                 goto cleanup;
1934
1935             MPI_CHK( mpi_add_int( X, X, 2 ) );
1936         }
1937     }
1938     else
1939     {
1940         MPI_CHK( mpi_sub_int( &Y, X, 1 ) );
1941         MPI_CHK( mpi_shift_r( &Y, 1 ) );
1942
1943         while( 1 )
1944         {
1945             if( ( ret = mpi_is_prime( X, f_rng, p_rng ) ) == 0 )
1946             {
1947                 if( ( ret = mpi_is_prime( &Y, f_rng, p_rng ) ) == 0 )
1948                     break;
1949
1950                 if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1951                     goto cleanup;
1952             }
1953
1954             if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1955                 goto cleanup;
1956
1957             MPI_CHK( mpi_add_int( &Y, X, 1 ) );
1958             MPI_CHK( mpi_add_int(  X, X, 2 ) );
1959             MPI_CHK( mpi_shift_r( &Y, 1 ) );
1960         }
1961     }
1962
1963 cleanup:
1964
1965     mpi_free( &Y );
1966
1967     return( ret );
1968 }
1969
1970 #endif
1971
1972 #if defined(POLARSSL_SELF_TEST)
1973
1974 #define GCD_PAIR_COUNT  3
1975
1976 static const int gcd_pairs[GCD_PAIR_COUNT][3] =
1977 {
1978     { 693, 609, 21 },
1979     { 1764, 868, 28 },
1980     { 768454923, 542167814, 1 }
1981 };
1982
1983 /*
1984  * Checkup routine
1985  */
1986 int mpi_self_test( int verbose )
1987 {
1988     int ret, i;
1989     mpi A, E, N, X, Y, U, V;
1990
1991     mpi_init( &A ); mpi_init( &E ); mpi_init( &N ); mpi_init( &X );
1992     mpi_init( &Y ); mpi_init( &U ); mpi_init( &V );
1993
1994     MPI_CHK( mpi_read_string( &A, 16,
1995         "EFE021C2645FD1DC586E69184AF4A31E" \
1996         "D5F53E93B5F123FA41680867BA110131" \
1997         "944FE7952E2517337780CB0DB80E61AA" \
1998         "E7C8DDC6C5C6AADEB34EB38A2F40D5E6" ) );
1999
2000     MPI_CHK( mpi_read_string( &E, 16,
2001         "B2E7EFD37075B9F03FF989C7C5051C20" \
2002         "34D2A323810251127E7BF8625A4F49A5" \
2003         "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
2004         "5B5C25763222FEFCCFC38B832366C29E" ) );
2005
2006     MPI_CHK( mpi_read_string( &N, 16,
2007         "0066A198186C18C10B2F5ED9B522752A" \
2008         "9830B69916E535C8F047518A889A43A5" \
2009         "94B6BED27A168D31D4A52F88925AA8F5" ) );
2010
2011     MPI_CHK( mpi_mul_mpi( &X, &A, &N ) );
2012
2013     MPI_CHK( mpi_read_string( &U, 16,
2014         "602AB7ECA597A3D6B56FF9829A5E8B85" \
2015         "9E857EA95A03512E2BAE7391688D264A" \
2016         "A5663B0341DB9CCFD2C4C5F421FEC814" \
2017         "8001B72E848A38CAE1C65F78E56ABDEF" \
2018         "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2019         "ECF677152EF804370C1A305CAF3B5BF1" \
2020         "30879B56C61DE584A0F53A2447A51E" ) );
2021
2022     if( verbose != 0 )
2023         printf( "  MPI test #1 (mul_mpi): " );
2024
2025     if( mpi_cmp_mpi( &X, &U ) != 0 )
2026     {
2027         if( verbose != 0 )
2028             printf( "failed\n" );
2029
2030         return( 1 );
2031     }
2032
2033     if( verbose != 0 )
2034         printf( "passed\n" );
2035
2036     MPI_CHK( mpi_div_mpi( &X, &Y, &A, &N ) );
2037
2038     MPI_CHK( mpi_read_string( &U, 16,
2039         "256567336059E52CAE22925474705F39A94" ) );
2040
2041     MPI_CHK( mpi_read_string( &V, 16,
2042         "6613F26162223DF488E9CD48CC132C7A" \
2043         "0AC93C701B001B092E4E5B9F73BCD27B" \
2044         "9EE50D0657C77F374E903CDFA4C642" ) );
2045
2046     if( verbose != 0 )
2047         printf( "  MPI test #2 (div_mpi): " );
2048
2049     if( mpi_cmp_mpi( &X, &U ) != 0 ||
2050         mpi_cmp_mpi( &Y, &V ) != 0 )
2051     {
2052         if( verbose != 0 )
2053             printf( "failed\n" );
2054
2055         return( 1 );
2056     }
2057
2058     if( verbose != 0 )
2059         printf( "passed\n" );
2060
2061     MPI_CHK( mpi_exp_mod( &X, &A, &E, &N, NULL ) );
2062
2063     MPI_CHK( mpi_read_string( &U, 16,
2064         "36E139AEA55215609D2816998ED020BB" \
2065         "BD96C37890F65171D948E9BC7CBAA4D9" \
2066         "325D24D6A3C12710F10A09FA08AB87" ) );
2067
2068     if( verbose != 0 )
2069         printf( "  MPI test #3 (exp_mod): " );
2070
2071     if( mpi_cmp_mpi( &X, &U ) != 0 )
2072     {
2073         if( verbose != 0 )
2074             printf( "failed\n" );
2075
2076         return( 1 );
2077     }
2078
2079     if( verbose != 0 )
2080         printf( "passed\n" );
2081
2082 #if defined(POLARSSL_GENPRIME)
2083     MPI_CHK( mpi_inv_mod( &X, &A, &N ) );
2084
2085     MPI_CHK( mpi_read_string( &U, 16,
2086         "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2087         "C3DBA76456363A10869622EAC2DD84EC" \
2088         "C5B8A74DAC4D09E03B5E0BE779F2DF61" ) );
2089
2090     if( verbose != 0 )
2091         printf( "  MPI test #4 (inv_mod): " );
2092
2093     if( mpi_cmp_mpi( &X, &U ) != 0 )
2094     {
2095         if( verbose != 0 )
2096             printf( "failed\n" );
2097
2098         return( 1 );
2099     }
2100
2101     if( verbose != 0 )
2102         printf( "passed\n" );
2103 #endif
2104
2105     if( verbose != 0 )
2106         printf( "  MPI test #5 (simple gcd): " );
2107
2108     for ( i = 0; i < GCD_PAIR_COUNT; i++)
2109     {
2110         MPI_CHK( mpi_lset( &X, gcd_pairs[i][0] ) );
2111         MPI_CHK( mpi_lset( &Y, gcd_pairs[i][1] ) );
2112
2113             MPI_CHK( mpi_gcd( &A, &X, &Y ) );
2114
2115             if( mpi_cmp_int( &A, gcd_pairs[i][2] ) != 0 )
2116             {
2117                     if( verbose != 0 )
2118                             printf( "failed at %d\n", i );
2119
2120                     return( 1 );
2121             }
2122     }
2123
2124     if( verbose != 0 )
2125         printf( "passed\n" );
2126
2127 cleanup:
2128
2129     if( ret != 0 && verbose != 0 )
2130         printf( "Unexpected error, return code = %08X\n", ret );
2131
2132     mpi_free( &A ); mpi_free( &E ); mpi_free( &N ); mpi_free( &X );
2133     mpi_free( &Y ); mpi_free( &U ); mpi_free( &V );
2134
2135     if( verbose != 0 )
2136         printf( "\n" );
2137
2138     return( ret );
2139 }
2140
2141 #endif
2142
2143 #endif