90 lines
1.6 KiB
C
90 lines
1.6 KiB
C
#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;
|
|
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++){
|
|
for(;;){
|
|
/* find x = random in [2, n-2] */
|
|
r = mprand(nbits, prng, nil);
|
|
mpmod(r, nm1, x);
|
|
mpfree(r);
|
|
if(mpcmp(x, mpone) > 0)
|
|
break;
|
|
}
|
|
|
|
/* y = x**q mod n */
|
|
mpexp(x, q, n, y);
|
|
|
|
if(mpcmp(y, mpone) == 0 || mpcmp(y, nm1) == 0)
|
|
continue;
|
|
|
|
for(j = 1;; j++){
|
|
if(j >= k) {
|
|
isprime = 0;
|
|
goto done;
|
|
}
|
|
mpmul(y, y, x);
|
|
mpmod(x, n, y); /* y = y*y mod n */
|
|
if(mpcmp(y, nm1) == 0)
|
|
break;
|
|
if(mpcmp(y, mpone) == 0){
|
|
isprime = 0;
|
|
goto done;
|
|
}
|
|
}
|
|
}
|
|
isprime = 1;
|
|
done:
|
|
mpfree(y);
|
|
mpfree(x);
|
|
mpfree(q);
|
|
mpfree(nm1);
|
|
return isprime;
|
|
}
|