freeswan/pluto/primegen.c

   1 /* primegen.c - prime number generator
   2  *      Copyright (C) 1998 Free Software Foundation, Inc.
   3  *
   4  * This file is part of GnuPG.
   5  *
   6  * GnuPG is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 2 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * GnuPG is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
  19  *
  20  * ***********************************************************************
  21  * The algorithm used to generate practically save primes is due to
  22  * Lim and Lee as described in the CRYPTO '97 proceedings (ISBN3540633847)
  23  * page 260.
  24  */
  25
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28 #include <string.h>
  29
  30 #ifdef PLUTO
  31 #include <gmp.h>
  32 #include <freeswan.h>
  33 #include "constants.h"
  34 #include "defs.h"
  35 #include "log.h"
  36 #include "rnd.h"
  37 #include "gcryptfix.h"
  38 #else /*! PLUTO */
  39 /* #include <assert.h> */
  40 /* #include <config.h> */
  41 /* #include "util.h" */
  42 /* #include "mpi.h" */
  43 /* #include "cipher.h" */
  44 #endif /* !PLUTO */
  45
  46 static int no_of_small_prime_numbers;
  47 static MPI gen_prime( unsigned  nbits, int mode, int randomlevel );
  48 static int check_prime( MPI prime, MPI val_2 );
  49 static int is_prime( MPI n, unsigned steps, int *count );
  50 static void m_out_of_n( char *array, int m, int n );
  51
  52
  53 static void
  54 progress( int c )
  55 {
  56     fputc( c, stderr );
  57 }
  58
  59
  60 /****************
  61  * Generate a prime number (stored in secure memory)
  62  */
  63 MPI
  64 generate_secret_prime( unsigned  nbits )
  65 {
  66     MPI prime;
  67
  68     prime = gen_prime( nbits, 1, 2 );
  69     progress('\n');
  70     return prime;
  71 }
  72
  73 MPI
  74 generate_public_prime( unsigned  nbits )
  75 {
  76     MPI prime;
  77
  78     prime = gen_prime( nbits, 0, 2 );
  79     progress('\n');
  80     return prime;
  81 }
  82
  83
  84 /****************
  85  * We do not need to use the strongest RNG because we gain no extra
  86  * security from it - The prime number is public and we could also
  87  * offer the factors for those who are willing to check that it is
  88  * indeed a strong prime.
  89  *
  90  * mode 0: Standard
  91  *      1: Make sure that at least one factor is of size qbits.
  92  */
  93 MPI
  94 generate_elg_prime( int mode, unsigned pbits, unsigned qbits,
  95                     MPI g, MPI **ret_factors )
  96 {
  97     int n;  /* number of factors */
  98     int m;  /* number of primes in pool */
  99     unsigned fbits; /* length of prime factors */
 100     MPI *factors; /* current factors */
 101     MPI *pool;  /* pool of primes */
 102     MPI q;      /* first prime factor (variable)*/
 103     MPI prime;  /* prime test value */
 104     MPI q_factor; /* used for mode 1 */
 105     byte *perms = NULL;
 106     int i, j;
 107     int count1, count2;
 108     unsigned nprime;
 109     unsigned req_qbits = qbits; /* the requested q bits size */
 110     MPI val_2  = mpi_alloc_set_ui( 2 );
 111
 112     /* find number of needed prime factors */
 113     for(n=1; (pbits - qbits - 1) / n  >= qbits; n++ )
 114         ;
 115     n--;
 116     if( !n || (mode==1 && n < 2) )
 117         log_fatal("can't gen prime with pbits=%u qbits=%u\n", pbits, qbits );
 118     if( mode == 1 ) {
 119         n--;
 120         fbits = (pbits - 2*req_qbits -1) / n;
 121         qbits =  pbits - req_qbits - n*fbits;
 122     }
 123     else {
 124         fbits = (pbits - req_qbits -1) / n;
 125         qbits = pbits - n*fbits;
 126     }
 127     if( DBG_CIPHER )
 128         log_debug("gen prime: pbits=%u qbits=%u fbits=%u/%u n=%d\n",
 129                     pbits, req_qbits, qbits, fbits, n  );
 130     prime = mpi_alloc( (pbits + BITS_PER_MPI_LIMB - 1) /  BITS_PER_MPI_LIMB );
 131     q = gen_prime( qbits, 0, 1 );
 132     q_factor = mode==1? gen_prime( req_qbits, 0, 1 ) : NULL;
 133
 134     /* allocate an array to hold the factors + 2 for later usage */
 135 #ifdef PLUTO
 136     m_alloc_ptrs_clear(factors, n+2);
 137 #else
 138     factors = m_alloc_clear( (n+2) * sizeof *factors );
 139 #endif
 140
 141     /* make a pool of 3n+5 primes (this is an arbitrary value) */
 142     m = n*3+5;
 143     if( mode == 1 )
 144         m += 5; /* need some more for DSA */
 145     if( m < 25 )
 146         m = 25;
 147 #ifdef PLUTO
 148     m_alloc_ptrs_clear(pool, m);
 149 #else
 150     pool = m_alloc_clear( m * sizeof *pool );
 151 #endif
 152
 153     /* permutate over the pool of primes */
 154     count1=count2=0;
 155     do {
 156       next_try:
 157         if( !perms ) {
 158             /* allocate new primes */
 159             for(i=0; i < m; i++ ) {
 160                 mpi_free(pool[i]);
 161                 pool[i] = NULL;
 162             }
 163             /* init m_out_of_n() */
 164 #ifdef PLUTO
 165             perms = alloc_bytes( m, "perms" );
 166 #else
 167             perms = m_alloc_clear( m );
 168 #endif
 169             for(i=0; i < n; i++ ) {
 170                 perms[i] = 1;
 171                 pool[i] = gen_prime( fbits, 0, 1 );
 172                 factors[i] = pool[i];
 173             }
 174         }
 175         else {
 176             m_out_of_n( perms, n, m );
 177             for(i=j=0; i < m && j < n ; i++ )
 178                 if( perms[i] ) {
 179                     if( !pool[i] )
 180                         pool[i] = gen_prime( fbits, 0, 1 );
 181                     factors[j++] = pool[i];
 182                 }
 183             if( i == n ) {
 184                 m_free(perms); perms = NULL;
 185                 progress('!');
 186                 goto next_try;  /* allocate new primes */
 187             }
 188         }
 189
 190         mpi_set( prime, q );
 191         mpi_mul_ui( prime, prime, 2 );
 192         if( mode == 1 )
 193             mpi_mul( prime, prime, q_factor );
 194         for(i=0; i < n; i++ )
 195             mpi_mul( prime, prime, factors[i] );
 196         mpi_add_ui( prime, prime, 1 );
 197         nprime = mpi_get_nbits(prime);
 198         if( nprime < pbits ) {
 199             if( ++count1 > 20 ) {
 200                 count1 = 0;
 201                 qbits++;
 202                 progress('>');
 203                 q = gen_prime( qbits, 0, 1 );
 204                 goto next_try;
 205             }
 206         }
 207         else
 208             count1 = 0;
 209         if( nprime > pbits ) {
 210             if( ++count2 > 20 ) {
 211                 count2 = 0;
 212                 qbits--;
 213                 progress('<');
 214                 q = gen_prime( qbits, 0, 1 );
 215                 goto next_try;
 216             }
 217         }
 218         else
 219             count2 = 0;
 220     } while( !(nprime == pbits && check_prime( prime, val_2 )) );
 221
 222     if( DBG_CIPHER ) {
 223         progress('\n');
 224         log_mpidump( "prime    : ", prime );
 225         log_mpidump( "factor  q: ", q );
 226         if( mode == 1 )
 227             log_mpidump( "factor q0: ", q_factor );
 228         for(i=0; i < n; i++ )
 229             log_mpidump( "factor pi: ", factors[i] );
 230         log_debug("bit sizes: prime=%u, q=%u", mpi_get_nbits(prime), mpi_get_nbits(q) );
 231         if( mode == 1 )
 232             fprintf(stderr, ", q0=%u", mpi_get_nbits(q_factor) );
 233         for(i=0; i < n; i++ )
 234             fprintf(stderr, ", p%d=%u", i, mpi_get_nbits(factors[i]) );
 235         progress('\n');
 236     }
 237
 238     if( ret_factors ) { /* caller wants the factors */
 239 #ifdef PLUTO
 240         m_alloc_ptrs_clear(*ret_factors, n+2);
 241 #else
 242         *ret_factors = m_alloc_clear( (n+2) * sizeof **ret_factors);
 243 #endif
 244         if( mode == 1 ) {
 245             i = 0;
 246             (*ret_factors)[i++] = mpi_copy( q_factor );
 247             for(; i <= n; i++ )
 248                 (*ret_factors)[i] = mpi_copy( factors[i] );
 249         }
 250         else {
 251             for(; i < n; i++ )
 252                 (*ret_factors)[i] = mpi_copy( factors[i] );
 253         }
 254     }
 255
 256     if( g ) { /* create a generator (start with 3)*/
 257         MPI tmp   = mpi_alloc( mpi_get_nlimbs(prime) );
 258         MPI b     = mpi_alloc( mpi_get_nlimbs(prime) );
 259         MPI pmin1 = mpi_alloc( mpi_get_nlimbs(prime) );
 260
 261         if( mode == 1 )
 262             BUG(); /* not yet implemented */
 263         factors[n] = q;
 264         factors[n+1] = mpi_alloc_set_ui(2);
 265         mpi_sub_ui( pmin1, prime, 1 );
 266         mpi_set_ui(g,2);
 267         do {
 268             mpi_add_ui(g, g, 1);
 269             if( DBG_CIPHER ) {
 270 #ifdef PLUTO
 271                 log_mpidump("checking g: ", g);
 272 #else
 273                 log_debug("checking g: ");
 274                 mpi_print( stderr, g, 1 );
 275 #endif
 276             }
 277             else
 278                 progress('^');
 279             for(i=0; i < n+2; i++ ) {
 280                 /*fputc('~', stderr);*/
 281                 mpi_fdiv_q(tmp, pmin1, factors[i] );
 282                 /* (no mpi_pow(), but it is okay to use this with mod prime) */
 283                 mpi_powm(b, g, tmp, prime );
 284                 if( !mpi_cmp_ui(b, 1) )
 285                     break;
 286             }
 287             if( DBG_CIPHER )
 288                 progress('\n');
 289         } while( i < n+2 );
 290         mpi_free(factors[n+1]);
 291         mpi_free(tmp);
 292         mpi_free(b);
 293         mpi_free(pmin1);
 294     }
 295     if( !DBG_CIPHER )
 296         progress('\n');
 297
 298     m_free( factors );  /* (factors are shallow copies) */
 299     for(i=0; i < m; i++ )
 300         mpi_free( pool[i] );
 301     m_free( pool );
 302     m_free(perms);
 303     mpi_free(val_2);
 304     return prime;
 305 }
 306
 307
 308
 309 static MPI
 310 gen_prime( unsigned  nbits, int secret, int randomlevel )
 311 {
 312     unsigned  nlimbs;
 313     MPI prime, ptest, pminus1, val_2, val_3, result;
 314     int i;
 315     unsigned x, step;
 316     unsigned count1, count2;
 317     int *mods;
 318
 319     if( 0 && DBG_CIPHER )
 320         log_debug("generate a prime of %u bits ", nbits );
 321
 322     if( !no_of_small_prime_numbers ) {
 323         for(i=0; small_prime_numbers[i]; i++ )
 324             no_of_small_prime_numbers++;
 325     }
 326     mods = m_alloc( no_of_small_prime_numbers * sizeof *mods );
 327     /* make nbits fit into MPI implementation */
 328     nlimbs = (nbits + BITS_PER_MPI_LIMB - 1) /  BITS_PER_MPI_LIMB;
 329     val_2  = mpi_alloc_set_ui( 2 );
 330     val_3 = mpi_alloc_set_ui( 3);
 331     prime  = secret? mpi_alloc_secure( nlimbs ): mpi_alloc( nlimbs );
 332     result = mpi_alloc_like( prime );
 333     pminus1= mpi_alloc_like( prime );
 334     ptest  = mpi_alloc_like( prime );
 335     count1 = count2 = 0;
 336     for(;;) {  /* try forvever */
 337         int dotcount=0;
 338
 339         /* generate a random number */
 340         {   char *p = get_random_bits( nbits, randomlevel, secret );
 341             mpi_set_buffer( prime, p, (nbits+7)/8, 0 );
 342             m_free(p);
 343         }
 344
 345         /* set high order bit to 1, set low order bit to 1 */
 346         mpi_set_highbit( prime, nbits-1 );
 347         mpi_set_bit( prime, 0 );
 348
 349         /* calculate all remainders */
 350         for(i=0; (x = small_prime_numbers[i]); i++ )
 351             mods[i] = mpi_fdiv_r_ui(NULL, prime, x);
 352
 353         /* now try some primes starting with prime */
 354         for(step=0; step < 20000; step += 2 ) {
 355             /* check against all the small primes we have in mods */
 356             count1++;
 357             for(i=0; (x = small_prime_numbers[i]); i++ ) {
 358                 while( mods[i] + step >= x )
 359                     mods[i] -= x;
 360                 if( !(mods[i] + step) )
 361                     break;
 362             }
 363             if( x )
 364                 continue;   /* found a multiple of an already known prime */
 365
 366             mpi_add_ui( ptest, prime, step );
 367
 368             /* do a faster Fermat test */
 369             count2++;
 370             mpi_sub_ui( pminus1, ptest, 1);
 371             mpi_powm( result, val_2, pminus1, ptest );
 372             if( !mpi_cmp_ui( result, 1 ) ) { /* not composite */
 373                 /* perform stronger tests */
 374                 if( is_prime(ptest, 5, &count2 ) ) {
 375                     if( !mpi_test_bit( ptest, nbits-1 ) ) {
 376                         progress('\n');
 377                         log_debug("overflow in prime generation\n");
 378                         break; /* step loop, continue with a new prime */
 379                     }
 380
 381                     mpi_free(val_2);
 382                     mpi_free(val_3);
 383                     mpi_free(result);
 384                     mpi_free(pminus1);
 385                     mpi_free(prime);
 386                     m_free(mods);
 387                     return ptest;
 388                 }
 389             }
 390             if( ++dotcount == 10 ) {
 391                 progress('.');
 392                 dotcount = 0;
 393             }
 394         }
 395         progress(':'); /* restart with a new random value */
 396     }
 397 }
 398
 399 /****************
 400  * Returns: true if this may be a prime
 401  */
 402 static int
 403 check_prime( MPI prime, MPI val_2 )
 404 {
 405     int i;
 406     unsigned x;
 407     int count=0;
 408
 409     /* check against small primes */
 410     for(i=0; (x = small_prime_numbers[i]); i++ ) {
 411         if( mpi_divisible_ui( prime, x ) )
 412             return 0;
 413     }
 414
 415     /* a quick fermat test */
 416     {
 417         MPI result = mpi_alloc_like( prime );
 418         MPI pminus1 = mpi_alloc_like( prime );
 419         mpi_sub_ui( pminus1, prime, 1);
 420         mpi_powm( result, val_2, pminus1, prime );
 421         mpi_free( pminus1 );
 422         if( mpi_cmp_ui( result, 1 ) ) { /* if composite */
 423             mpi_free( result );
 424             progress('.');
 425             return 0;
 426         }
 427         mpi_free( result );
 428     }
 429
 430     /* perform stronger tests */
 431     if( is_prime(prime, 5, &count ) )
 432         return 1; /* is probably a prime */
 433     progress('.');
 434     return 0;
 435 }
 436
 437
 438 /****************
 439  * Return true if n is probably a prime
 440  */
 441 static int
 442 is_prime( MPI n, unsigned steps, int *count )
 443 {
 444     MPI x = mpi_alloc( mpi_get_nlimbs( n ) );
 445     MPI y = mpi_alloc( mpi_get_nlimbs( n ) );
 446     MPI z = mpi_alloc( mpi_get_nlimbs( n ) );
 447     MPI nminus1 = mpi_alloc( mpi_get_nlimbs( n ) );
 448     MPI a2 = mpi_alloc_set_ui( 2 );
 449     MPI q;
 450     unsigned i, j, k;
 451     int rc = 0;
 452     unsigned nbits = mpi_get_nbits( n );
 453
 454     mpi_sub_ui( nminus1, n, 1 );
 455
 456     /* find q and k, so that n = 1 + 2^k * q */
 457     q = mpi_copy( nminus1 );
 458     k = mpi_trailing_zeros( q );
 459     mpi_tdiv_q_2exp(q, q, k);
 460
 461     for(i=0 ; i < steps; i++ ) {
 462         ++*count;
 463         if( !i ) {
 464             mpi_set_ui( x, 2 );
 465         }
 466         else {
 467             /*mpi_set_bytes( x, nbits-1, get_random_byte, 0 );*/
 468             {   char *p = get_random_bits( nbits, 0, 0 );
 469                 mpi_set_buffer( x, p, (nbits+7)/8, 0 );
 470                 m_free(p);
 471             }
 472             /* make sure that the number is smaller than the prime
 473              * and keep the randomness of the high bit */
 474             if( mpi_test_bit( x, nbits-2 ) ) {
 475                 mpi_set_highbit( x, nbits-2 ); /* clear all higher bits */
 476             }
 477             else {
 478                 mpi_set_highbit( x, nbits-2 );
 479                 mpi_clear_bit( x, nbits-2 );
 480             }
 481             assert( mpi_cmp( x, nminus1 ) < 0 && mpi_cmp_ui( x, 1 ) > 0 );
 482         }
 483         mpi_powm( y, x, q, n);
 484         if( mpi_cmp_ui(y, 1) && mpi_cmp( y, nminus1 ) ) {
 485             for( j=1; j < k && mpi_cmp( y, nminus1 ); j++ ) {
 486                 mpi_powm(y, y, a2, n);
 487                 if( !mpi_cmp_ui( y, 1 ) )
 488                     goto leave; /* not a prime */
 489             }
 490             if( mpi_cmp( y, nminus1 ) )
 491                 goto leave; /* not a prime */
 492         }
 493         progress('+');
 494     }
 495     rc = 1; /* may be a prime */
 496
 497   leave:
 498     mpi_free( x );
 499     mpi_free( y );
 500     mpi_free( z );
 501     mpi_free( nminus1 );
 502     mpi_free( q );
 503
 504     return rc;
 505 }
 506
 507
 508 static void
 509 m_out_of_n( char *array, int m, int n )
 510 {
 511     int i=0, i1=0, j=0, jp=0,  j1=0, k1=0, k2=0;
 512
 513     if( !m || m >= n )
 514         return;
 515
 516     if( m == 1 ) { /* special case */
 517         for(i=0; i < n; i++ )
 518             if( array[i] ) {
 519                 array[i++] = 0;
 520                 if( i >= n )
 521                     i = 0;
 522                 array[i] = 1;
 523                 return;
 524             }
 525         BUG();
 526     }
 527
 528     for(j=1; j < n; j++ ) {
 529         if( array[n-1] == array[n-j-1] )
 530             continue;
 531         j1 = j;
 532         break;
 533     }
 534
 535     if( m & 1 ) { /* m is odd */
 536         if( array[n-1] ) {
 537             if( j1 & 1 ) {
 538                 k1 = n - j1;
 539                 k2 = k1+2;
 540                 if( k2 > n )
 541                     k2 = n;
 542                 goto leave;
 543             }
 544             goto scan;
 545         }
 546         k2 = n - j1 - 1;
 547         if( k2 == 0 ) {
 548             k1 = i;
 549             k2 = n - j1;
 550         }
 551         else if( array[k2] && array[k2-1] )
 552             k1 = n;
 553         else
 554             k1 = k2 + 1;
 555     }
 556     else { /* m is even */
 557         if( !array[n-1] ) {
 558             k1 = n - j1;
 559             k2 = k1 + 1;
 560             goto leave;
 561         }
 562
 563         if( !(j1 & 1) ) {
 564             k1 = n - j1;
 565             k2 = k1+2;
 566             if( k2 > n )
 567                 k2 = n;
 568             goto leave;
 569         }
 570       scan:
 571         jp = n - j1 - 1;
 572         for(i=1; i <= jp; i++ ) {
 573             i1 = jp + 2 - i;
 574             if( array[i1-1]  ) {
 575                 if( array[i1-2] ) {
 576                     k1 = i1 - 1;
 577                     k2 = n - j1;
 578                 }
 579                 else {
 580                     k1 = i1 - 1;
 581                     k2 = n + 1 - j1;
 582                 }
 583                 goto leave;
 584             }
 585         }
 586         k1 = 1;
 587         k2 = n + 1 - m;
 588     }
 589   leave:
 590     array[k1-1] = !array[k1-1];
 591     array[k2-1] = !array[k2-1];
 592 }
 593