sys/src/libsec/port/probably_prime.c

   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)               /* 1 is not prime */
  25                 return 0;
  26         if(k == 2 || k == 3)    /* 2, 3 is prime */
  27                 return 1;
  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 }