Warning: Attempt to read property "date" on null in /usr/local/www/websvn.planix.org/blame.php on line 247

Warning: Attempt to read property "msg" on null in /usr/local/www/websvn.planix.org/blame.php on line 247
WebSVN – planix.SVN – Blame – /os/branches/feature_posix/sys/src/cmd/unix/drawterm/libsec/probably_prime.c – Rev 2

Subversion Repositories planix.SVN

Rev

Go to most recent revision | Details | 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 
// 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 
}