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