src/src/pdkim/bignum.c

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