123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384 |
- #include "os.h"
- #include <mp.h>
- #include <libsec.h>
- // Miller-Rabin probabilistic primality testing
- // Knuth (1981) Seminumerical Algorithms, p.379
- // Menezes et al () Handbook, p.39
- // 0 if composite; 1 if almost surely prime, Pr(err)<1/4**nrep
- int
- probably_prime(mpint *n, int nrep)
- {
- int j, k, rep, nbits, isprime = 1;
- mpint *nm1, *q, *x, *y, *r;
- if(n->sign < 0)
- sysfatal("negative prime candidate");
- if(nrep <= 0)
- nrep = 18;
- k = mptoi(n);
- if(k == 2) // 2 is prime
- return 1;
- if(k < 2) // 1 is not prime
- return 0;
- if((n->p[0] & 1) == 0) // even is not prime
- return 0;
- // test against small prime numbers
- if(smallprimetest(n) < 0)
- return 0;
- // fermat test, 2^n mod n == 2 if p is prime
- x = uitomp(2, nil);
- y = mpnew(0);
- mpexp(x, n, n, y);
- k = mptoi(y);
- if(k != 2){
- mpfree(x);
- mpfree(y);
- return 0;
- }
- nbits = mpsignif(n);
- nm1 = mpnew(nbits);
- mpsub(n, mpone, nm1); // nm1 = n - 1 */
- k = mplowbits0(nm1);
- q = mpnew(0);
- mpright(nm1, k, q); // q = (n-1)/2**k
- for(rep = 0; rep < nrep; rep++){
-
- // x = random in [2, n-2]
- r = mprand(nbits, prng, nil);
- mpmod(r, nm1, x);
- mpfree(r);
- if(mpcmp(x, mpone) <= 0)
- continue;
- // y = x**q mod n
- mpexp(x, q, n, y);
- if(mpcmp(y, mpone) == 0 || mpcmp(y, nm1) == 0)
- goto done;
- for(j = 1; j < k; j++){
- mpmul(y, y, x);
- mpmod(x, n, y); // y = y*y mod n
- if(mpcmp(y, nm1) == 0)
- goto done;
- if(mpcmp(y, mpone) == 0){
- isprime = 0;
- goto done;
- }
- }
- isprime = 0;
- }
- done:
- mpfree(y);
- mpfree(x);
- mpfree(q);
- mpfree(nm1);
- return isprime;
- }
|