* rndegd.c (rndegd_constructor): Fixed name of register function
[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 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 static void (*progress_cb) ( void *, 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
129
130 void
131 _gcry_register_primegen_progress ( void (*cb)( void *, 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, c );
143     else
144         fputc( c, stderr );
145 }
146
147
148 /****************
149  * Generate a prime number (stored in secure memory)
150  */
151 MPI
152 _gcry_generate_secret_prime( unsigned  nbits )
153 {
154     MPI prime;
155
156     prime = gen_prime( nbits, 1, 2 );
157     progress('\n');
158     return prime;
159 }
160
161 MPI
162 _gcry_generate_public_prime( unsigned  nbits )
163 {
164     MPI prime;
165
166     prime = gen_prime( nbits, 0, 2 );
167     progress('\n');
168     return prime;
169 }
170
171
172 /****************
173  * We do not need to use the strongest RNG because we gain no extra
174  * security from it - The prime number is public and we could also
175  * offer the factors for those who are willing to check that it is
176  * indeed a strong prime.
177  *
178  * mode 0: Standard
179  *      1: Make sure that at least one factor is of size qbits.
180  */
181 MPI
182 _gcry_generate_elg_prime( int mode, unsigned pbits, unsigned qbits,
183                     MPI g, MPI **ret_factors )
184 {
185     int n;  /* number of factors */
186     int m;  /* number of primes in pool */
187     unsigned fbits; /* length of prime factors */
188     MPI *factors; /* current factors */
189     MPI *pool;  /* pool of primes */
190     MPI q;      /* first prime factor (variable)*/
191     MPI prime;  /* prime test value */
192     MPI q_factor; /* used for mode 1 */
193     byte *perms = NULL;
194     int i, j;
195     int count1, count2;
196     unsigned nprime;
197     unsigned req_qbits = qbits; /* the requested q bits size */
198     MPI val_2  = mpi_alloc_set_ui( 2 );
199
200     /* find number of needed prime factors */
201     for(n=1; (pbits - qbits - 1) / n  >= qbits; n++ )
202         ;
203     n--;
204     if( !n || (mode==1 && n < 2) )
205         log_fatal("can't gen prime with pbits=%u qbits=%u\n", pbits, qbits );
206     if( mode == 1 ) {
207         n--;
208         fbits = (pbits - 2*req_qbits -1) / n;
209         qbits =  pbits - req_qbits - n*fbits;
210     }
211     else {
212         fbits = (pbits - req_qbits -1) / n;
213         qbits = pbits - n*fbits;
214     }
215     if( DBG_CIPHER )
216         log_debug("gen prime: pbits=%u qbits=%u fbits=%u/%u n=%d\n",
217                     pbits, req_qbits, qbits, fbits, n  );
218     prime = gcry_mpi_new ( pbits  );
219     q = gen_prime( qbits, 0, 0 );
220     q_factor = mode==1? gen_prime( req_qbits, 0, 0 ) : NULL;
221
222     /* allocate an array to hold the factors + 2 for later usage */
223     factors = gcry_xcalloc( n+2, sizeof *factors );
224
225     /* make a pool of 3n+5 primes (this is an arbitrary value) */
226     m = n*3+5;
227     if( mode == 1 )
228         m += 5; /* need some more for DSA */
229     if( m < 25 )
230         m = 25;
231     pool = gcry_xcalloc( m , sizeof *pool );
232
233     /* permutate over the pool of primes */
234     count1=count2=0;
235     do {
236       next_try:
237         if( !perms ) {
238             /* allocate new primes */
239             for(i=0; i < m; i++ ) {
240                 mpi_free(pool[i]);
241                 pool[i] = NULL;
242             }
243             /* init m_out_of_n() */
244             perms = gcry_xcalloc( 1, m );
245             for(i=0; i < n; i++ ) {
246                 perms[i] = 1;
247                 pool[i] = gen_prime( fbits, 0, 1 );
248                 factors[i] = pool[i];
249             }
250         }
251         else {
252             m_out_of_n( perms, n, m );
253             for(i=j=0; i < m && j < n ; i++ )
254                 if( perms[i] ) {
255                     if( !pool[i] )
256                         pool[i] = gen_prime( fbits, 0, 1 );
257                     factors[j++] = pool[i];
258                 }
259             if( i == n ) {
260                 gcry_free(perms); perms = NULL;
261                 progress('!');
262                 goto next_try;  /* allocate new primes */
263             }
264         }
265
266         mpi_set( prime, q );
267         mpi_mul_ui( prime, prime, 2 );
268         if( mode == 1 )
269             mpi_mul( prime, prime, q_factor );
270         for(i=0; i < n; i++ )
271             mpi_mul( prime, prime, factors[i] );
272         mpi_add_ui( prime, prime, 1 );
273         nprime = mpi_get_nbits(prime);
274         if( nprime < pbits ) {
275             if( ++count1 > 20 ) {
276                 count1 = 0;
277                 qbits++;
278                 progress('>');
279                 mpi_free (q);
280                 q = gen_prime( qbits, 0, 0 );
281                 goto next_try;
282             }
283         }
284         else
285             count1 = 0;
286         if( nprime > pbits ) {
287             if( ++count2 > 20 ) {
288                 count2 = 0;
289                 qbits--;
290                 progress('<');
291                 mpi_free (q);
292                 q = gen_prime( qbits, 0, 0 );
293                 goto next_try;
294             }
295         }
296         else
297             count2 = 0;
298     } while( !(nprime == pbits && check_prime( prime, val_2 )) );
299
300     if( DBG_CIPHER ) {
301         progress('\n');
302         log_mpidump( "prime    : ", prime );
303         log_mpidump( "factor  q: ", q );
304         if( mode == 1 )
305             log_mpidump( "factor q0: ", q_factor );
306         for(i=0; i < n; i++ )
307             log_mpidump( "factor pi: ", factors[i] );
308         log_debug("bit sizes: prime=%u, q=%u", mpi_get_nbits(prime), mpi_get_nbits(q) );
309         if( mode == 1 )
310             fprintf(stderr, ", q0=%u", mpi_get_nbits(q_factor) );
311         for(i=0; i < n; i++ )
312             fprintf(stderr, ", p%d=%u", i, mpi_get_nbits(factors[i]) );
313         progress('\n');
314     }
315
316     if( ret_factors ) { /* caller wants the factors */
317         *ret_factors = gcry_xcalloc( n+2 , sizeof **ret_factors);
318         i = 0;
319         if( mode == 1 ) {
320             (*ret_factors)[i++] = mpi_copy( q_factor );
321             for(; i <= n; i++ )
322                 (*ret_factors)[i] = mpi_copy( factors[i] );
323         }
324         else {
325             for(; i < n; i++ )
326                 (*ret_factors)[i] = mpi_copy( factors[i] );
327         }
328     }
329
330     if( g ) { /* create a generator (start with 3)*/
331         MPI tmp   = mpi_alloc( mpi_get_nlimbs(prime) );
332         MPI b     = mpi_alloc( mpi_get_nlimbs(prime) );
333         MPI pmin1 = mpi_alloc( mpi_get_nlimbs(prime) );
334
335         if( mode == 1 )
336             BUG(); /* not yet implemented */
337         factors[n] = q;
338         factors[n+1] = mpi_alloc_set_ui(2);
339         mpi_sub_ui( pmin1, prime, 1 );
340         mpi_set_ui(g,2);
341         do {
342             mpi_add_ui(g, g, 1);
343             if( DBG_CIPHER ) {
344                 log_debug("checking g: ");
345                 /*mpi_print( stderr, g, 1 );*/
346 #if __GNUC__ >= 2
347 #   warning we need an internal mpi_print for debugging
348 #endif
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 MPI
384 gen_prime( unsigned  nbits, int secret, int randomlevel )
385 {
386     MPI prime, ptest, pminus1, val_2, val_3, result;
387     int i;
388     unsigned x, step;
389     unsigned count1, count2;
390     int *mods;
391
392     if( 0 && DBG_CIPHER )
393         log_debug("generate a prime of %u bits ", nbits );
394
395     if( !no_of_small_prime_numbers ) {
396         for(i=0; small_prime_numbers[i]; i++ )
397             no_of_small_prime_numbers++;
398     }
399     mods = gcry_xmalloc( no_of_small_prime_numbers * sizeof *mods );
400     /* make nbits fit into MPI implementation */
401     val_2  = mpi_alloc_set_ui( 2 );
402     val_3 = mpi_alloc_set_ui( 3);
403     prime  = secret? gcry_mpi_snew ( nbits ): gcry_mpi_new ( nbits );
404     result = mpi_alloc_like( prime );
405     pminus1= mpi_alloc_like( prime );
406     ptest  = mpi_alloc_like( prime );
407     count1 = count2 = 0;
408     for(;;) {  /* 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             /* check against all the small primes we have in mods */
430             count1++;
431             for(i=0; (x = small_prime_numbers[i]); i++ ) {
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 ) ) { /* not composite */
447                 /* perform stronger tests */
448                 if( is_prime(ptest, 5, &count2 ) ) {
449                     if( !mpi_test_bit( ptest, nbits-1-secret ) ) {
450                         progress('\n');
451                         log_debug("overflow in prime generation\n");
452                         break; /* step loop, continue with a new prime */
453                     }
454
455                     mpi_free(val_2);
456                     mpi_free(val_3);
457                     mpi_free(result);
458                     mpi_free(pminus1);
459                     mpi_free(prime);
460                     gcry_free(mods);
461                     return ptest;
462                 }
463             }
464             if( ++dotcount == 10 ) {
465                 progress('.');
466                 dotcount = 0;
467             }
468         }
469         progress(':'); /* restart with a new random value */
470     }
471 }
472
473 /****************
474  * Returns: true if this may be a prime
475  */
476 static int
477 check_prime( MPI prime, MPI val_2 )
478 {
479     int i;
480     unsigned x;
481     int count=0;
482
483     /* check against small primes */
484     for(i=0; (x = small_prime_numbers[i]); i++ ) {
485         if( mpi_divisible_ui( prime, x ) )
486             return 0;
487     }
488
489     /* a quick fermat test */
490     {
491         MPI result = mpi_alloc_like( prime );
492         MPI pminus1 = mpi_alloc_like( prime );
493         mpi_sub_ui( pminus1, prime, 1);
494         gcry_mpi_powm( result, val_2, pminus1, prime );
495         mpi_free( pminus1 );
496         if( mpi_cmp_ui( result, 1 ) ) { /* if composite */
497             mpi_free( result );
498             progress('.');
499             return 0;
500         }
501         mpi_free( result );
502     }
503
504     /* perform stronger tests */
505     if( is_prime(prime, 5, &count ) )
506         return 1; /* is probably a prime */
507     progress('.');
508     return 0;
509 }
510
511
512 /****************
513  * Return true if n is probably a prime
514  */
515 static int
516 is_prime( MPI n, int steps, int *count )
517 {
518     MPI x = mpi_alloc( mpi_get_nlimbs( n ) );
519     MPI y = mpi_alloc( mpi_get_nlimbs( n ) );
520     MPI z = mpi_alloc( mpi_get_nlimbs( n ) );
521     MPI nminus1 = mpi_alloc( mpi_get_nlimbs( n ) );
522     MPI a2 = mpi_alloc_set_ui( 2 );
523     MPI q;
524     unsigned i, j, k;
525     int rc = 0;
526     unsigned nbits = mpi_get_nbits( n );
527
528     mpi_sub_ui( nminus1, n, 1 );
529
530     /* find q and k, so that n = 1 + 2^k * q */
531     q = mpi_copy( nminus1 );
532     k = mpi_trailing_zeros( q );
533     mpi_tdiv_q_2exp(q, q, k);
534
535     for(i=0 ; i < steps; i++ ) {
536         ++*count;
537         if( !i ) {
538             mpi_set_ui( x, 2 );
539         }
540         else {
541             gcry_mpi_randomize( x, nbits, GCRY_WEAK_RANDOM );
542
543             /* make sure that the number is smaller than the prime
544              * and keep the randomness of the high bit */
545             if( mpi_test_bit( x, nbits-2 ) ) {
546                 mpi_set_highbit( x, nbits-2 ); /* clear all higher bits */
547             }
548             else {
549                 mpi_set_highbit( x, nbits-2 );
550                 mpi_clear_bit( x, nbits-2 );
551             }
552             assert( mpi_cmp( x, nminus1 ) < 0 && mpi_cmp_ui( x, 1 ) > 0 );
553         }
554         gcry_mpi_powm( y, x, q, n);
555         if( mpi_cmp_ui(y, 1) && mpi_cmp( y, nminus1 ) ) {
556             for( j=1; j < k && mpi_cmp( y, nminus1 ); j++ ) {
557                 gcry_mpi_powm(y, y, a2, n);
558                 if( !mpi_cmp_ui( y, 1 ) )
559                     goto leave; /* not a prime */
560             }
561             if( mpi_cmp( y, nminus1 ) )
562                 goto leave; /* not a prime */
563         }
564         progress('+');
565     }
566     rc = 1; /* may be a prime */
567
568   leave:
569     mpi_free( x );
570     mpi_free( y );
571     mpi_free( z );
572     mpi_free( nminus1 );
573     mpi_free( q );
574
575     return rc;
576 }
577
578
579 static void
580 m_out_of_n( char *array, int m, int n )
581 {
582     int i=0, i1=0, j=0, jp=0,  j1=0, k1=0, k2=0;
583
584     if( !m || m >= n )
585         return;
586
587     if( m == 1 ) { /* special case */
588         for(i=0; i < n; i++ )
589             if( array[i] ) {
590                 array[i++] = 0;
591                 if( i >= n )
592                     i = 0;
593                 array[i] = 1;
594                 return;
595             }
596         BUG();
597     }
598
599     for(j=1; j < n; j++ ) {
600         if( array[n-1] == array[n-j-1] )
601             continue;
602         j1 = j;
603         break;
604     }
605
606     if( m & 1 ) { /* m is odd */
607         if( array[n-1] ) {
608             if( j1 & 1 ) {
609                 k1 = n - j1;
610                 k2 = k1+2;
611                 if( k2 > n )
612                     k2 = n;
613                 goto leave;
614             }
615             goto scan;
616         }
617         k2 = n - j1 - 1;
618         if( k2 == 0 ) {
619             k1 = i;
620             k2 = n - j1;
621         }
622         else if( array[k2] && array[k2-1] )
623             k1 = n;
624         else
625             k1 = k2 + 1;
626     }
627     else { /* m is even */
628         if( !array[n-1] ) {
629             k1 = n - j1;
630             k2 = k1 + 1;
631             goto leave;
632         }
633
634         if( !(j1 & 1) ) {
635             k1 = n - j1;
636             k2 = k1+2;
637             if( k2 > n )
638                 k2 = n;
639             goto leave;
640         }
641       scan:
642         jp = n - j1 - 1;
643         for(i=1; i <= jp; i++ ) {
644             i1 = jp + 2 - i;
645             if( array[i1-1]  ) {
646                 if( array[i1-2] ) {
647                     k1 = i1 - 1;
648                     k2 = n - j1;
649                 }
650                 else {
651                     k1 = i1 - 1;
652                     k2 = n + 1 - j1;
653                 }
654                 goto leave;
655             }
656         }
657         k1 = 1;
658         k2 = n + 1 - m;
659     }
660   leave:
661     array[k1-1] = !array[k1-1];
662     array[k2-1] = !array[k2-1];
663 }
664