2003-07-14 Moritz Schulte <moritz@g10code.com>
[libgcrypt.git] / cipher / primegen.c
1 /* primegen.c - prime number generator
2  * Copyright (C) 1998, 2000, 2001, 2002, 2003 Free Software Foundation, Inc.
3  *
4  * This file is part of Libgcrypt.
5  *
6  * Libgcrypt is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser general Public License as
8  * published by the Free Software Foundation; either version 2.1 of
9  * the License, or (at your option) any later version.
10  *
11  * Libgcrypt 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 Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License 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 "g10lib.h"
32 #include "mpi.h"
33 #include "cipher.h"
34
35 static gcry_mpi_t gen_prime (unsigned int nbits, int secret, int randomlevel, 
36                       int (*extra_check)(void *, gcry_mpi_t), void *extra_check_arg);
37 static int check_prime( gcry_mpi_t prime, gcry_mpi_t val_2 );
38 static int is_prime( gcry_mpi_t n, int steps, int *count );
39 static void m_out_of_n( char *array, int m, int n );
40
41 static void (*progress_cb) (void *,const char*,int,int, int );
42 static void *progress_cb_data;
43
44 /* Note: 2 is not included because it can be tested more easily by
45    looking at bit 0. The last entry in this list is marked by a zero */
46 static ushort small_prime_numbers[] = {
47     3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
48     47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101,
49     103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
50     157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
51     211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
52     269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
53     331, 337, 347, 349, 353, 359, 367, 373, 379, 383,
54     389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
55     449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
56     509, 521, 523, 541, 547, 557, 563, 569, 571, 577,
57     587, 593, 599, 601, 607, 613, 617, 619, 631, 641,
58     643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
59     709, 719, 727, 733, 739, 743, 751, 757, 761, 769,
60     773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
61     853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
62     919, 929, 937, 941, 947, 953, 967, 971, 977, 983,
63     991, 997, 1009, 1013, 1019, 1021, 1031, 1033,
64     1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091,
65     1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
66     1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213,
67     1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,
68     1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307,
69     1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
70     1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
71     1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
72     1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
73     1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609,
74     1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667,
75     1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
76     1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
77     1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871,
78     1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931,
79     1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
80     1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
81     2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111,
82     2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
83     2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243,
84     2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
85     2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
86     2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411,
87     2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473,
88     2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551,
89     2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633,
90     2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
91     2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729,
92     2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791,
93     2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
94     2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917,
95     2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
96     3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061,
97     3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137,
98     3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
99     3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271,
100     3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
101     3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391,
102     3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467,
103     3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533,
104     3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,
105     3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
106     3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,
107     3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779,
108     3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851,
109     3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917,
110     3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
111     4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049,
112     4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111,
113     4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177,
114     4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243,
115     4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
116     4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391,
117     4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457,
118     4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,
119     4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597,
120     4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
121     4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729,
122     4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799,
123     4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889,
124     4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951,
125     4957, 4967, 4969, 4973, 4987, 4993, 4999,
126     0
127 };
128 static int no_of_small_prime_numbers = DIM (small_prime_numbers) - 1;
129
130 void
131 _gcry_register_primegen_progress ( void (*cb)(void *,const char*,int,int,int), void *cb_data )
132 {
133     progress_cb = cb;
134     progress_cb_data = cb_data;
135 }
136
137
138 static void
139 progress( int c )
140 {
141   if ( progress_cb )
142     progress_cb ( progress_cb_data, "primegen", c, 0, 0 );
143 }
144
145
146 /****************
147  * Generate a prime number (stored in secure memory)
148  */
149 gcry_mpi_t
150 _gcry_generate_secret_prime (unsigned int nbits,
151                              int (*extra_check)(void*, gcry_mpi_t),
152                              void *extra_check_arg)
153 {
154     gcry_mpi_t prime;
155
156     prime = gen_prime( nbits, 1, 2, extra_check, extra_check_arg);
157     progress('\n');
158     return prime;
159 }
160
161 gcry_mpi_t
162 _gcry_generate_public_prime( unsigned int nbits,
163                              int (*extra_check)(void*, gcry_mpi_t),
164                              void *extra_check_arg)
165 {
166     gcry_mpi_t prime;
167
168     prime = gen_prime( nbits, 0, 2, extra_check, extra_check_arg );
169     progress('\n');
170     return prime;
171 }
172
173
174 /****************
175  * We do not need to use the strongest RNG because we gain no extra
176  * security from it - The prime number is public and we could also
177  * offer the factors for those who are willing to check that it is
178  * indeed a strong prime.
179  *
180  * mode 0: Standard
181  *      1: Make sure that at least one factor is of size qbits.
182  */
183 gcry_mpi_t
184 _gcry_generate_elg_prime( int mode, unsigned pbits, unsigned qbits,
185                     gcry_mpi_t g, gcry_mpi_t **ret_factors )
186 {
187     int n;  /* number of factors */
188     int m;  /* number of primes in pool */
189     unsigned fbits; /* length of prime factors */
190     gcry_mpi_t *factors; /* current factors */
191     gcry_mpi_t *pool;   /* pool of primes */
192     gcry_mpi_t q;       /* first prime factor (variable)*/
193     gcry_mpi_t prime;   /* prime test value */
194     gcry_mpi_t q_factor; /* used for mode 1 */
195     byte *perms = NULL;
196     int i, j;
197     int count1, count2;
198     unsigned nprime;
199     unsigned req_qbits = qbits; /* the requested q bits size */
200     gcry_mpi_t val_2  = mpi_alloc_set_ui( 2 );
201
202     /* find number of needed prime factors */
203     for(n=1; (pbits - qbits - 1) / n  >= qbits; n++ )
204         ;
205     n--;
206     if( !n || (mode==1 && n < 2) )
207         log_fatal("can't gen prime with pbits=%u qbits=%u\n", pbits, qbits );
208     if( mode == 1 ) {
209         n--;
210         fbits = (pbits - 2*req_qbits -1) / n;
211         qbits =  pbits - req_qbits - n*fbits;
212     }
213     else {
214         fbits = (pbits - req_qbits -1) / n;
215         qbits = pbits - n*fbits;
216     }
217     if( DBG_CIPHER )
218         log_debug("gen prime: pbits=%u qbits=%u fbits=%u/%u n=%d\n",
219                     pbits, req_qbits, qbits, fbits, n  );
220     prime = gcry_mpi_new ( pbits  );
221     q = gen_prime( qbits, 0, 0, NULL, NULL );
222     q_factor = mode==1? gen_prime( req_qbits, 0, 0, NULL, NULL ) : NULL;
223
224     /* allocate an array to hold the factors + 2 for later usage */
225     factors = gcry_xcalloc( n+2, sizeof *factors );
226
227     /* make a pool of 3n+5 primes (this is an arbitrary value) */
228     m = n*3+5;
229     if( mode == 1 )
230         m += 5; /* need some more for DSA */
231     if( m < 25 )
232         m = 25;
233     pool = gcry_xcalloc( m , sizeof *pool );
234
235     /* permutate over the pool of primes */
236     count1=count2=0;
237     do {
238       next_try:
239         if( !perms ) {
240             /* allocate new primes */
241             for(i=0; i < m; i++ ) {
242                 mpi_free(pool[i]);
243                 pool[i] = NULL;
244             }
245             /* init m_out_of_n() */
246             perms = gcry_xcalloc( 1, m );
247             for(i=0; i < n; i++ ) {
248                 perms[i] = 1;
249                 pool[i] = gen_prime( fbits, 0, 1, NULL, NULL );
250                 factors[i] = pool[i];
251             }
252         }
253         else {
254             m_out_of_n( perms, n, m );
255             for(i=j=0; i < m && j < n ; i++ )
256                 if( perms[i] ) {
257                     if( !pool[i] )
258                         pool[i] = gen_prime( fbits, 0, 1, NULL, NULL );
259                     factors[j++] = pool[i];
260                 }
261             if( i == n ) {
262                 gcry_free(perms); perms = NULL;
263                 progress('!');
264                 goto next_try;  /* allocate new primes */
265             }
266         }
267
268         mpi_set( prime, q );
269         mpi_mul_ui( prime, prime, 2 );
270         if( mode == 1 )
271             mpi_mul( prime, prime, q_factor );
272         for(i=0; i < n; i++ )
273             mpi_mul( prime, prime, factors[i] );
274         mpi_add_ui( prime, prime, 1 );
275         nprime = mpi_get_nbits(prime);
276         if( nprime < pbits ) {
277             if( ++count1 > 20 ) {
278                 count1 = 0;
279                 qbits++;
280                 progress('>');
281                 mpi_free (q);
282                 q = gen_prime( qbits, 0, 0, NULL, NULL );
283                 goto next_try;
284             }
285         }
286         else
287             count1 = 0;
288         if( nprime > pbits ) {
289             if( ++count2 > 20 ) {
290                 count2 = 0;
291                 qbits--;
292                 progress('<');
293                 mpi_free (q);
294                 q = gen_prime( qbits, 0, 0, NULL, NULL );
295                 goto next_try;
296             }
297         }
298         else
299             count2 = 0;
300     } while( !(nprime == pbits && check_prime( prime, val_2 )) );
301
302     if( DBG_CIPHER ) {
303         progress('\n');
304         log_mpidump( "prime    : ", prime );
305         log_mpidump( "factor  q: ", q );
306         if( mode == 1 )
307             log_mpidump( "factor q0: ", q_factor );
308         for(i=0; i < n; i++ )
309             log_mpidump( "factor pi: ", factors[i] );
310         log_debug("bit sizes: prime=%u, q=%u", mpi_get_nbits(prime), mpi_get_nbits(q) );
311         if( mode == 1 )
312             log_debug (", q0=%u", mpi_get_nbits(q_factor) );
313         for(i=0; i < n; i++ )
314             log_debug (", p%d=%u", i, mpi_get_nbits(factors[i]) );
315         progress('\n');
316     }
317
318     if( ret_factors ) { /* caller wants the factors */
319         *ret_factors = gcry_xcalloc( n+2 , sizeof **ret_factors);
320         i = 0;
321         if( mode == 1 ) {
322             (*ret_factors)[i++] = mpi_copy( q_factor );
323             for(; i <= n; i++ )
324                 (*ret_factors)[i] = mpi_copy( factors[i] );
325         }
326         else {
327             for(; i < n; i++ )
328                 (*ret_factors)[i] = mpi_copy( factors[i] );
329         }
330     }
331
332     if( g ) { /* create a generator (start with 3)*/
333         gcry_mpi_t tmp   = mpi_alloc( mpi_get_nlimbs(prime) );
334         gcry_mpi_t b      = mpi_alloc( mpi_get_nlimbs(prime) );
335         gcry_mpi_t pmin1 = mpi_alloc( mpi_get_nlimbs(prime) );
336
337         if( mode == 1 )
338             BUG(); /* not yet implemented */
339         factors[n] = q;
340         factors[n+1] = mpi_alloc_set_ui(2);
341         mpi_sub_ui( pmin1, prime, 1 );
342         mpi_set_ui(g,2);
343         do {
344             mpi_add_ui(g, g, 1);
345             if( DBG_CIPHER ) {
346                 log_debug("checking g: ");
347                 /*mpi_print( stderr, g, 1 );*/
348 #if __GNUC__ >= 2
349 #warning we need an internal mpi_print for debugging
350 #endif
351             }
352             else
353                 progress('^');
354             for(i=0; i < n+2; i++ ) {
355                 /*fputc('~', stderr);*/
356                 mpi_fdiv_q(tmp, pmin1, factors[i] );
357                 /* (no mpi_pow(), but it is okay to use this with mod prime) */
358                 gcry_mpi_powm(b, g, tmp, prime );
359                 if( !mpi_cmp_ui(b, 1) )
360                     break;
361             }
362             if( DBG_CIPHER )
363                 progress('\n');
364         } while( i < n+2 );
365         mpi_free(factors[n+1]);
366         mpi_free(tmp);
367         mpi_free(b);
368         mpi_free(pmin1);
369     }
370     if( !DBG_CIPHER )
371         progress('\n');
372
373     gcry_free( factors );  /* (factors are shallow copies) */
374     for(i=0; i < m; i++ )
375         mpi_free( pool[i] );
376     gcry_free( pool );
377     gcry_free(perms);
378     mpi_free(val_2);
379     mpi_free (q);
380     return prime;
381 }
382
383
384
385 static gcry_mpi_t
386 gen_prime (unsigned int nbits, int secret, int randomlevel, 
387            int (*extra_check)(void *, gcry_mpi_t), void *extra_check_arg)
388 {
389   gcry_mpi_t prime, ptest, pminus1, val_2, val_3, result;
390   int i;
391   unsigned x, step;
392   unsigned count1, count2;
393   int *mods;
394   
395   if( 0 && DBG_CIPHER )
396     log_debug("generate a prime of %u bits ", nbits );
397
398   mods = gcry_xmalloc( no_of_small_prime_numbers * sizeof *mods );
399   /* make nbits fit into gcry_mpi_t implementation */
400   val_2  = mpi_alloc_set_ui( 2 );
401   val_3 = mpi_alloc_set_ui( 3);
402   prime  = secret? gcry_mpi_snew ( nbits ): gcry_mpi_new ( nbits );
403   result = mpi_alloc_like( prime );
404   pminus1= mpi_alloc_like( prime );
405   ptest  = mpi_alloc_like( prime );
406   count1 = count2 = 0;
407   for (;;)
408     {  /* try forvever */
409       int dotcount=0;
410       
411       /* generate a random number */
412       gcry_mpi_randomize( prime, nbits, randomlevel );
413       
414       /* Set high order bit to 1, set low order bit to 0.  If we are
415          generating a secret prime we are most probably doing that
416          for RSA, to make sure that the modulus does have the
417          requested keysize we set the 2 high order bits */
418       mpi_set_highbit (prime, nbits-1);
419       if (secret)
420         mpi_set_bit (prime, nbits-2);
421       mpi_set_bit(prime, 0);
422       
423       /* calculate all remainders */
424       for (i=0; (x = small_prime_numbers[i]); i++ )
425         mods[i] = mpi_fdiv_r_ui(NULL, prime, x);
426       
427       /* now try some primes starting with prime */
428       for(step=0; step < 20000; step += 2 ) 
429         {
430           /* check against all the small primes we have in mods */
431           count1++;
432           for (i=0; (x = small_prime_numbers[i]); i++ ) 
433             {
434               while ( mods[i] + step >= x )
435                 mods[i] -= x;
436               if ( !(mods[i] + step) )
437                 break;
438             }
439           if ( x )
440             continue;   /* found a multiple of an already known prime */
441
442           mpi_add_ui( ptest, prime, step );
443
444           /* do a faster Fermat test */
445           count2++;
446           mpi_sub_ui( pminus1, ptest, 1);
447           gcry_mpi_powm( result, val_2, pminus1, ptest );
448           if ( !mpi_cmp_ui( result, 1 ) )
449             { /* not composite, perform stronger tests */
450                 if (is_prime(ptest, 5, &count2 ))
451                   {
452                     if (!mpi_test_bit( ptest, nbits-1-secret ))
453                       {
454                         progress('\n');
455                         log_debug("overflow in prime generation\n");
456                         break; /* stop loop, continue with a new prime */
457                       }
458
459                     if (extra_check && extra_check (extra_check_arg, ptest))
460                       { /* The extra check told us that this prime is
461                            not of the caller's taste. */
462                         progress ('/');
463                       }
464                     else
465                       { /* got it */
466                         mpi_free(val_2);
467                         mpi_free(val_3);
468                         mpi_free(result);
469                         mpi_free(pminus1);
470                         mpi_free(prime);
471                         gcry_free(mods);
472                         return ptest; 
473                       }
474                   }
475             }
476           if (++dotcount == 10 )
477             {
478               progress('.');
479               dotcount = 0;
480             }
481         }
482       progress(':'); /* restart with a new random value */
483     }
484 }
485
486 /****************
487  * Returns: true if this may be a prime
488  */
489 static int
490 check_prime( gcry_mpi_t prime, gcry_mpi_t val_2 )
491 {
492     int i;
493     unsigned x;
494     int count=0;
495
496     /* check against small primes */
497     for(i=0; (x = small_prime_numbers[i]); i++ ) {
498         if( mpi_divisible_ui( prime, x ) )
499             return 0;
500     }
501
502     /* a quick fermat test */
503     {
504         gcry_mpi_t result = mpi_alloc_like( prime );
505         gcry_mpi_t pminus1 = mpi_alloc_like( prime );
506         mpi_sub_ui( pminus1, prime, 1);
507         gcry_mpi_powm( result, val_2, pminus1, prime );
508         mpi_free( pminus1 );
509         if( mpi_cmp_ui( result, 1 ) ) { /* if composite */
510             mpi_free( result );
511             progress('.');
512             return 0;
513         }
514         mpi_free( result );
515     }
516
517     /* perform stronger tests */
518     if( is_prime(prime, 5, &count ) )
519         return 1; /* is probably a prime */
520     progress('.');
521     return 0;
522 }
523
524
525 /****************
526  * Return true if n is probably a prime
527  */
528 static int
529 is_prime( gcry_mpi_t n, int steps, int *count )
530 {
531     gcry_mpi_t x = mpi_alloc( mpi_get_nlimbs( n ) );
532     gcry_mpi_t y = mpi_alloc( mpi_get_nlimbs( n ) );
533     gcry_mpi_t z = mpi_alloc( mpi_get_nlimbs( n ) );
534     gcry_mpi_t nminus1 = mpi_alloc( mpi_get_nlimbs( n ) );
535     gcry_mpi_t a2 = mpi_alloc_set_ui( 2 );
536     gcry_mpi_t q;
537     unsigned i, j, k;
538     int rc = 0;
539     unsigned nbits = mpi_get_nbits( n );
540
541     mpi_sub_ui( nminus1, n, 1 );
542
543     /* find q and k, so that n = 1 + 2^k * q */
544     q = mpi_copy( nminus1 );
545     k = mpi_trailing_zeros( q );
546     mpi_tdiv_q_2exp(q, q, k);
547
548     for(i=0 ; i < steps; i++ ) {
549         ++*count;
550         if( !i ) {
551             mpi_set_ui( x, 2 );
552         }
553         else {
554             gcry_mpi_randomize( x, nbits, GCRY_WEAK_RANDOM );
555
556             /* make sure that the number is smaller than the prime
557              * and keep the randomness of the high bit */
558             if( mpi_test_bit( x, nbits-2 ) ) {
559                 mpi_set_highbit( x, nbits-2 ); /* clear all higher bits */
560             }
561             else {
562                 mpi_set_highbit( x, nbits-2 );
563                 mpi_clear_bit( x, nbits-2 );
564             }
565             assert( mpi_cmp( x, nminus1 ) < 0 && mpi_cmp_ui( x, 1 ) > 0 );
566         }
567         gcry_mpi_powm( y, x, q, n);
568         if( mpi_cmp_ui(y, 1) && mpi_cmp( y, nminus1 ) ) {
569             for( j=1; j < k && mpi_cmp( y, nminus1 ); j++ ) {
570                 gcry_mpi_powm(y, y, a2, n);
571                 if( !mpi_cmp_ui( y, 1 ) )
572                     goto leave; /* not a prime */
573             }
574             if( mpi_cmp( y, nminus1 ) )
575                 goto leave; /* not a prime */
576         }
577         progress('+');
578     }
579     rc = 1; /* may be a prime */
580
581   leave:
582     mpi_free( x );
583     mpi_free( y );
584     mpi_free( z );
585     mpi_free( nminus1 );
586     mpi_free( q );
587
588     return rc;
589 }
590
591
592 static void
593 m_out_of_n( char *array, int m, int n )
594 {
595     int i=0, i1=0, j=0, jp=0,  j1=0, k1=0, k2=0;
596
597     if( !m || m >= n )
598         return;
599
600     if( m == 1 ) { /* special case */
601         for(i=0; i < n; i++ )
602             if( array[i] ) {
603                 array[i++] = 0;
604                 if( i >= n )
605                     i = 0;
606                 array[i] = 1;
607                 return;
608             }
609         BUG();
610     }
611
612     for(j=1; j < n; j++ ) {
613         if( array[n-1] == array[n-j-1] )
614             continue;
615         j1 = j;
616         break;
617     }
618
619     if( m & 1 ) { /* m is odd */
620         if( array[n-1] ) {
621             if( j1 & 1 ) {
622                 k1 = n - j1;
623                 k2 = k1+2;
624                 if( k2 > n )
625                     k2 = n;
626                 goto leave;
627             }
628             goto scan;
629         }
630         k2 = n - j1 - 1;
631         if( k2 == 0 ) {
632             k1 = i;
633             k2 = n - j1;
634         }
635         else if( array[k2] && array[k2-1] )
636             k1 = n;
637         else
638             k1 = k2 + 1;
639     }
640     else { /* m is even */
641         if( !array[n-1] ) {
642             k1 = n - j1;
643             k2 = k1 + 1;
644             goto leave;
645         }
646
647         if( !(j1 & 1) ) {
648             k1 = n - j1;
649             k2 = k1+2;
650             if( k2 > n )
651                 k2 = n;
652             goto leave;
653         }
654       scan:
655         jp = n - j1 - 1;
656         for(i=1; i <= jp; i++ ) {
657             i1 = jp + 2 - i;
658             if( array[i1-1]  ) {
659                 if( array[i1-2] ) {
660                     k1 = i1 - 1;
661                     k2 = n - j1;
662                 }
663                 else {
664                     k1 = i1 - 1;
665                     k2 = n + 1 - j1;
666                 }
667                 goto leave;
668             }
669         }
670         k1 = 1;
671         k2 = n + 1 - m;
672     }
673   leave:
674     array[k1-1] = !array[k1-1];
675     array[k2-1] = !array[k2-1];
676 }
677