package/utils/px5g/src/library/bignum.c

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