math/strongprime.c: Choose the smaller primes' sizes more carefully.

The old code would indeed, as the warning in the comment said, produce
numbers which are larger than requested because the component primes'
sizes were chosen in a naïve manner.  I've now (eventually!) thought
about the issue some more and come up with a better approach.

The `BITSLOP' macro is now gone, replaced by a carefully chosen value
supported by some actual mathematics.  As a result, the warning comments
have been removed.  Also, `strongprime' will fail if it actually returns
a number of the wrong size.

diff --git a/math/strongprime.c b/math/strongprime.c
index a82bfad090b44bc0bc7c660f5fdbb43e13634927..fc20bfea8077714b391b37dad61a74a64e8fa518 100644 (file)
--- 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);

diff --git a/math/strongprime.h b/math/strongprime.h
index 6f39f500030ce073883c674ca38c3ad1214ab2d8..90102362b4337a8f40131f5781fe185ee59417f1 100644 (file)
--- a/math/strongprime.h
+++ b/math/strongprime.h
@@ -85,9 +85,6 @@ extern mp *strongprime_setup(const char */*name*/, mp */*d*/, pfilt */*f*/,
  *               * %$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.
  */
 
 extern mp *strongprime(const char */*name*/, mp */*d*/, unsigned /*nbits*/,