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 */
10 probably_prime(mpint *n, int nrep)
12 int j, k, rep, nbits, isprime;
13 mpint *nm1, *q, *x, *y, *r;
16 sysfatal("negative prime candidate");
22 if(k == 2) /* 2 is prime */
24 if(k < 2) /* 1 is not prime */
26 if((n->p[0] & 1) == 0) /* even is not prime */
29 /* test against small prime numbers */
30 if(smallprimetest(n) < 0)
33 /* fermat test, 2^n mod n == 2 if p is prime */
46 mpsub(n, mpone, nm1); /* nm1 = n - 1 */
49 mpright(nm1, k, q); /* q = (n-1)/2**k */
51 for(rep = 0; rep < nrep; rep++){
53 /* find x = random in [2, n-2] */
54 r = mprand(nbits, prng, nil);
57 if(mpcmp(x, mpone) > 0)
64 if(mpcmp(y, mpone) == 0 || mpcmp(y, nm1) == 0)
73 mpmod(x, n, y); /* y = y*y mod n */
74 if(mpcmp(y, nm1) == 0)
76 if(mpcmp(y, mpone) == 0){