probably_prime.c (1662B)
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; 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 for(;;){ 53 /* find 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 break; 59 } 60 61 /* y = x**q mod n */ 62 mpexp(x, q, n, y); 63 64 if(mpcmp(y, mpone) == 0 || mpcmp(y, nm1) == 0) 65 continue; 66 67 for(j = 1;; j++){ 68 if(j >= k) { 69 isprime = 0; 70 goto done; 71 } 72 mpmul(y, y, x); 73 mpmod(x, n, y); /* y = y*y mod n */ 74 if(mpcmp(y, nm1) == 0) 75 break; 76 if(mpcmp(y, mpone) == 0){ 77 isprime = 0; 78 goto done; 79 } 80 } 81 } 82 isprime = 1; 83 done: 84 mpfree(y); 85 mpfree(x); 86 mpfree(q); 87 mpfree(nm1); 88 return isprime; 89 }