X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~mdw/git/catacomb/blobdiff_plain/285bf989997b8dc94a0783e260fe73787c7ae767..32bd36cff950788ec62bd36a6d437da12aa4fd0f:/math/strongprime.c diff --git a/math/strongprime.c b/math/strongprime.c index a82bfad0..fc20bfea 100644 --- a/math/strongprime.c +++ b/math/strongprime.c @@ -63,27 +63,43 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits, { mp *s, *t, *q; dstr dn = DSTR_INIT; - size_t nb; + unsigned slop, nb, u, i; mp *rr = d; pgen_filterctx c; pgen_jumpctx j; rabin rb; - /* --- The bitslop parameter --- * + /* --- Figure out how large the smaller primes should be --- * * - * There's quite a lot of prime searching to be done. The constant - * @BITSLOP@ is a (low) approximation to the base-2 log of the expected - * number of steps to find a prime number. Experimentation shows that - * numbers around 10 seem to be good. + * We want them to be `as large as possible', subject to the constraint + * that we produce a number of the requested size at the end. This is + * tricky, because the final prime search is going to involve quite large + * jumps from its starting point; the size of the jumps are basically + * determined by our choice here, and if they're too big then we won't find + * a prime in time. + * + * Let's suppose we're trying to make an %$N$%-bit prime. The expected + * number of steps tends to increase linearly with size, i.e., we need to + * take about %2^k N$% steps for some %$k$%. If we're jumping by a + * %$J$%-bit quantity each time, from an %$N$%-bit starting point, then we + * will only be able to find a match if %$2^k N 2^{J-1} \le 2^{N-1}$%, + * i.e., if %$J \le N - (k + \log_2 N)$%. + * + * Experimentation shows that taking %$k + \log_2 N = 12$% works well for + * %$N = 1024$%, so %$k = 2$%. */ -#define BITSLOP 12 + for (i = 1; i && nbits >> i; i <<= 1); assert(i); + for (slop = 2, nb = nbits; nb > 1; i >>= 1) { + u = nb >> i; + if (u) { slop += i; nb = u; } + } + if (nbits/2 <= slop) return (0); /* --- Choose two primes %$s$% and %$t$% of half the required size --- */ - if (nbits/2 <= BITSLOP) return (0); - nb = nbits/2 - BITSLOP; + nb = nbits/2 - slop; c.step = 1; rr = mprand(rr, nb, r, 1); @@ -102,12 +118,12 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits, rr = mp_lsl(rr, t, 1); pfilt_create(&c.f, rr); - rr = mp_lsl(rr, rr, BITSLOP - 1); + rr = mp_lsl(rr, rr, slop - 1); rr = mp_add(rr, rr, MP_ONE); DRESET(&dn); dstr_putf(&dn, "%s [r]", name); j.j = &c.f; q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &j, - rabin_iters(nb + BITSLOP), pgen_test, &rb); + rabin_iters(nb + slop), pgen_test, &rb); pfilt_destroy(&c.f); if (!q) goto fail_r; @@ -159,8 +175,6 @@ fail_s: mp_drop(rr); dstr_destroy(&dn); return (0); - -#undef BITSLOP } /* --- @strongprime@ --- * @@ -180,9 +194,6 @@ fail_s: * * %$p - 1$% has a large prime factor %$r$%, * * %$p + 1$% has a large prime factor %$s$%, and * * %$r - 1$% has a large prime factor %$t$%. - * - * The numbers produced may be slightly larger than requested, - * by a few bits. */ mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r, @@ -199,6 +210,7 @@ mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r, j.j = &f; p = pgen(name, p, p, event, ectx, n, pgen_jump, &j, rabin_iters(nbits), pgen_test, &rb); + if (mp_bits(p) != nbits) { mp_drop(p); return (0); } pfilt_destroy(&f); mp_drop(d); return (p);