math/strongprime.c: Clamp the starting point.

Now the result will be in the upper quarter of the `obvious' range, and
the product of two such values is guaranteed to have the desired number
of bits.  This saves callers from doing stupid things like trying to
clamp one of the factors by hand, which ends up significantly biasing
the second factor.  (This isn't very bad, because there's a /lot/ of
randomness in the chosen congruence class, but it's good to fix this
sort of thing.)

diff --git a/math/strongprime.c b/math/strongprime.c
index e883492..9eab8b2 100644 (file)
--- a/math/strongprime.c
+++ b/math/strongprime.c
@@ -55,7 +55,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,
@@ -128,7 +135,7 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
   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$%.
    */
@@ -145,14 +152,19 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
     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, *y;
     x = mp_mul(MP_NEW, q, s);
     x = mp_lsl(x, x, 1);
-    pfilt_create(f, x);
-    y = mp_lsl(MP_NEW, MP_ONE, nbits - 1);
+    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);
   }

diff --git a/math/strongprime.h b/math/strongprime.h
index 9010236..14ad39e 100644 (file)
--- a/math/strongprime.h
+++ b/math/strongprime.h
@@ -60,7 +60,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.
  */
 
 extern mp *strongprime_setup(const char */*name*/, mp */*d*/, pfilt */*f*/,