src/xyssl/library/bignum.c

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