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