src/src/pdkim/bignum.c

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