Cleaned up the public key module calling conventions.
[libgcrypt.git] / cipher / primegen.c
1 /* primegen.c - prime number generator
2  * Copyright (C) 1998, 2000, 2001, 2002, 2003
3  *               2004, 2008 Free Software Foundation, Inc.
4  *
5  * This file is part of Libgcrypt.
6  *
7  * Libgcrypt is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU Lesser general Public License as
9  * published by the Free Software Foundation; either version 2.1 of
10  * the License, or (at your option) any later version.
11  *
12  * Libgcrypt is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  * GNU Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
20  */
21
22 #include <config.h>
23
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include <errno.h>
28
29 #include "g10lib.h"
30 #include "mpi.h"
31 #include "cipher.h"
32 #include "ath.h"
33
34 static gcry_mpi_t gen_prime (unsigned int nbits, int secret, int randomlevel, 
35                              int (*extra_check)(void *, gcry_mpi_t),
36                              void *extra_check_arg);
37 static int check_prime( gcry_mpi_t prime, gcry_mpi_t val_2, int rm_rounds,
38                         gcry_prime_check_func_t cb_func, void *cb_arg );
39 static int is_prime (gcry_mpi_t n, int steps, unsigned int *count);
40 static void m_out_of_n( char *array, int m, int n );
41
42 static void (*progress_cb) (void *,const char*,int,int, int );
43 static void *progress_cb_data;
44
45 /* Note: 2 is not included because it can be tested more easily by
46    looking at bit 0. The last entry in this list is marked by a zero */
47 static ushort small_prime_numbers[] = {
48     3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
49     47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101,
50     103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
51     157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
52     211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
53     269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
54     331, 337, 347, 349, 353, 359, 367, 373, 379, 383,
55     389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
56     449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
57     509, 521, 523, 541, 547, 557, 563, 569, 571, 577,
58     587, 593, 599, 601, 607, 613, 617, 619, 631, 641,
59     643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
60     709, 719, 727, 733, 739, 743, 751, 757, 761, 769,
61     773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
62     853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
63     919, 929, 937, 941, 947, 953, 967, 971, 977, 983,
64     991, 997, 1009, 1013, 1019, 1021, 1031, 1033,
65     1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091,
66     1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
67     1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213,
68     1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,
69     1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307,
70     1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
71     1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
72     1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
73     1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
74     1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609,
75     1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667,
76     1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
77     1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
78     1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871,
79     1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931,
80     1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
81     1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
82     2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111,
83     2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
84     2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243,
85     2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
86     2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
87     2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411,
88     2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473,
89     2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551,
90     2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633,
91     2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
92     2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729,
93     2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791,
94     2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
95     2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917,
96     2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
97     3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061,
98     3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137,
99     3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
100     3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271,
101     3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
102     3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391,
103     3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467,
104     3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533,
105     3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,
106     3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
107     3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,
108     3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779,
109     3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851,
110     3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917,
111     3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
112     4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049,
113     4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111,
114     4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177,
115     4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243,
116     4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
117     4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391,
118     4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457,
119     4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,
120     4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597,
121     4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
122     4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729,
123     4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799,
124     4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889,
125     4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951,
126     4957, 4967, 4969, 4973, 4987, 4993, 4999,
127     0
128 };
129 static int no_of_small_prime_numbers = DIM (small_prime_numbers) - 1;
130
131
132 \f
133 /* An object and a list to build up a global pool of primes.  See
134    save_pool_prime and get_pool_prime. */
135 struct primepool_s 
136 {
137   struct primepool_s *next;
138   gcry_mpi_t prime;      /* If this is NULL the entry is not used. */
139   unsigned int nbits;
140   gcry_random_level_t randomlevel;
141 };
142 struct primepool_s *primepool;
143 /* Mutex used to protect access to the primepool.  */
144 static ath_mutex_t primepool_lock = ATH_MUTEX_INITIALIZER;
145
146
147
148 /* Save PRIME which has been generated at RANDOMLEVEL for later
149    use. Needs to be called while primepool_lock is being hold.  Note
150    that PRIME should be considered released after calling this
151    function. */
152 static void
153 save_pool_prime (gcry_mpi_t prime, gcry_random_level_t randomlevel)
154 {
155   struct primepool_s *item, *item2;
156   size_t n;
157
158   for (n=0, item = primepool; item; item = item->next, n++)
159     if (!item->prime)
160       break;
161   if (!item && n > 100)
162     {
163       /* Remove some of the entries.  Our strategy is removing
164          the last third from the list. */
165       int i;
166       
167       for (i=0, item2 = primepool; item2; item2 = item2->next)
168         {
169           if (i >= n/3*2)
170             {
171               gcry_mpi_release (item2->prime);
172               item2->prime = NULL;
173               if (!item)
174                 item = item2;
175             }
176         }
177     }
178   if (!item)
179     {
180       item = gcry_calloc (1, sizeof *item);
181       if (!item)
182         {
183           /* Out of memory.  Silently giving up. */
184           gcry_mpi_release (prime);
185           return; 
186         }
187       item->next = primepool;
188       primepool = item;
189     }
190   item->prime = prime;
191   item->nbits = mpi_get_nbits (prime);
192   item->randomlevel = randomlevel;
193 }
194
195
196 /* Return a prime for the prime pool or NULL if none has been found.
197    The prime needs to match NBITS and randomlevel. This function needs
198    to be called why the primepool_look is being hold. */
199 static gcry_mpi_t
200 get_pool_prime (unsigned int nbits, gcry_random_level_t randomlevel)
201 {
202   struct primepool_s *item;
203
204   for (item = primepool; item; item = item->next)
205     if (item->prime
206         && item->nbits == nbits && item->randomlevel == randomlevel)
207       {
208         gcry_mpi_t prime = item->prime;
209         item->prime = NULL;
210         gcry_assert (nbits == mpi_get_nbits (prime));
211         return prime;
212       }
213   return NULL;
214 }
215
216
217
218
219
220 \f
221 void
222 _gcry_register_primegen_progress ( void (*cb)(void *,const char*,int,int,int),
223                                    void *cb_data )
224 {
225   progress_cb = cb;
226   progress_cb_data = cb_data;
227 }
228
229
230 static void
231 progress( int c )
232 {
233   if ( progress_cb )
234     progress_cb ( progress_cb_data, "primegen", c, 0, 0 );
235 }
236
237
238 /****************
239  * Generate a prime number (stored in secure memory)
240  */
241 gcry_mpi_t
242 _gcry_generate_secret_prime (unsigned int nbits,
243                              gcry_random_level_t random_level,
244                              int (*extra_check)(void*, gcry_mpi_t),
245                              void *extra_check_arg)
246 {
247   gcry_mpi_t prime;
248
249   prime = gen_prime (nbits, 1, random_level, extra_check, extra_check_arg);
250   progress('\n');
251   return prime;
252 }
253
254
255 /* Generate a prime number which may be public, i.e. not allocated in
256    secure memory.  */
257 gcry_mpi_t
258 _gcry_generate_public_prime (unsigned int nbits,
259                              gcry_random_level_t random_level,
260                              int (*extra_check)(void*, gcry_mpi_t),
261                              void *extra_check_arg)
262 {
263   gcry_mpi_t prime;
264
265   prime = gen_prime (nbits, 0, random_level, extra_check, extra_check_arg);
266   progress('\n');
267   return prime;
268 }
269
270
271 /* Core prime generation function.  The algorithm used to generate
272    practically save primes is due to Lim and Lee as described in the
273    CRYPTO '97 proceedings (ISBN3540633847) page 260.
274
275    NEED_Q_FACTOR: If true make sure that at least one factor is of
276                   size qbits.  This is for example required for DSA.
277    PRIME_GENERATED: Adresss of a variable where the resulting prime
278                     number will be stored.
279    PBITS: Requested size of the prime number.  At least 48.
280    QBITS: One factor of the prime needs to be of this size.  Maybe 0
281           if this is not required.  See also MODE.
282    G: If not NULL an MPI which will receive a generator for the prime
283       for use with Elgamal.
284    RET_FACTORS: if not NULL, an array with all factors are stored at
285                 that address.
286    ALL_FACTORS: If set to true all factors of prime-1 are returned.
287    RANDOMLEVEL:  How strong should the random numers be.
288    FLAGS: Prime generation bit flags. Currently supported:
289           GCRY_PRIME_FLAG_SECRET - The prime needs to be kept secret.
290    CB_FUNC, CB_ARG:  Callback to be used for extra checks.
291
292  */
293 static gcry_err_code_t
294 prime_generate_internal (int need_q_factor,
295                          gcry_mpi_t *prime_generated, unsigned int pbits,
296                          unsigned int qbits, gcry_mpi_t g,
297                          gcry_mpi_t **ret_factors,
298                          gcry_random_level_t randomlevel, unsigned int flags,
299                          int all_factors,
300                          gcry_prime_check_func_t cb_func, void *cb_arg)
301 {
302   gcry_err_code_t err = 0;
303   gcry_mpi_t *factors_new = NULL; /* Factors to return to the
304                                      caller.  */
305   gcry_mpi_t *factors = NULL;   /* Current factors.  */
306   gcry_random_level_t poolrandomlevel; /* Random level used for pool primes. */
307   gcry_mpi_t *pool = NULL;      /* Pool of primes.  */
308   int *pool_in_use = NULL;      /* Array with currently used POOL elements. */
309   unsigned char *perms = NULL;  /* Permutations of POOL.  */
310   gcry_mpi_t q_factor = NULL;   /* Used if QBITS is non-zero.  */
311   unsigned int fbits = 0;       /* Length of prime factors.  */
312   unsigned int n = 0;           /* Number of factors.  */
313   unsigned int m = 0;           /* Number of primes in pool.  */
314   gcry_mpi_t q = NULL;          /* First prime factor.  */
315   gcry_mpi_t prime = NULL;      /* Prime candidate.  */
316   unsigned int nprime = 0;      /* Bits of PRIME.  */
317   unsigned int req_qbits;       /* The original QBITS value.  */
318   gcry_mpi_t val_2;             /* For check_prime().  */
319   int is_locked = 0;            /* Flag to help unlocking the primepool. */
320   unsigned int is_secret = (flags & GCRY_PRIME_FLAG_SECRET);
321   unsigned int count1 = 0, count2 = 0;
322   unsigned int i = 0, j = 0;
323
324   if (pbits < 48)
325     return GPG_ERR_INV_ARG;
326
327   /* We won't use a too strong random elvel for the pooled subprimes. */
328   poolrandomlevel = (randomlevel > GCRY_STRONG_RANDOM?
329                      GCRY_STRONG_RANDOM : randomlevel);
330
331
332   /* If QBITS is not given, assume a reasonable value. */
333   if (!qbits)
334     qbits = pbits / 3;
335
336   req_qbits = qbits;
337
338   /* Find number of needed prime factors N.  */
339   for (n = 1; (pbits - qbits - 1) / n  >= qbits; n++)
340     ;
341   n--;
342
343   val_2 = mpi_alloc_set_ui (2);
344
345   if ((! n) || ((need_q_factor) && (n < 2)))
346     {
347       err = GPG_ERR_INV_ARG;
348       goto leave;
349     }
350
351   if (need_q_factor)
352     {
353       n--;  /* Need one factor less because we want a specific Q-FACTOR. */
354       fbits = (pbits - 2 * req_qbits -1) / n;
355       qbits =  pbits - req_qbits - n * fbits;
356     }
357   else
358     {
359       fbits = (pbits - req_qbits -1) / n;
360       qbits = pbits - n * fbits;
361     }
362   
363   if (DBG_CIPHER)
364     log_debug ("gen prime: pbits=%u qbits=%u fbits=%u/%u n=%d\n",
365                pbits, req_qbits, qbits, fbits, n);
366
367   /* Allocate an integer to old the new prime. */
368   prime = gcry_mpi_new (pbits);
369
370   /* Generate first prime factor.  */
371   q = gen_prime (qbits, is_secret, randomlevel, NULL, NULL);
372
373   /* Generate a specific Q-Factor if requested. */
374   if (need_q_factor)
375     q_factor = gen_prime (req_qbits, is_secret, randomlevel, NULL, NULL);
376   
377   /* Allocate an array to hold all factors + 2 for later usage.  */
378   factors = gcry_calloc (n + 2, sizeof (*factors));
379   if (!factors)
380     {
381       err = gpg_err_code_from_errno (errno);
382       goto leave;
383     }
384
385   /* Allocate an array to track pool usage. */
386   pool_in_use = gcry_malloc (n * sizeof *pool_in_use);
387   if (!pool_in_use)
388     {
389       err = gpg_err_code_from_errno (errno);
390       goto leave;
391     }
392   for (i=0; i < n; i++)
393     pool_in_use[i] = -1;
394       
395   /* Make a pool of 3n+5 primes (this is an arbitrary value).  We
396      require at least 30 primes for are useful selection process. 
397      
398      Fixme: We need to research the best formula for sizing the pool.
399   */
400   m = n * 3 + 5;
401   if (need_q_factor) /* Need some more in this case. */
402     m += 5;
403   if (m < 30)
404     m = 30;
405   pool = gcry_calloc (m , sizeof (*pool));
406   if (! pool)
407     {
408       err = gpg_err_code_from_errno (errno);
409       goto leave;
410     }
411
412   /* Permutate over the pool of primes until we find a prime of the
413      requested length.  */
414   do
415     {
416     next_try:
417       for (i=0; i < n; i++)
418         pool_in_use[i] = -1;
419
420       if (!perms)
421         {
422           /* Allocate new primes.  This is done right at the beginning
423              of the loop and if we have later run out of primes. */
424           for (i = 0; i < m; i++)
425             {
426               mpi_free (pool[i]);
427               pool[i] = NULL;
428             }
429
430           /* Init m_out_of_n().  */
431           perms = gcry_calloc (1, m);
432           if (!perms)
433             {
434               err = gpg_err_code_from_errno (errno);
435               goto leave;
436             }
437
438           if (ath_mutex_lock (&primepool_lock))
439             {
440               err = GPG_ERR_INTERNAL;
441               goto leave;
442             }
443           is_locked = 1;
444           for (i = 0; i < n; i++)
445             {
446               perms[i] = 1; 
447               /* At a maximum we use strong random for the factors.
448                  This saves us a lot of entropy. Given that Q and
449                  possible Q-factor are also used in the final prime
450                  this should be acceptable.  We also don't allocate in
451                  secure memory to save on that scare resource too.  If
452                  Q has been allocated in secure memory, the final
453                  prime will be saved there anyway.  This is because
454                  our MPI routines take care of that.  GnuPG has worked
455                  this way ever since.  */
456               pool[i] = NULL;
457               if (is_locked)
458                 {
459                   pool[i] = get_pool_prime (fbits, poolrandomlevel);
460                   if (!pool[i])
461                     {
462                       if (ath_mutex_unlock (&primepool_lock))
463                         {
464                           err = GPG_ERR_INTERNAL;
465                           goto leave;
466                         }
467                       is_locked = 0;
468                     }
469                 }
470               if (!pool[i])
471                 pool[i] = gen_prime (fbits, 0, poolrandomlevel, NULL, NULL);
472               pool_in_use[i] = i;
473               factors[i] = pool[i];
474             }
475           if (is_locked && ath_mutex_unlock (&primepool_lock))
476             {
477               err = GPG_ERR_INTERNAL;
478               goto leave;
479             }
480           is_locked = 0;
481         }
482       else
483         {
484           /* Get next permutation. */
485           m_out_of_n ( (char*)perms, n, m);
486           if (ath_mutex_lock (&primepool_lock))
487             {
488               err = GPG_ERR_INTERNAL;
489               goto leave;
490             }
491           is_locked = 1;
492           for (i = j = 0; (i < m) && (j < n); i++)
493             if (perms[i])
494               {
495                 /* If the subprime has not yet beed generated do it now. */
496                 if (!pool[i] && is_locked)
497                   {
498                     pool[i] = get_pool_prime (fbits, poolrandomlevel);
499                     if (!pool[i])
500                       {
501                         if (ath_mutex_unlock (&primepool_lock))
502                           {
503                             err = GPG_ERR_INTERNAL;
504                             goto leave;
505                           }
506                         is_locked = 0;
507                       }
508                   }
509                 if (!pool[i])
510                   pool[i] = gen_prime (fbits, 0, poolrandomlevel, NULL, NULL);
511                 pool_in_use[j] = i;
512                 factors[j++] = pool[i];
513               }
514           if (is_locked && ath_mutex_unlock (&primepool_lock))
515             {
516               err = GPG_ERR_INTERNAL;
517               goto leave;
518             }
519           is_locked = 0;
520           if (i == n)
521             {
522               /* Ran out of permutations: Allocate new primes.  */
523               gcry_free (perms);
524               perms = NULL;
525               progress ('!');
526               goto next_try;    
527             }
528         }
529
530         /* Generate next prime candidate:
531            p = 2 * q [ * q_factor] * factor_0 * factor_1 * ... * factor_n + 1. 
532          */
533         mpi_set (prime, q);
534         mpi_mul_ui (prime, prime, 2);
535         if (need_q_factor)
536           mpi_mul (prime, prime, q_factor);
537         for(i = 0; i < n; i++)
538           mpi_mul (prime, prime, factors[i]);
539         mpi_add_ui (prime, prime, 1);
540         nprime = mpi_get_nbits (prime);
541
542         if (nprime < pbits)
543           {
544             if (++count1 > 20)
545               {
546                 count1 = 0;
547                 qbits++;
548                 progress('>');
549                 mpi_free (q);
550                 q = gen_prime (qbits, is_secret, randomlevel, NULL, NULL);
551                 goto next_try;
552               }
553           }
554         else
555           count1 = 0;
556         
557         if (nprime > pbits)
558           {
559             if (++count2 > 20)
560               {
561                 count2 = 0;
562                 qbits--;
563                 progress('<');
564                 mpi_free (q);
565                 q = gen_prime (qbits, is_secret, randomlevel, NULL, NULL);
566                 goto next_try;
567               }
568           }
569         else
570           count2 = 0;
571     }
572   while (! ((nprime == pbits) && check_prime (prime, val_2, 5,
573                                               cb_func, cb_arg)));
574
575   if (DBG_CIPHER)
576     {
577       progress ('\n');
578       log_mpidump ("prime    : ", prime);
579       log_mpidump ("factor  q: ", q);
580       if (need_q_factor)
581         log_mpidump ("factor q0: ", q_factor);
582       for (i = 0; i < n; i++)
583         log_mpidump ("factor pi: ", factors[i]);
584       log_debug ("bit sizes: prime=%u, q=%u",
585                  mpi_get_nbits (prime), mpi_get_nbits (q));
586       if (need_q_factor)
587         log_debug (", q0=%u", mpi_get_nbits (q_factor));
588       for (i = 0; i < n; i++)
589         log_debug (", p%d=%u", i, mpi_get_nbits (factors[i]));
590       progress('\n');
591     }
592
593   if (ret_factors)
594     {
595       /* Caller wants the factors.  */
596       factors_new = gcry_calloc (n + 4, sizeof (*factors_new));
597       if (! factors_new)
598         {
599           err = gpg_err_code_from_errno (errno);
600           goto leave;
601         }
602
603       if (all_factors)
604         {
605           i = 0;
606           factors_new[i++] = gcry_mpi_set_ui (NULL, 2);
607           factors_new[i++] = mpi_copy (q);
608           if (need_q_factor)
609             factors_new[i++] = mpi_copy (q_factor);
610           for(j=0; j < n; j++)
611             factors_new[i++] = mpi_copy (factors[j]);
612         }
613       else
614         {
615           i = 0;
616           if (need_q_factor)
617             {
618               factors_new[i++] = mpi_copy (q_factor);
619               for (; i <= n; i++)
620                 factors_new[i] = mpi_copy (factors[i]);
621             }
622           else
623             for (; i < n; i++ )
624               factors_new[i] = mpi_copy (factors[i]);
625         }
626     }
627   
628   if (g)
629     {
630       /* Create a generator (start with 3).  */
631       gcry_mpi_t tmp = mpi_alloc (mpi_get_nlimbs (prime));
632       gcry_mpi_t b = mpi_alloc (mpi_get_nlimbs (prime));
633       gcry_mpi_t pmin1 = mpi_alloc (mpi_get_nlimbs (prime));
634       
635       if (need_q_factor)
636         err = GPG_ERR_NOT_IMPLEMENTED;
637       else
638         {
639           factors[n] = q;
640           factors[n + 1] = mpi_alloc_set_ui (2);
641           mpi_sub_ui (pmin1, prime, 1);
642           mpi_set_ui (g, 2);
643           do
644             {
645               mpi_add_ui (g, g, 1);
646               if (DBG_CIPHER)
647                 {
648                   log_debug ("checking g:");
649                   gcry_mpi_dump (g);
650                   log_printf ("\n");
651                 }
652               else
653                 progress('^');
654               for (i = 0; i < n + 2; i++)
655                 {
656                   mpi_fdiv_q (tmp, pmin1, factors[i]);
657                   /* No mpi_pow(), but it is okay to use this with mod
658                      prime.  */
659                   gcry_mpi_powm (b, g, tmp, prime);
660                   if (! mpi_cmp_ui (b, 1))
661                     break;
662                 }
663               if (DBG_CIPHER)
664                 progress('\n');
665             } 
666           while (i < n + 2);
667
668           mpi_free (factors[n+1]);
669           mpi_free (tmp);
670           mpi_free (b);
671           mpi_free (pmin1);
672         }
673     }
674   
675   if (! DBG_CIPHER)
676     progress ('\n');
677
678
679  leave:
680   if (pool)
681     {
682       is_locked = !ath_mutex_lock (&primepool_lock);
683       for(i = 0; i < m; i++)
684         {
685           if (pool[i])
686             {
687               for (j=0; j < n; j++)
688                 if (pool_in_use[j] == i)
689                   break;
690               if (j == n && is_locked)
691                 {
692                   /* This pooled subprime has not been used. */
693                   save_pool_prime (pool[i], poolrandomlevel);
694                 }
695               else
696                 mpi_free (pool[i]);
697             }
698         }
699       if (is_locked && ath_mutex_unlock (&primepool_lock))
700         err = GPG_ERR_INTERNAL;
701       is_locked = 0;
702       gcry_free (pool);
703     }
704   gcry_free (pool_in_use);
705   if (factors)
706     gcry_free (factors);  /* Factors are shallow copies.  */
707   if (perms)
708     gcry_free (perms);
709
710   mpi_free (val_2);
711   mpi_free (q);
712   mpi_free (q_factor);
713
714   if (! err)
715     {
716       *prime_generated = prime;
717       if (ret_factors)
718         *ret_factors = factors_new;
719     }
720   else
721     {
722       if (factors_new)
723         {
724           for (i = 0; factors_new[i]; i++)
725             mpi_free (factors_new[i]);
726           gcry_free (factors_new);
727         }
728       mpi_free (prime);
729     }
730
731   return err;
732 }
733
734
735 /* Generate a prime used for discrete logarithm algorithms; i.e. this
736    prime will be public and no strong random is required.  */
737 gcry_mpi_t
738 _gcry_generate_elg_prime (int mode, unsigned pbits, unsigned qbits,
739                           gcry_mpi_t g, gcry_mpi_t **ret_factors)
740 {
741   gcry_err_code_t err = GPG_ERR_NO_ERROR;
742   gcry_mpi_t prime = NULL;
743   
744   err = prime_generate_internal ((mode == 1), &prime, pbits, qbits, g,
745                                  ret_factors, GCRY_WEAK_RANDOM, 0, 0,
746                                  NULL, NULL);
747
748   return prime;
749 }
750
751
752 static gcry_mpi_t
753 gen_prime (unsigned int nbits, int secret, int randomlevel, 
754            int (*extra_check)(void *, gcry_mpi_t), void *extra_check_arg)
755 {
756   gcry_mpi_t prime, ptest, pminus1, val_2, val_3, result;
757   int i;
758   unsigned int x, step;
759   unsigned int count1, count2;
760   int *mods;
761   
762 /*   if (  DBG_CIPHER ) */
763 /*     log_debug ("generate a prime of %u bits ", nbits ); */
764
765   if (nbits < 16)
766     log_fatal ("can't generate a prime with less than %d bits\n", 16);
767
768   mods = gcry_xmalloc( no_of_small_prime_numbers * sizeof *mods );
769   /* Make nbits fit into gcry_mpi_t implementation. */
770   val_2  = mpi_alloc_set_ui( 2 );
771   val_3 = mpi_alloc_set_ui( 3);
772   prime  = secret? gcry_mpi_snew ( nbits ): gcry_mpi_new ( nbits );
773   result = mpi_alloc_like( prime );
774   pminus1= mpi_alloc_like( prime );
775   ptest  = mpi_alloc_like( prime );
776   count1 = count2 = 0;
777   for (;;)
778     {  /* try forvever */
779       int dotcount=0;
780       
781       /* generate a random number */
782       gcry_mpi_randomize( prime, nbits, randomlevel );
783       
784       /* Set high order bit to 1, set low order bit to 1.  If we are
785          generating a secret prime we are most probably doing that
786          for RSA, to make sure that the modulus does have the
787          requested key size we set the 2 high order bits. */
788       mpi_set_highbit (prime, nbits-1);
789       if (secret)
790         mpi_set_bit (prime, nbits-2);
791       mpi_set_bit(prime, 0);
792       
793       /* Calculate all remainders. */
794       for (i=0; (x = small_prime_numbers[i]); i++ )
795         mods[i] = mpi_fdiv_r_ui(NULL, prime, x);
796       
797       /* Now try some primes starting with prime. */
798       for(step=0; step < 20000; step += 2 ) 
799         {
800           /* Check against all the small primes we have in mods. */
801           count1++;
802           for (i=0; (x = small_prime_numbers[i]); i++ ) 
803             {
804               while ( mods[i] + step >= x )
805                 mods[i] -= x;
806               if ( !(mods[i] + step) )
807                 break;
808             }
809           if ( x )
810             continue;   /* Found a multiple of an already known prime. */
811           
812           mpi_add_ui( ptest, prime, step );
813
814           /* Do a fast Fermat test now. */
815           count2++;
816           mpi_sub_ui( pminus1, ptest, 1);
817           gcry_mpi_powm( result, val_2, pminus1, ptest );
818           if ( !mpi_cmp_ui( result, 1 ) )
819             { 
820               /* Not composite, perform stronger tests */
821               if (is_prime(ptest, 5, &count2 ))
822                 {
823                   if (!mpi_test_bit( ptest, nbits-1-secret ))
824                     {
825                       progress('\n');
826                       log_debug ("overflow in prime generation\n");
827                       break; /* Stop loop, continue with a new prime. */
828                     }
829
830                   if (extra_check && extra_check (extra_check_arg, ptest))
831                     { 
832                       /* The extra check told us that this prime is
833                          not of the caller's taste. */
834                       progress ('/');
835                     }
836                   else
837                     { 
838                       /* Got it. */
839                       mpi_free(val_2);
840                       mpi_free(val_3);
841                       mpi_free(result);
842                       mpi_free(pminus1);
843                       mpi_free(prime);
844                       gcry_free(mods);
845                       return ptest; 
846                     }
847                 }
848             }
849           if (++dotcount == 10 )
850             {
851               progress('.');
852               dotcount = 0;
853             }
854         }
855       progress(':'); /* restart with a new random value */
856     }
857 }
858
859 /****************
860  * Returns: true if this may be a prime
861  * RM_ROUNDS gives the number of Rabin-Miller tests to run.
862  */
863 static int
864 check_prime( gcry_mpi_t prime, gcry_mpi_t val_2, int rm_rounds,
865              gcry_prime_check_func_t cb_func, void *cb_arg)
866 {
867   int i;
868   unsigned int x;
869   unsigned int count=0;
870
871   /* Check against small primes. */
872   for (i=0; (x = small_prime_numbers[i]); i++ )
873     {
874       if ( mpi_divisible_ui( prime, x ) )
875         return 0;
876     }
877
878   /* A quick Fermat test. */
879   {
880     gcry_mpi_t result = mpi_alloc_like( prime );
881     gcry_mpi_t pminus1 = mpi_alloc_like( prime );
882     mpi_sub_ui( pminus1, prime, 1);
883     gcry_mpi_powm( result, val_2, pminus1, prime );
884     mpi_free( pminus1 );
885     if ( mpi_cmp_ui( result, 1 ) )
886       { 
887         /* Is composite. */
888         mpi_free( result );
889         progress('.');
890         return 0;
891       }
892     mpi_free( result );
893   }
894
895   if (!cb_func || cb_func (cb_arg, GCRY_PRIME_CHECK_AT_MAYBE_PRIME, prime))
896     {
897       /* Perform stronger tests. */
898       if ( is_prime( prime, rm_rounds, &count ) )
899         {
900           if (!cb_func
901               || cb_func (cb_arg, GCRY_PRIME_CHECK_AT_GOT_PRIME, prime))
902             return 1; /* Probably a prime. */
903         }
904     }
905   progress('.');
906   return 0;
907 }
908
909
910 /*
911  * Return true if n is probably a prime
912  */
913 static int
914 is_prime (gcry_mpi_t n, int steps, unsigned int *count)
915 {
916   gcry_mpi_t x = mpi_alloc( mpi_get_nlimbs( n ) );
917   gcry_mpi_t y = mpi_alloc( mpi_get_nlimbs( n ) );
918   gcry_mpi_t z = mpi_alloc( mpi_get_nlimbs( n ) );
919   gcry_mpi_t nminus1 = mpi_alloc( mpi_get_nlimbs( n ) );
920   gcry_mpi_t a2 = mpi_alloc_set_ui( 2 );
921   gcry_mpi_t q;
922   unsigned i, j, k;
923   int rc = 0;
924   unsigned nbits = mpi_get_nbits( n );
925
926   if (steps < 5) /* Make sure that we do at least 5 rounds. */
927     steps = 5; 
928
929   mpi_sub_ui( nminus1, n, 1 );
930
931   /* Find q and k, so that n = 1 + 2^k * q . */
932   q = mpi_copy ( nminus1 );
933   k = mpi_trailing_zeros ( q );
934   mpi_tdiv_q_2exp (q, q, k);
935
936   for (i=0 ; i < steps; i++ )
937     {
938       ++*count;
939       if( !i )
940         {
941           mpi_set_ui( x, 2 );
942         }
943       else
944         {
945           gcry_mpi_randomize( x, nbits, GCRY_WEAK_RANDOM );
946
947           /* Make sure that the number is smaller than the prime and
948              keep the randomness of the high bit. */
949           if ( mpi_test_bit ( x, nbits-2) )
950             {
951               mpi_set_highbit ( x, nbits-2); /* Clear all higher bits. */
952             }
953           else
954             {
955               mpi_set_highbit( x, nbits-2 );
956               mpi_clear_bit( x, nbits-2 );
957             }
958           gcry_assert (mpi_cmp (x, nminus1) < 0 && mpi_cmp_ui (x, 1) > 0);
959         }
960       gcry_mpi_powm ( y, x, q, n);
961       if ( mpi_cmp_ui(y, 1) && mpi_cmp( y, nminus1 ) )
962         {
963           for ( j=1; j < k && mpi_cmp( y, nminus1 ); j++ )
964             {
965               gcry_mpi_powm(y, y, a2, n);
966               if( !mpi_cmp_ui( y, 1 ) )
967                 goto leave; /* Not a prime. */
968             }
969           if (mpi_cmp( y, nminus1 ) )
970             goto leave; /* Not a prime. */
971         }
972       progress('+');
973     }
974   rc = 1; /* May be a prime. */
975
976  leave:
977   mpi_free( x );
978   mpi_free( y );
979   mpi_free( z );
980   mpi_free( nminus1 );
981   mpi_free( q );
982   mpi_free( a2 );
983
984   return rc;
985 }
986
987
988 /* Given ARRAY of size N with M elements set to true produce a
989    modified array with the next permutation of M elements.  Note, that
990    ARRAY is used in a one-bit-per-byte approach.  To detected the last
991    permutation it is useful to intialize the array with the first M
992    element set to true and use this test:
993        m_out_of_n (array, m, n);
994        for (i = j = 0; i < n && j < m; i++)
995          if (array[i])
996            j++;
997        if (j == m)
998          goto ready;
999      
1000    This code is based on the algorithm 452 from the "Collected
1001    Algorithms From ACM, Volume II" by C. N. Liu and D. T. Tang.
1002 */
1003 static void
1004 m_out_of_n ( char *array, int m, int n )
1005 {
1006   int i=0, i1=0, j=0, jp=0,  j1=0, k1=0, k2=0;
1007
1008   if( !m || m >= n )
1009     return;
1010
1011   /* Need to handle this simple case separately. */
1012   if( m == 1 )
1013     { 
1014       for (i=0; i < n; i++ )
1015         {
1016           if ( array[i] )
1017             {
1018               array[i++] = 0;
1019               if( i >= n )
1020                 i = 0;
1021               array[i] = 1;
1022               return;
1023             }
1024         }
1025       BUG();
1026     }
1027
1028
1029   for (j=1; j < n; j++ )
1030     {
1031       if ( array[n-1] == array[n-j-1])
1032         continue;
1033       j1 = j;
1034       break;
1035     }
1036
1037   if ( (m & 1) )
1038     {
1039       /* M is odd. */
1040       if( array[n-1] )
1041         {
1042           if( j1 & 1 )
1043             {
1044               k1 = n - j1;
1045               k2 = k1+2;
1046               if( k2 > n )
1047                 k2 = n;
1048               goto leave;
1049             }
1050           goto scan;
1051         }
1052       k2 = n - j1 - 1;
1053       if( k2 == 0 )
1054         {
1055           k1 = i;
1056           k2 = n - j1;
1057         }
1058       else if( array[k2] && array[k2-1] )
1059         k1 = n;
1060       else
1061         k1 = k2 + 1;
1062     }
1063   else 
1064     {
1065       /* M is even. */
1066       if( !array[n-1] )
1067         {
1068           k1 = n - j1;
1069           k2 = k1 + 1;
1070           goto leave;
1071         }
1072         
1073       if( !(j1 & 1) )
1074         {
1075           k1 = n - j1;
1076           k2 = k1+2;
1077           if( k2 > n )
1078             k2 = n;
1079           goto leave;
1080         }
1081     scan:
1082       jp = n - j1 - 1;
1083       for (i=1; i <= jp; i++ ) 
1084         {
1085           i1 = jp + 2 - i;
1086           if( array[i1-1]  )
1087             {
1088               if( array[i1-2] )
1089                 {
1090                   k1 = i1 - 1;
1091                   k2 = n - j1;
1092                 }
1093               else
1094                 {
1095                   k1 = i1 - 1;
1096                   k2 = n + 1 - j1;
1097                 }
1098               goto leave;
1099             }
1100         }
1101       k1 = 1;
1102       k2 = n + 1 - m;
1103     }
1104  leave:
1105   /* Now complement the two selected bits. */
1106   array[k1-1] = !array[k1-1];
1107   array[k2-1] = !array[k2-1];
1108 }
1109
1110
1111 /* Generate a new prime number of PRIME_BITS bits and store it in
1112    PRIME.  If FACTOR_BITS is non-zero, one of the prime factors of
1113    (prime - 1) / 2 must be FACTOR_BITS bits long.  If FACTORS is
1114    non-zero, allocate a new, NULL-terminated array holding the prime
1115    factors and store it in FACTORS.  FLAGS might be used to influence
1116    the prime number generation process.  */
1117 gcry_error_t
1118 gcry_prime_generate (gcry_mpi_t *prime, unsigned int prime_bits,
1119                      unsigned int factor_bits, gcry_mpi_t **factors,
1120                      gcry_prime_check_func_t cb_func, void *cb_arg,
1121                      gcry_random_level_t random_level,
1122                      unsigned int flags)
1123 {
1124   gcry_err_code_t err = GPG_ERR_NO_ERROR;
1125   gcry_mpi_t *factors_generated = NULL;
1126   gcry_mpi_t prime_generated = NULL;
1127   unsigned int mode = 0;
1128
1129   if (!prime)
1130     return gpg_error (GPG_ERR_INV_ARG);
1131   *prime = NULL; 
1132
1133   if (flags & GCRY_PRIME_FLAG_SPECIAL_FACTOR)
1134     mode = 1;
1135
1136   /* Generate.  */
1137   err = prime_generate_internal ((mode==1), &prime_generated, prime_bits,
1138                                  factor_bits, NULL,
1139                                  factors? &factors_generated : NULL,
1140                                  random_level, flags, 1,
1141                                  cb_func, cb_arg);
1142
1143   if (! err)
1144     if (cb_func)
1145       {
1146         /* Additional check. */
1147         if ( !cb_func (cb_arg, GCRY_PRIME_CHECK_AT_FINISH, prime_generated))
1148           {
1149             /* Failed, deallocate resources.  */
1150             unsigned int i;
1151
1152             mpi_free (prime_generated);
1153             if (factors)
1154               {
1155                 for (i = 0; factors_generated[i]; i++)
1156                   mpi_free (factors_generated[i]);
1157                 gcry_free (factors_generated);
1158               }
1159             err = GPG_ERR_GENERAL; 
1160           }
1161       }
1162
1163   if (! err)
1164     {
1165       if (factors)
1166         *factors = factors_generated;
1167       *prime = prime_generated;
1168     }
1169
1170   return gcry_error (err);
1171 }
1172
1173 /* Check wether the number X is prime.  */
1174 gcry_error_t
1175 gcry_prime_check (gcry_mpi_t x, unsigned int flags)
1176 {
1177   gcry_err_code_t err = GPG_ERR_NO_ERROR;
1178   gcry_mpi_t val_2 = mpi_alloc_set_ui (2); /* Used by the Fermat test. */
1179
1180   (void)flags;
1181
1182   /* We use 64 rounds because the prime we are going to test is not
1183      guaranteed to be a random one. */
1184   if (! check_prime (x, val_2, 64, NULL, NULL))
1185     err = GPG_ERR_NO_PRIME;
1186
1187   mpi_free (val_2);
1188
1189   return gcry_error (err);
1190 }
1191
1192 /* Find a generator for PRIME where the factorization of (prime-1) is
1193    in the NULL terminated array FACTORS. Return the generator as a
1194    newly allocated MPI in R_G.  If START_G is not NULL, use this as s
1195    atart for the search. Returns 0 on success.*/
1196 gcry_error_t
1197 gcry_prime_group_generator (gcry_mpi_t *r_g,
1198                             gcry_mpi_t prime, gcry_mpi_t *factors,
1199                             gcry_mpi_t start_g)
1200 {
1201   gcry_mpi_t tmp = gcry_mpi_new (0);
1202   gcry_mpi_t b = gcry_mpi_new (0);
1203   gcry_mpi_t pmin1 = gcry_mpi_new (0);
1204   gcry_mpi_t g = start_g? gcry_mpi_copy (start_g) : gcry_mpi_set_ui (NULL, 3);
1205   int first = 1;
1206   int i, n;
1207
1208   if (!factors || !r_g || !prime)
1209     return gpg_error (GPG_ERR_INV_ARG);
1210   *r_g = NULL; 
1211
1212   for (n=0; factors[n]; n++)
1213     ;
1214   if (n < 2)
1215     return gpg_error (GPG_ERR_INV_ARG);
1216
1217   /* Extra sanity check - usually disabled. */  
1218 /*   mpi_set (tmp, factors[0]); */
1219 /*   for(i = 1; i < n; i++) */
1220 /*     mpi_mul (tmp, tmp, factors[i]); */
1221 /*   mpi_add_ui (tmp, tmp, 1); */
1222 /*   if (mpi_cmp (prime, tmp)) */
1223 /*     return gpg_error (GPG_ERR_INV_ARG); */
1224   
1225   gcry_mpi_sub_ui (pmin1, prime, 1);      
1226   do         
1227     {
1228       if (first)
1229         first = 0;
1230       else
1231         gcry_mpi_add_ui (g, g, 1);
1232       
1233       if (DBG_CIPHER)
1234         {
1235           log_debug ("checking g:");
1236           gcry_mpi_dump (g);
1237           log_debug ("\n");
1238         }
1239       else
1240         progress('^');
1241       
1242       for (i = 0; i < n; i++)
1243         {
1244           mpi_fdiv_q (tmp, pmin1, factors[i]);
1245           gcry_mpi_powm (b, g, tmp, prime);
1246           if (! mpi_cmp_ui (b, 1))
1247             break;
1248         }
1249       if (DBG_CIPHER)
1250         progress('\n');
1251     }
1252   while (i < n);
1253   
1254   gcry_mpi_release (tmp);
1255   gcry_mpi_release (b); 
1256   gcry_mpi_release (pmin1); 
1257   *r_g = g; 
1258
1259   return 0; 
1260 }
1261
1262 /* Convenience function to release the factors array. */
1263 void
1264 gcry_prime_release_factors (gcry_mpi_t *factors)
1265 {
1266   if (factors)
1267     {
1268       int i;
1269       
1270       for (i=0; factors[i]; i++)
1271         mpi_free (factors[i]);
1272       gcry_free (factors);
1273     }
1274 }
1275
1276
1277 /* Helper for _gcry_generate_x931_prime.  */
1278 static gcry_mpi_t
1279 find_x931_prime (const gcry_mpi_t pfirst)
1280 {
1281   gcry_mpi_t val_2 = mpi_alloc_set_ui (2); 
1282   gcry_mpi_t prime;
1283   
1284   prime = gcry_mpi_copy (pfirst);
1285   /* If P is even add 1.  */ 
1286   mpi_set_bit (prime, 0);
1287
1288   /* We use 64 Rabin-Miller rounds which is better and thus
1289      sufficient.  We do not have a Lucas test implementaion thus we
1290      can't do it in the X9.31 preferred way of running a few
1291      Rabin-Miller followed by one Lucas test.  */
1292   while ( !check_prime (prime, val_2, 64, NULL, NULL) )
1293     mpi_add_ui (prime, prime, 2);
1294
1295   mpi_free (val_2);
1296
1297   return prime;
1298 }
1299
1300
1301 /* Generate a prime using the algorithm from X9.31 appendix B.4. 
1302
1303    This function requires that the provided public exponent E is odd.
1304    XP, XP1 and XP2 are the seed values.  All values are mandatory.
1305
1306    On success the prime is returned.  If R_P1 or R_P2 are given the
1307    internal values P1 and P2 are saved at these addresses.  On error
1308    NULL is returned.  */
1309 gcry_mpi_t
1310 _gcry_derive_x931_prime (const gcry_mpi_t xp, 
1311                          const gcry_mpi_t xp1, const gcry_mpi_t xp2,
1312                          const gcry_mpi_t e,
1313                          gcry_mpi_t *r_p1, gcry_mpi_t *r_p2)
1314 {
1315   gcry_mpi_t p1, p2, p1p2, yp0;
1316
1317   if (!xp || !xp1 || !xp2)
1318     return NULL;
1319   if (!e || !mpi_test_bit (e, 0))
1320     return NULL;  /* We support only odd values for E.  */
1321
1322   p1 = find_x931_prime (xp1);
1323   p2 = find_x931_prime (xp2);
1324   p1p2 = mpi_alloc_like (xp);
1325   mpi_mul (p1p2, p1, p2);
1326
1327   {
1328     gcry_mpi_t r1, tmp;
1329   
1330     /* r1 = (p2^{-1} mod p1)p2 - (p1^{-1} mod p2) */
1331     tmp = mpi_alloc_like (p1);
1332     mpi_invm (tmp, p2, p1);
1333     mpi_mul (tmp, tmp, p2);
1334     r1 = tmp;
1335     
1336     tmp = mpi_alloc_like (p2);
1337     mpi_invm (tmp, p1, p2);
1338     mpi_mul (tmp, tmp, p1);
1339     mpi_sub (r1, r1, tmp);
1340
1341     /* Fixup a negative value.  */
1342     if (mpi_is_neg (r1)) 
1343       mpi_add (r1, r1, p1p2);
1344
1345     /* yp0 = xp + (r1 - xp mod p1*p2)  */
1346     yp0 = tmp; tmp = NULL;
1347     mpi_subm (yp0, r1, xp, p1p2);
1348     mpi_add (yp0, yp0, xp);
1349     mpi_free (r1);
1350
1351     /* Fixup a negative value.  */
1352     if (mpi_cmp (yp0, xp) < 0 ) 
1353       mpi_add (yp0, yp0, p1p2);
1354   }
1355
1356   /* yp0 is now the first integer greater than xp with p1 being a
1357      large prime factor of yp0-1 and p2 a large prime factor of yp0+1.  */
1358
1359   /* Note that the first example from X9.31 (D.1.1) which uses
1360        (Xq1 #1A5CF72EE770DE50CB09ACCEA9#)
1361        (Xq2 #134E4CAA16D2350A21D775C404#)
1362        (Xq  #CC1092495D867E64065DEE3E7955F2EBC7D47A2D
1363              7C9953388F97DDDC3E1CA19C35CA659EDC2FC325
1364              6D29C2627479C086A699A49C4C9CEE7EF7BD1B34
1365              321DE34A#))))
1366      returns an yp0 of
1367             #CC1092495D867E64065DEE3E7955F2EBC7D47A2D
1368              7C9953388F97DDDC3E1CA19C35CA659EDC2FC4E3
1369              BF20CB896EE37E098A906313271422162CB6C642
1370              75C1201F#
1371      and not
1372             #CC1092495D867E64065DEE3E7955F2EBC7D47A2D
1373              7C9953388F97DDDC3E1CA19C35CA659EDC2FC2E6
1374              C88FE299D52D78BE405A97E01FD71DD7819ECB91
1375              FA85A076#
1376      as stated in the standard.  This seems to be a bug in X9.31.
1377    */
1378
1379   {
1380     gcry_mpi_t val_2 = mpi_alloc_set_ui (2); 
1381     gcry_mpi_t gcdtmp = mpi_alloc_like (yp0);
1382     int gcdres;
1383   
1384     mpi_sub_ui (p1p2, p1p2, 1); /* Adjust for loop body.  */
1385     mpi_sub_ui (yp0, yp0, 1);   /* Ditto.  */
1386     for (;;)
1387       {
1388         gcdres = gcry_mpi_gcd (gcdtmp, e, yp0);
1389         mpi_add_ui (yp0, yp0, 1);
1390         if (!gcdres)
1391           progress ('/');  /* gcd (e, yp0-1) != 1  */
1392         else if (check_prime (yp0, val_2, 64, NULL, NULL))
1393           break; /* Found.  */
1394         /* We add p1p2-1 because yp0 is incremented after the gcd test.  */
1395         mpi_add (yp0, yp0, p1p2);
1396       }
1397     mpi_free (gcdtmp);
1398     mpi_free (val_2);
1399   }
1400
1401   mpi_free (p1p2);
1402
1403   progress('\n');
1404   if (r_p1)
1405     *r_p1 = p1;
1406   else
1407     mpi_free (p1);
1408   if (r_p2)
1409     *r_p2 = p2;
1410   else
1411     mpi_free (p2);
1412   return yp0;
1413 }