Subversion Repositories planix.SVN

Rev

Rev 33 | Details | Compare with Previous | Last modification | View Log | RSS feed

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