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