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