math/, pub/: Generate primes with exactly the right size.

[catacomb] / math / strongprime.c

diff --git a/math/strongprime.c b/math/strongprime.c
index a2dd43858cb1d6d1ef1c624a76d97e425f335c5f..60a7cc75310b9c926aadb4afc6c3a206a04689bb 100644 (file)
--- a/math/strongprime.c
+++ b/math/strongprime.c
@@ -63,6 +63,7 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
 {
   mp *s, *t, *q;
   dstr dn = DSTR_INIT;
+  size_t nb;
 
   mp *rr = d;
   pgen_filterctx c;
@@ -82,19 +83,19 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
   /* --- 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 - BITSLOP;
   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$% --- */
@@ -105,9 +106,8 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
   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 + BITSLOP), pgen_test, &rb);
   pfilt_destroy(&c.f);
   if (!q)
     goto fail_r;
@@ -132,13 +132,13 @@ mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
   /* --- Now find %$p = p_0 + 2jrs$% for some %$j$% --- */
 
   {
-    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);
+    y = mp_lsl(MP_NEW, MP_ONE, nbits - 1);
+    rr = mp_leastcongruent(rr, y, rr, x);
+    mp_drop(x); mp_drop(y);
   }
 
   /* --- Return the result --- */