1 0fc65b37 2004-03-21 devnull #include "os.h"
 2 0fc65b37 2004-03-21 devnull #include <mp.h>
 3 0fc65b37 2004-03-21 devnull #include <libsec.h>
 4 0fc65b37 2004-03-21 devnull 
 5 cbeb0b26 2006-04-01 devnull /* Miller-Rabin probabilistic primality testing */
 6 cbeb0b26 2006-04-01 devnull /*	Knuth (1981) Seminumerical Algorithms, p.379 */
 7 cbeb0b26 2006-04-01 devnull /*	Menezes et al () Handbook, p.39 */
 8 cbeb0b26 2006-04-01 devnull /* 0 if composite; 1 if almost surely prime, Pr(err)<1/4**nrep */
 9 0fc65b37 2004-03-21 devnull int
10 0fc65b37 2004-03-21 devnull probably_prime(mpint *n, int nrep)
11 0fc65b37 2004-03-21 devnull {
12 1619f52c 2010-03-12 rsc 	int j, k, rep, nbits, isprime;
13 0fc65b37 2004-03-21 devnull 	mpint *nm1, *q, *x, *y, *r;
14 0fc65b37 2004-03-21 devnull 
15 0fc65b37 2004-03-21 devnull 	if(n->sign < 0)
16 0fc65b37 2004-03-21 devnull 		sysfatal("negative prime candidate");
17 0fc65b37 2004-03-21 devnull 
18 0fc65b37 2004-03-21 devnull 	if(nrep <= 0)
19 0fc65b37 2004-03-21 devnull 		nrep = 18;
20 0fc65b37 2004-03-21 devnull 
21 0fc65b37 2004-03-21 devnull 	k = mptoi(n);
22 cbeb0b26 2006-04-01 devnull 	if(k == 2)		/* 2 is prime */
23 0fc65b37 2004-03-21 devnull 		return 1;
24 cbeb0b26 2006-04-01 devnull 	if(k < 2)		/* 1 is not prime */
25 0fc65b37 2004-03-21 devnull 		return 0;
26 cbeb0b26 2006-04-01 devnull 	if((n->p[0] & 1) == 0)	/* even is not prime */
27 0fc65b37 2004-03-21 devnull 		return 0;
28 0fc65b37 2004-03-21 devnull 
29 cbeb0b26 2006-04-01 devnull 	/* test against small prime numbers */
30 0fc65b37 2004-03-21 devnull 	if(smallprimetest(n) < 0)
31 0fc65b37 2004-03-21 devnull 		return 0;
32 0fc65b37 2004-03-21 devnull 
33 cbeb0b26 2006-04-01 devnull 	/* fermat test, 2^n mod n == 2 if p is prime */
34 0fc65b37 2004-03-21 devnull 	x = uitomp(2, nil);
35 0fc65b37 2004-03-21 devnull 	y = mpnew(0);
36 0fc65b37 2004-03-21 devnull 	mpexp(x, n, n, y);
37 0fc65b37 2004-03-21 devnull 	k = mptoi(y);
38 0fc65b37 2004-03-21 devnull 	if(k != 2){
39 0fc65b37 2004-03-21 devnull 		mpfree(x);
40 0fc65b37 2004-03-21 devnull 		mpfree(y);
41 0fc65b37 2004-03-21 devnull 		return 0;
42 0fc65b37 2004-03-21 devnull 	}
43 0fc65b37 2004-03-21 devnull 
44 0fc65b37 2004-03-21 devnull 	nbits = mpsignif(n);
45 0fc65b37 2004-03-21 devnull 	nm1 = mpnew(nbits);
46 cbeb0b26 2006-04-01 devnull 	mpsub(n, mpone, nm1);	/* nm1 = n - 1 */
47 0fc65b37 2004-03-21 devnull 	k = mplowbits0(nm1);
48 0fc65b37 2004-03-21 devnull 	q = mpnew(0);
49 cbeb0b26 2006-04-01 devnull 	mpright(nm1, k, q);	/* q = (n-1)/2**k */
50 0fc65b37 2004-03-21 devnull 
51 0fc65b37 2004-03-21 devnull 	for(rep = 0; rep < nrep; rep++){
52 1619f52c 2010-03-12 rsc 		for(;;){
53 1619f52c 2010-03-12 rsc 			/* find x = random in [2, n-2] */
54 1619f52c 2010-03-12 rsc 			r = mprand(nbits, prng, nil);
55 1619f52c 2010-03-12 rsc 			mpmod(r, nm1, x);
56 1619f52c 2010-03-12 rsc 			mpfree(r);
57 1619f52c 2010-03-12 rsc 			if(mpcmp(x, mpone) > 0)
58 1619f52c 2010-03-12 rsc 				break;
59 1619f52c 2010-03-12 rsc 		}
60 0fc65b37 2004-03-21 devnull 
61 cbeb0b26 2006-04-01 devnull 		/* y = x**q mod n */
62 0fc65b37 2004-03-21 devnull 		mpexp(x, q, n, y);
63 0fc65b37 2004-03-21 devnull 
64 0fc65b37 2004-03-21 devnull 		if(mpcmp(y, mpone) == 0 || mpcmp(y, nm1) == 0)
65 1619f52c 2010-03-12 rsc 			continue;
66 0fc65b37 2004-03-21 devnull 
67 1619f52c 2010-03-12 rsc 		for(j = 1;; j++){
68 1619f52c 2010-03-12 rsc 			if(j >= k) {
69 1619f52c 2010-03-12 rsc 				isprime = 0;
70 1619f52c 2010-03-12 rsc 				goto done;
71 1619f52c 2010-03-12 rsc 			}
72 0fc65b37 2004-03-21 devnull 			mpmul(y, y, x);
73 cbeb0b26 2006-04-01 devnull 			mpmod(x, n, y);	/* y = y*y mod n */
74 0fc65b37 2004-03-21 devnull 			if(mpcmp(y, nm1) == 0)
75 1619f52c 2010-03-12 rsc 				break;
76 0fc65b37 2004-03-21 devnull 			if(mpcmp(y, mpone) == 0){
77 0fc65b37 2004-03-21 devnull 				isprime = 0;
78 0fc65b37 2004-03-21 devnull 				goto done;
79 0fc65b37 2004-03-21 devnull 			}
80 0fc65b37 2004-03-21 devnull 		}
81 0fc65b37 2004-03-21 devnull 	}
82 1619f52c 2010-03-12 rsc 	isprime = 1;
83 0fc65b37 2004-03-21 devnull done:
84 0fc65b37 2004-03-21 devnull 	mpfree(y);
85 0fc65b37 2004-03-21 devnull 	mpfree(x);
86 0fc65b37 2004-03-21 devnull 	mpfree(q);
87 0fc65b37 2004-03-21 devnull 	mpfree(nm1);
88 0fc65b37 2004-03-21 devnull 	return isprime;
89 0fc65b37 2004-03-21 devnull }