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
probably_prime(mpint * n,int nrep)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) /* 2 is prime */
25 return 1;
26 if(k < 2) /* 1 is not prime */
27 return 0;
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 }
92