#include #include #include #include int miller(mpint *n, int k) { mpint *d, *a, *t, *x; int r, i, rc; d = mpnew(0); t = mpnew(0); a = mpnew(0); x = mpnew(0); mpsub(n, mpone, t); r = mplowbits0(t); mpright(t, r, d); if(r == 0) return 0; rc = 1; while(k-- > 0){ itomp(3, t); mpsub(n, t, t); mpnrand(t, genrandom, a); mpadd(a, mptwo, a); mpsub(n, mpone, t); mpexp(a, d, n, x); if(mpcmp(x, mpone) == 0 || mpcmp(x, t) == 0) continue; for(i = 0; i < r - 1; i++){ mpexp(x, mptwo, n, x); if(mpcmp(x, mpone) == 0) break; if(mpcmp(x, t) == 0) goto next; } rc = 0; break; next: ; } out: mpfree(d); mpfree(t); mpfree(a); mpfree(x); return rc; } void main() { mpint *n; fmtinstall('B', mpfmt); n = mpnew(0); while(!miller(mprand(64, genrandom, n), 10000)) ; print("%.10B\n", n); }