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) /* 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