#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 */intprobably_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;}
#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;
}