X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/6687eff59d78369e4d19fb0b3a9a0bb6035d1031..HEAD:/math/strongprime.c diff --git a/math/strongprime.c b/math/strongprime.c index a2dd4385..4ea62537 100644 --- a/math/strongprime.c +++ b/math/strongprime.c @@ -28,6 +28,7 @@ /*----- Header files ------------------------------------------------------*/ #include +#include #include "grand.h" #include "mp.h" @@ -39,6 +40,9 @@ /*----- Main code ---------------------------------------------------------*/ +/* Oh, just shut up. */ +CLANG_WARNING("-Wempty-body") + /* --- @strongprime_setup@ --- * * * Arguments: @const char *name@ = pointer to name root @@ -55,7 +59,14 @@ * Use: Sets up for a strong prime search, so that primes with * particular properties can be found. It's probably important * to note that the number left in the filter context @f@ is - * congruent to 2 (mod 4). + * congruent to 2 (mod 4); that the jump value is twice the + * product of two large primes; and that the starting point is + * at least %$3 \cdot 2^{N-2}$%. (Hence, if you multiply two + * such numbers, the product is at least + * + * %$9 \cdot 2^{2N-4} > 2^{2N-1}$% + * + * i.e., it will be (at least) a %$2 N$%-bit value. */ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits, @@ -63,82 +74,104 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits, { mp *s, *t, *q; dstr dn = DSTR_INIT; + 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 --- * + * + * 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)$%. * - * 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. + * Experimentation shows that taking %$k + \log_2 N = 12$% works well for + * %$N = 1024$%, so %$k = 2$%. Add a few extra bits for luck. */ -#define BITSLOP 12 + for (i = 1; i && nbits >> i; i <<= 1); assert(i); + for (slop = 6, 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 --- */ - assert(((void)"nbits too small in strongprime_setup", nbits/2 > BITSLOP)); - nbits = nbits/2 - BITSLOP; + nb = nbits/2 - slop; c.step = 1; - rr = mprand(rr, nbits, r, 1); + rr = mprand(rr, nb, r, 1); DRESET(&dn); dstr_putf(&dn, "%s [s]", name); if ((s = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c, - rabin_iters(nbits), pgen_test, &rb)) == 0) + rabin_iters(nb), pgen_test, &rb)) == 0) goto fail_s; - rr = mprand(rr, nbits, r, 1); + rr = mprand(rr, nb, r, 1); DRESET(&dn); dstr_putf(&dn, "%s [t]", name); if ((t = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c, - rabin_iters(nbits), pgen_test, &rb)) == 0) + rabin_iters(nb), pgen_test, &rb)) == 0) goto fail_t; - /* --- Choose a suitable value for %$r = 2it + 1$% for some %$i$% --- */ + /* --- Choose a suitable value for %$r = 2it + 1$% for some %$i$% --- * + * + * Then %$r \equiv 1 \pmod{t}$%, i.e., %$r - 1$% is a multiple of %$t$%. + */ 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; - nbits += BITSLOP; q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &j, - rabin_iters(nbits), pgen_test, &rb); + rabin_iters(nb + slop), pgen_test, &rb); pfilt_destroy(&c.f); if (!q) goto fail_r; - /* --- Select a suitable starting-point for finding %$p$% --- * + /* --- Select a suitable congruence class for %$p$% --- * * - * This computes %$p_0 = 2 s (s^{r - 2} \bmod r) - 1$%. + * This computes %$p_0 = 2 s (s^{-1} \bmod r) - 1$%. Then %$p_0 + 1$% is + * clearly a multiple of %$s$%, and + * + * %$p_0 - 1 \equiv 2 s s^{-1} - 2 \equiv 0 \pmod{r}$% + * + * is a multiple of %$r$%. */ - { - mpmont mm; - - mpmont_create(&mm, q); - rr = mp_sub(rr, q, MP_TWO); - rr = mpmont_exp(&mm, rr, s, rr); - mpmont_destroy(&mm); - rr = mp_mul(rr, rr, s); - rr = mp_lsl(rr, rr, 1); - rr = mp_sub(rr, rr, MP_ONE); - } + rr = mp_modinv(rr, s, q); + rr = mp_mul(rr, rr, s); + rr = mp_lsl(rr, rr, 1); + rr = mp_sub(rr, rr, MP_ONE); - /* --- Now find %$p = p_0 + 2jrs$% for some %$j$% --- */ + /* --- Pick a starting point for the search --- * + * + * Select %$3 \cdot 2^{N-2} < p_1 < 2^N$% at random, only with + * %$p_1 \equiv p_0 \pmod{2 r s}$. + */ { - mp *x; + mp *x, *y; x = mp_mul(MP_NEW, q, s); x = mp_lsl(x, x, 1); - pfilt_create(f, x); - x = mp_lsl(x, x, BITSLOP - 1); - rr = mp_add(rr, rr, x); - mp_drop(x); + pfilt_create(f, x); /* %$2 r s$% */ + y = mprand(MP_NEW, nbits, r, 0); + y = mp_setbit(y, y, nbits - 2); + rr = mp_leastcongruent(rr, y, rr, x); + mp_drop(x); mp_drop(y); } /* --- Return the result --- */ @@ -159,8 +192,6 @@ fail_s: mp_drop(rr); dstr_destroy(&dn); return (0); - -#undef BITSLOP } /* --- @strongprime@ --- * @@ -180,24 +211,26 @@ 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, unsigned n, pgen_proc *event, void *ectx) { + mp *p; pfilt f; pgen_jumpctx j; rabin rb; - d = strongprime_setup(name, d, &f, nbits, r, n, event, ectx); + if (d) mp_copy(d); + p = strongprime_setup(name, d, &f, nbits, r, n, event, ectx); + if (!p) { mp_drop(d); return (0); } j.j = &f; - d = pgen(name, d, d, event, ectx, n, pgen_jump, &j, + 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); - return (d); + mp_drop(d); + return (p); } /*----- That's all, folks -------------------------------------------------*/