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