* gcrypt.h (gcry_mpi_scan, gcry_mpi_print): API change.
[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                 gcry_mpi_dump (g);
348                 log_debug ("\n");
349             }
350             else
351                 progress('^');
352             for(i=0; i < n+2; i++ ) {
353                 /*fputc('~', stderr);*/
354                 mpi_fdiv_q(tmp, pmin1, factors[i] );
355                 /* (no mpi_pow(), but it is okay to use this with mod prime) */
356                 gcry_mpi_powm(b, g, tmp, prime );
357                 if( !mpi_cmp_ui(b, 1) )
358                     break;
359             }
360             if( DBG_CIPHER )
361                 progress('\n');
362         } while( i < n+2 );
363         mpi_free(factors[n+1]);
364         mpi_free(tmp);
365         mpi_free(b);
366         mpi_free(pmin1);
367     }
368     if( !DBG_CIPHER )
369         progress('\n');
370
371     gcry_free( factors );  /* (factors are shallow copies) */
372     for(i=0; i < m; i++ )
373         mpi_free( pool[i] );
374     gcry_free( pool );
375     gcry_free(perms);
376     mpi_free(val_2);
377     mpi_free (q);
378     return prime;
379 }
380
381
382
383 static gcry_mpi_t
384 gen_prime (unsigned int nbits, int secret, int randomlevel, 
385            int (*extra_check)(void *, gcry_mpi_t), void *extra_check_arg)
386 {
387   gcry_mpi_t prime, ptest, pminus1, val_2, val_3, result;
388   int i;
389   unsigned x, step;
390   unsigned count1, count2;
391   int *mods;
392   
393   if( 0 && DBG_CIPHER )
394     log_debug("generate a prime of %u bits ", nbits );
395
396   mods = gcry_xmalloc( no_of_small_prime_numbers * sizeof *mods );
397   /* make nbits fit into gcry_mpi_t implementation */
398   val_2  = mpi_alloc_set_ui( 2 );
399   val_3 = mpi_alloc_set_ui( 3);
400   prime  = secret? gcry_mpi_snew ( nbits ): gcry_mpi_new ( nbits );
401   result = mpi_alloc_like( prime );
402   pminus1= mpi_alloc_like( prime );
403   ptest  = mpi_alloc_like( prime );
404   count1 = count2 = 0;
405   for (;;)
406     {  /* try forvever */
407       int dotcount=0;
408       
409       /* generate a random number */
410       gcry_mpi_randomize( prime, nbits, randomlevel );
411       
412       /* Set high order bit to 1, set low order bit to 0.  If we are
413          generating a secret prime we are most probably doing that
414          for RSA, to make sure that the modulus does have the
415          requested keysize we set the 2 high order bits */
416       mpi_set_highbit (prime, nbits-1);
417       if (secret)
418         mpi_set_bit (prime, nbits-2);
419       mpi_set_bit(prime, 0);
420       
421       /* calculate all remainders */
422       for (i=0; (x = small_prime_numbers[i]); i++ )
423         mods[i] = mpi_fdiv_r_ui(NULL, prime, x);
424       
425       /* now try some primes starting with prime */
426       for(step=0; step < 20000; step += 2 ) 
427         {
428           /* check against all the small primes we have in mods */
429           count1++;
430           for (i=0; (x = small_prime_numbers[i]); i++ ) 
431             {
432               while ( mods[i] + step >= x )
433                 mods[i] -= x;
434               if ( !(mods[i] + step) )
435                 break;
436             }
437           if ( x )
438             continue;   /* found a multiple of an already known prime */
439
440           mpi_add_ui( ptest, prime, step );
441
442           /* do a faster Fermat test */
443           count2++;
444           mpi_sub_ui( pminus1, ptest, 1);
445           gcry_mpi_powm( result, val_2, pminus1, ptest );
446           if ( !mpi_cmp_ui( result, 1 ) )
447             { /* not composite, perform stronger tests */
448                 if (is_prime(ptest, 5, &count2 ))
449                   {
450                     if (!mpi_test_bit( ptest, nbits-1-secret ))
451                       {
452                         progress('\n');
453                         log_debug("overflow in prime generation\n");
454                         break; /* stop loop, continue with a new prime */
455                       }
456
457                     if (extra_check && extra_check (extra_check_arg, ptest))
458                       { /* The extra check told us that this prime is
459                            not of the caller's taste. */
460                         progress ('/');
461                       }
462                     else
463                       { /* got it */
464                         mpi_free(val_2);
465                         mpi_free(val_3);
466                         mpi_free(result);
467                         mpi_free(pminus1);
468                         mpi_free(prime);
469                         gcry_free(mods);
470                         return ptest; 
471                       }
472                   }
473             }
474           if (++dotcount == 10 )
475             {
476               progress('.');
477               dotcount = 0;
478             }
479         }
480       progress(':'); /* restart with a new random value */
481     }
482 }
483
484 /****************
485  * Returns: true if this may be a prime
486  */
487 static int
488 check_prime( gcry_mpi_t prime, gcry_mpi_t val_2 )
489 {
490     int i;
491     unsigned x;
492     int count=0;
493
494     /* check against small primes */
495     for(i=0; (x = small_prime_numbers[i]); i++ ) {
496         if( mpi_divisible_ui( prime, x ) )
497             return 0;
498     }
499
500     /* a quick fermat test */
501     {
502         gcry_mpi_t result = mpi_alloc_like( prime );
503         gcry_mpi_t pminus1 = mpi_alloc_like( prime );
504         mpi_sub_ui( pminus1, prime, 1);
505         gcry_mpi_powm( result, val_2, pminus1, prime );
506         mpi_free( pminus1 );
507         if( mpi_cmp_ui( result, 1 ) ) { /* if composite */
508             mpi_free( result );
509             progress('.');
510             return 0;
511         }
512         mpi_free( result );
513     }
514
515     /* perform stronger tests */
516     if( is_prime(prime, 5, &count ) )
517         return 1; /* is probably a prime */
518     progress('.');
519     return 0;
520 }
521
522
523 /****************
524  * Return true if n is probably a prime
525  */
526 static int
527 is_prime( gcry_mpi_t n, int steps, int *count )
528 {
529     gcry_mpi_t x = mpi_alloc( mpi_get_nlimbs( n ) );
530     gcry_mpi_t y = mpi_alloc( mpi_get_nlimbs( n ) );
531     gcry_mpi_t z = mpi_alloc( mpi_get_nlimbs( n ) );
532     gcry_mpi_t nminus1 = mpi_alloc( mpi_get_nlimbs( n ) );
533     gcry_mpi_t a2 = mpi_alloc_set_ui( 2 );
534     gcry_mpi_t q;
535     unsigned i, j, k;
536     int rc = 0;
537     unsigned nbits = mpi_get_nbits( n );
538
539     mpi_sub_ui( nminus1, n, 1 );
540
541     /* find q and k, so that n = 1 + 2^k * q */
542     q = mpi_copy( nminus1 );
543     k = mpi_trailing_zeros( q );
544     mpi_tdiv_q_2exp(q, q, k);
545
546     for(i=0 ; i < steps; i++ ) {
547         ++*count;
548         if( !i ) {
549             mpi_set_ui( x, 2 );
550         }
551         else {
552             gcry_mpi_randomize( x, nbits, GCRY_WEAK_RANDOM );
553
554             /* make sure that the number is smaller than the prime
555              * and keep the randomness of the high bit */
556             if( mpi_test_bit( x, nbits-2 ) ) {
557                 mpi_set_highbit( x, nbits-2 ); /* clear all higher bits */
558             }
559             else {
560                 mpi_set_highbit( x, nbits-2 );
561                 mpi_clear_bit( x, nbits-2 );
562             }
563             assert( mpi_cmp( x, nminus1 ) < 0 && mpi_cmp_ui( x, 1 ) > 0 );
564         }
565         gcry_mpi_powm( y, x, q, n);
566         if( mpi_cmp_ui(y, 1) && mpi_cmp( y, nminus1 ) ) {
567             for( j=1; j < k && mpi_cmp( y, nminus1 ); j++ ) {
568                 gcry_mpi_powm(y, y, a2, n);
569                 if( !mpi_cmp_ui( y, 1 ) )
570                     goto leave; /* not a prime */
571             }
572             if( mpi_cmp( y, nminus1 ) )
573                 goto leave; /* not a prime */
574         }
575         progress('+');
576     }
577     rc = 1; /* may be a prime */
578
579   leave:
580     mpi_free( x );
581     mpi_free( y );
582     mpi_free( z );
583     mpi_free( nminus1 );
584     mpi_free( q );
585
586     return rc;
587 }
588
589
590 static void
591 m_out_of_n( char *array, int m, int n )
592 {
593     int i=0, i1=0, j=0, jp=0,  j1=0, k1=0, k2=0;
594
595     if( !m || m >= n )
596         return;
597
598     if( m == 1 ) { /* special case */
599         for(i=0; i < n; i++ )
600             if( array[i] ) {
601                 array[i++] = 0;
602                 if( i >= n )
603                     i = 0;
604                 array[i] = 1;
605                 return;
606             }
607         BUG();
608     }
609
610     for(j=1; j < n; j++ ) {
611         if( array[n-1] == array[n-j-1] )
612             continue;
613         j1 = j;
614         break;
615     }
616
617     if( m & 1 ) { /* m is odd */
618         if( array[n-1] ) {
619             if( j1 & 1 ) {
620                 k1 = n - j1;
621                 k2 = k1+2;
622                 if( k2 > n )
623                     k2 = n;
624                 goto leave;
625             }
626             goto scan;
627         }
628         k2 = n - j1 - 1;
629         if( k2 == 0 ) {
630             k1 = i;
631             k2 = n - j1;
632         }
633         else if( array[k2] && array[k2-1] )
634             k1 = n;
635         else
636             k1 = k2 + 1;
637     }
638     else { /* m is even */
639         if( !array[n-1] ) {
640             k1 = n - j1;
641             k2 = k1 + 1;
642             goto leave;
643         }
644
645         if( !(j1 & 1) ) {
646             k1 = n - j1;
647             k2 = k1+2;
648             if( k2 > n )
649                 k2 = n;
650             goto leave;
651         }
652       scan:
653         jp = n - j1 - 1;
654         for(i=1; i <= jp; i++ ) {
655             i1 = jp + 2 - i;
656             if( array[i1-1]  ) {
657                 if( array[i1-2] ) {
658                     k1 = i1 - 1;
659                     k2 = n - j1;
660                 }
661                 else {
662                     k1 = i1 - 1;
663                     k2 = n + 1 - j1;
664                 }
665                 goto leave;
666             }
667         }
668         k1 = 1;
669         k2 = n + 1 - m;
670     }
671   leave:
672     array[k1-1] = !array[k1-1];
673     array[k2-1] = !array[k2-1];
674 }
675