tprobably_prime.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tprobably_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 }