rndsort{Miller, Rabin} primality test.
[openssl.git] / crypto / bn / bn_prime.c
1 /* crypto/bn/bn_prime.c */
2 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
3  * All rights reserved.
4  *
5  * This package is an SSL implementation written
6  * by Eric Young (eay@cryptsoft.com).
7  * The implementation was written so as to conform with Netscapes SSL.
8  * 
9  * This library is free for commercial and non-commercial use as long as
10  * the following conditions are aheared to.  The following conditions
11  * apply to all code found in this distribution, be it the RC4, RSA,
12  * lhash, DES, etc., code; not just the SSL code.  The SSL documentation
13  * included with this distribution is covered by the same copyright terms
14  * except that the holder is Tim Hudson (tjh@cryptsoft.com).
15  * 
16  * Copyright remains Eric Young's, and as such any Copyright notices in
17  * the code are not to be removed.
18  * If this package is used in a product, Eric Young should be given attribution
19  * as the author of the parts of the library used.
20  * This can be in the form of a textual message at program startup or
21  * in documentation (online or textual) provided with the package.
22  * 
23  * Redistribution and use in source and binary forms, with or without
24  * modification, are permitted provided that the following conditions
25  * are met:
26  * 1. Redistributions of source code must retain the copyright
27  *    notice, this list of conditions and the following disclaimer.
28  * 2. Redistributions in binary form must reproduce the above copyright
29  *    notice, this list of conditions and the following disclaimer in the
30  *    documentation and/or other materials provided with the distribution.
31  * 3. All advertising materials mentioning features or use of this software
32  *    must display the following acknowledgement:
33  *    "This product includes cryptographic software written by
34  *     Eric Young (eay@cryptsoft.com)"
35  *    The word 'cryptographic' can be left out if the rouines from the library
36  *    being used are not cryptographic related :-).
37  * 4. If you include any Windows specific code (or a derivative thereof) from 
38  *    the apps directory (application code) you must include an acknowledgement:
39  *    "This product includes software written by Tim Hudson (tjh@cryptsoft.com)"
40  * 
41  * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND
42  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
43  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
44  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
45  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
46  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
47  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
48  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
49  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
50  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
51  * SUCH DAMAGE.
52  * 
53  * The licence and distribution terms for any publically available version or
54  * derivative of this code cannot be changed.  i.e. this code cannot simply be
55  * copied and put under another distribution licence
56  * [including the GNU Public Licence.]
57  */
58
59 #include <stdio.h>
60 #include <time.h>
61 #include "cryptlib.h"
62 #include "bn_lcl.h"
63 #include <openssl/rand.h>
64
65 /* The quick sieve algorithm approach to weeding out primes is
66  * Philip Zimmermann's, as implemented in PGP.  I have had a read of
67  * his comments and implemented my own version.
68  */
69 #include "bn_prime.h"
70
71 static int witness(BIGNUM *w, BIGNUM *a, BIGNUM *a1, BIGNUM *a1_odd, int k,
72         BN_CTX *ctx, BN_MONT_CTX *mont);
73 static int probable_prime(BIGNUM *rnd, int bits);
74 static int probable_prime_dh(BIGNUM *rnd, int bits,
75         BIGNUM *add, BIGNUM *rem, BN_CTX *ctx);
76 static int probable_prime_dh_safe(BIGNUM *rnd, int bits,
77         BIGNUM *add, BIGNUM *rem, BN_CTX *ctx);
78
79 BIGNUM *BN_generate_prime(BIGNUM *ret, int bits, int safe, BIGNUM *add,
80              BIGNUM *rem, void (*callback)(int,int,void *), void *cb_arg)
81         {
82         BIGNUM *rnd=NULL;
83         BIGNUM t;
84         int found=0;
85         int i,j,c1=0;
86         BN_CTX *ctx;
87         int checks = BN_prime_checks_for_size(bits);
88
89         ctx=BN_CTX_new();
90         if (ctx == NULL) goto err;
91         if (ret == NULL)
92                 {
93                 if ((rnd=BN_new()) == NULL) goto err;
94                 }
95         else
96                 rnd=ret;
97         BN_init(&t);
98 loop: 
99         /* make a random number and set the top and bottom bits */
100         if (add == NULL)
101                 {
102                 if (!probable_prime(rnd,bits)) goto err;
103                 }
104         else
105                 {
106                 if (safe)
107                         {
108                         if (!probable_prime_dh_safe(rnd,bits,add,rem,ctx))
109                                  goto err;
110                         }
111                 else
112                         {
113                         if (!probable_prime_dh(rnd,bits,add,rem,ctx))
114                                 goto err;
115                         }
116                 }
117         /* if (BN_mod_word(rnd,(BN_ULONG)3) == 1) goto loop; */
118         if (callback != NULL) callback(0,c1++,cb_arg);
119
120         if (!safe)
121                 {
122                 i=BN_is_prime_fasttest(rnd,checks,callback,ctx,cb_arg,0);
123                 if (i == -1) goto err;
124                 if (i == 0) goto loop;
125                 }
126         else
127                 {
128                 /* for "safe prime" generation,
129                  * check that (p-1)/2 is prime.
130                  * Since a prime is odd, We just
131                  * need to divide by 2 */
132                 if (!BN_rshift1(&t,rnd)) goto err;
133
134                 for (i=0; i<checks; i++)
135                         {
136                         j=BN_is_prime_fasttest(rnd,1,callback,ctx,cb_arg,0);
137                         if (j == -1) goto err;
138                         if (j == 0) goto loop;
139
140                         j=BN_is_prime_fasttest(&t,1,callback,ctx,cb_arg,0);
141                         if (j == -1) goto err;
142                         if (j == 0) goto loop;
143
144                         if (callback != NULL) callback(2,c1-1,cb_arg);
145                         /* We have a safe prime test pass */
146                         }
147                 }
148         /* we have a prime :-) */
149         found = 1;
150 err:
151         if (!found && (ret == NULL) && (rnd != NULL)) BN_free(rnd);
152         BN_free(&t);
153         if (ctx != NULL) BN_CTX_free(ctx);
154         return(found ? rnd : NULL);
155         }
156
157 int BN_is_prime(BIGNUM *a, int checks, void (*callback)(int,int,void *),
158         BN_CTX *ctx_passed, void *cb_arg)
159         {
160         return BN_is_prime_fasttest(a, checks, callback, ctx_passed, cb_arg, 0);
161         }
162
163 int BN_is_prime_fasttest(BIGNUM *a, int checks,
164                 void (*callback)(int,int,void *),
165                 BN_CTX *ctx_passed, void *cb_arg,
166                 int do_trial_division)
167         {
168         int i, j, ret = -1;
169         int k;
170         BN_CTX *ctx = NULL;
171         BIGNUM *a1, *a1_odd, *check; /* taken from ctx */
172         BN_MONT_CTX *mont = NULL;
173
174         if (checks == BN_prime_checks)
175                 checks = BN_prime_checks_for_size(BN_num_bits(a));
176
177         if (a->neg) /* for now, refuse to handle negative numbers */
178                 return -1;
179
180         /* first look for small factors */
181         if (!BN_is_odd(a))
182                 return(0);
183         if (do_trial_division)
184                 {
185                 for (i = 1; i < NUMPRIMES; i++)
186                         if (BN_mod_word(a, primes[i]) == 0) 
187                                 return 0;
188                 if (callback != NULL) callback(1, -1, cb_arg);
189                 }
190
191         if (ctx_passed != NULL)
192                 ctx = ctx_passed;
193         else
194                 if ((ctx=BN_CTX_new()) == NULL)
195                         goto err;
196         a1 = &(ctx->bn[ctx->tos++]);
197         a1_odd = &(ctx->bn[ctx->tos++]);
198         check = &(ctx->bn[ctx->tos++]);;
199
200         /* compute a1 := a - 1 */
201         if (!BN_copy(a1, a))
202                 goto err;
203         if (!BN_sub_word(a1, 1))
204                 goto err;
205         if (BN_is_zero(a1))
206                 {
207                 ret = 0;
208                 goto err;
209                 }
210
211         /* write  a1  as  a1_odd * 2^k */
212         k = 1;
213         while (!BN_is_bit_set(a1, k))
214                 k++;
215         if (!BN_rshift(a1_odd, a1, k))
216                 goto err;
217
218         /* Montgomery setup for computations mod a */
219         mont = BN_MONT_CTX_new();
220         if (mont == NULL)
221                 goto err;
222         if (!BN_MONT_CTX_set(mont, a, ctx))
223                 goto err;
224         
225         for (i = 0; i < checks; i++)
226                 {
227                 if (!BN_pseudo_rand(check, BN_num_bits(a1), 0, 0))
228                         goto err;
229                 if (BN_cmp(check, a1) >= 0)
230                         if (!BN_sub(check, check, a1))
231                                 goto err;
232                 if (!BN_add_word(check, 1))
233                         goto err;
234                 /* now 1 <= check < a */
235
236                 j = witness(check, a, a1, a1_odd, k, ctx, mont);
237                 if (j == -1) goto err;
238                 if (j)
239                         {
240                         ret=0;
241                         goto err;
242                         }
243                 if (callback != NULL) callback(1,i,cb_arg);
244                 }
245         ret=1;
246 err:
247         if (ctx_passed != NULL)
248                 ctx_passed->tos -= 3; /* a1, a1_odd, check */
249         else if (ctx != NULL)
250                 BN_CTX_free(ctx);
251         if (mont != NULL)
252                 BN_MONT_CTX_free(mont);
253
254         return(ret);
255         }
256
257 static int witness(BIGNUM *w, BIGNUM *a, BIGNUM *a1, BIGNUM *a1_odd, int k,
258         BN_CTX *ctx, BN_MONT_CTX *mont)
259         {
260         if (!BN_mod_exp_mont(w, w, a1_odd, a, ctx, mont)) /* w := w^a1_odd mod a */
261                 return -1;
262         if (BN_is_one(w))
263                 return 0; /* probably prime */
264         if (BN_cmp(w, a1) == 0)
265                 return 0; /* w == -1 (mod a),  'a' is probably prime */
266         while (--k)
267                 {
268                 if (!BN_mod_mul(w, w, w, a, ctx)) /* w := w^2 mod a */
269                         return -1;
270                 if (BN_is_one(w))
271                         return 1; /* 'a' is composite, otherwise a previous 'w' would
272                                    * have been == -1 (mod 'a') */
273                 if (BN_cmp(w, a1) == 0)
274                         return 0; /* w == -1 (mod a), 'a' is probably prime */
275                 }
276         /* If we get here, 'w' is the (a-1)/2-th power of the original 'w',
277          * and it is neither -1 nor +1 -- so 'a' cannot be prime */
278         return 1;
279         }
280
281 static int probable_prime(BIGNUM *rnd, int bits)
282         {
283         int i;
284         BN_ULONG mods[NUMPRIMES];
285         BN_ULONG delta,d;
286
287 again:
288         if (!BN_rand(rnd,bits,1,1)) return(0);
289         /* we now have a random number 'rand' to test. */
290         for (i=1; i<NUMPRIMES; i++)
291                 mods[i]=BN_mod_word(rnd,(BN_ULONG)primes[i]);
292         delta=0;
293         loop: for (i=1; i<NUMPRIMES; i++)
294                 {
295                 /* check that rnd is not a prime and also
296                  * that gcd(rnd-1,primes) == 1 (except for 2) */
297                 if (((mods[i]+delta)%primes[i]) <= 1)
298                         {
299                         d=delta;
300                         delta+=2;
301                         /* perhaps need to check for overflow of
302                          * delta (but delta can be upto 2^32)
303                          * 21-May-98 eay - added overflow check */
304                         if (delta < d) goto again;
305                         goto loop;
306                         }
307                 }
308         if (!BN_add_word(rnd,delta)) return(0);
309         return(1);
310         }
311
312 static int probable_prime_dh(BIGNUM *rnd, int bits, BIGNUM *add, BIGNUM *rem,
313              BN_CTX *ctx)
314         {
315         int i,ret=0;
316         BIGNUM *t1;
317
318         t1= &(ctx->bn[ctx->tos++]);
319
320         if (!BN_rand(rnd,bits,0,1)) goto err;
321
322         /* we need ((rnd-rem) % add) == 0 */
323
324         if (!BN_mod(t1,rnd,add,ctx)) goto err;
325         if (!BN_sub(rnd,rnd,t1)) goto err;
326         if (rem == NULL)
327                 { if (!BN_add_word(rnd,1)) goto err; }
328         else
329                 { if (!BN_add(rnd,rnd,rem)) goto err; }
330
331         /* we now have a random number 'rand' to test. */
332
333         loop: for (i=1; i<NUMPRIMES; i++)
334                 {
335                 /* check that rnd is a prime */
336                 if (BN_mod_word(rnd,(BN_ULONG)primes[i]) <= 1)
337                         {
338                         if (!BN_add(rnd,rnd,add)) goto err;
339                         goto loop;
340                         }
341                 }
342         ret=1;
343 err:
344         ctx->tos--;
345         return(ret);
346         }
347
348 static int probable_prime_dh_safe(BIGNUM *p, int bits, BIGNUM *padd,
349              BIGNUM *rem, BN_CTX *ctx)
350         {
351         int i,ret=0;
352         BIGNUM *t1,*qadd=NULL,*q=NULL;
353
354         bits--;
355         t1= &(ctx->bn[ctx->tos++]);
356         q= &(ctx->bn[ctx->tos++]);
357         qadd= &(ctx->bn[ctx->tos++]);
358
359         if (!BN_rshift1(qadd,padd)) goto err;
360                 
361         if (!BN_rand(q,bits,0,1)) goto err;
362
363         /* we need ((rnd-rem) % add) == 0 */
364         if (!BN_mod(t1,q,qadd,ctx)) goto err;
365         if (!BN_sub(q,q,t1)) goto err;
366         if (rem == NULL)
367                 { if (!BN_add_word(q,1)) goto err; }
368         else
369                 {
370                 if (!BN_rshift1(t1,rem)) goto err;
371                 if (!BN_add(q,q,t1)) goto err;
372                 }
373
374         /* we now have a random number 'rand' to test. */
375         if (!BN_lshift1(p,q)) goto err;
376         if (!BN_add_word(p,1)) goto err;
377
378         loop: for (i=1; i<NUMPRIMES; i++)
379                 {
380                 /* check that p and q are prime */
381                 /* check that for p and q
382                  * gcd(p-1,primes) == 1 (except for 2) */
383                 if (    (BN_mod_word(p,(BN_ULONG)primes[i]) == 0) ||
384                         (BN_mod_word(q,(BN_ULONG)primes[i]) == 0))
385                         {
386                         if (!BN_add(p,p,padd)) goto err;
387                         if (!BN_add(q,q,qadd)) goto err;
388                         goto loop;
389                         }
390                 }
391         ret=1;
392 err:
393         ctx->tos-=3;
394         return(ret);
395         }