progs/perftest.c: Use from Glibc syscall numbers.

[catacomb] / math / strongprime.c

diff --git a/math/strongprime.c b/math/strongprime.c
index fc20bfe..4ea6253 100644 (file)
--- a/math/strongprime.c
+++ b/math/strongprime.c
@@ -28,6 +28,7 @@
 /*----- Header files ------------------------------------------------------*/
 
 #include <mLib/dstr.h>
+#include <mLib/macros.h>
 
 #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,
@@ -87,11 +98,11 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
    * 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$%.
+   * %$N = 1024$%, so %$k = 2$%.  Add a few extra bits for luck.
    */
 
   for (i = 1; i && nbits >> i; i <<= 1); assert(i);
-  for (slop = 2, nb = nbits; nb > 1; i >>= 1) {
+  for (slop = 6, nb = nbits; nb > 1; i >>= 1) {
     u = nb >> i;
     if (u) { slop += i; nb = u; }
   }
@@ -114,7 +125,10 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
                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);
@@ -128,31 +142,34 @@ 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^{-1} \bmod r) - 1$%.  Then %$p_0 + 1$% is
+   * clearly a multiple of %$s$%, and
    *
-   * This computes %$p_0 = 2 s (s^{r - 2} \bmod r) - 1$%.
+   *   %$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, *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);
   }