From 32bd36cff950788ec62bd36a6d437da12aa4fd0f Mon Sep 17 00:00:00 2001
From: Mark Wooding <mdw@distorted.org.uk>
Date: Thu, 26 May 2016 09:26:09 +0100
Subject: [PATCH] math/strongprime.c: Choose the smaller primes' sizes more
 carefully.
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit

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.
---
 math/strongprime.c | 44 ++++++++++++++++++++++++++++----------------
 math/strongprime.h |  3 ---
 2 files changed, 28 insertions(+), 19 deletions(-)

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);
diff --git a/math/strongprime.h b/math/strongprime.h
index 6f39f500..90102362 100644
--- 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*/,
-- 
2.11.0