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