chiark / gitweb /
math/, pub/: Generate primes with exactly the right size.
authorMark Wooding <mdw@distorted.org.uk>
Thu, 26 May 2016 08:26:09 +0000 (09:26 +0100)
committerMark Wooding <mdw@distorted.org.uk>
Fri, 24 Jun 2016 00:17:49 +0000 (01:17 +0100)
Previously, the `strongprime' and `limlee' machinery had a tendency to
generate primes which were a few bits shorter than actually requested.
Fix this unfortunate state of affairs by using a more careful analysis
of sizes of things.

The Lim--Lee prime generation has been changed quite a bit.  The large
factor now might not be suitable, so there's some new machinery for
tracking how useful it's being for generating numbers of the right size
and choosing a different one if it gets out of whack.  Unfortunately,
this means applying a rather unpleasant hack to the structure layout to
preserve ABI compatibility.

math/limlee.c
math/limlee.h
math/strongprime.c
pub/rsa-gen.c

index b20965a90550d69a7d8c2fb89b93ee5856a550f6..13bbc2128465f05346e4f8a1cc54e43fd1d6c69e 100644 (file)
@@ -117,7 +117,7 @@ static void llgen(limlee_factor *f, unsigned pl, limlee_stepctx *l)
 again:
   p = mprand(l->newp, pl, l->r, 1);
   pf.step = 2;
-  p = pgen(l->d.buf, p, p, l->iev, l->iec, 0, pgen_filter, &pf,
+  p = pgen(l->u.s.name, p, p, l->iev, l->iec, 0, pgen_filter, &pf,
           rabin_iters(pl), pgen_test, &r);
   if (!p)
     goto again;
@@ -146,18 +146,12 @@ static const limlee_primeops primeops_simple = { llgen, llfree };
 static int init(pgen_event *ev, limlee_stepctx *l)
 {
   size_t i;
-  unsigned qql;
 
   /* --- First of all, decide on a number of factors to make --- */
 
   l->nf = l->pl / l->ql;
-  qql = l->pl % l->ql;
-  if (!l->nf)
-    return (PGEN_ABORT);
-  else if (qql && l->nf > 1) {
-    l->nf--;
-    qql += l->ql;
-  }
+  if (l->nf < 2) return (PGEN_ABORT);
+  l->nf--;
 
   /* --- Now decide on how many primes I'll actually generate --- *
    *
@@ -180,20 +174,14 @@ static int init(pgen_event *ev, limlee_stepctx *l)
   /* --- Other bits of initialization --- */
 
   l->seq = 0;
-  dstr_create(&l->d);
   if (!l->pops) {
     l->pops = &primeops_simple;
     l->pc = 0;
   }
 
-  /* --- Find a big prime --- */
+  /* --- Find a big prime later --- */
 
-  if (!qql)
-    l->qq.p = 0;
-  else {
-    dstr_putf(&l->d, "%s*", ev->name);
-    l->pops->pgen(&l->qq, qql, l);
-  }
+  l->qq.p = 0;
 
   return (PGEN_TRY);
 }
@@ -211,8 +199,11 @@ static int init(pgen_event *ev, limlee_stepctx *l)
 
 static int next(int rq, pgen_event *ev, limlee_stepctx *l)
 {
+  dstr d = DSTR_INIT;
   mp *p;
   int rc;
+  int dist;
+  unsigned nb;
 
   mp_drop(ev->m);
 
@@ -230,32 +221,64 @@ static int next(int rq, pgen_event *ev, limlee_stepctx *l)
     }
     rq = PGEN_TRY; /* For next time through */
 
+    /* --- If the large factor is performing badly, make a new one --- */
+
+    if (l->qq.p) {
+      dist = l->u.s.disp < 0 ? -l->u.s.disp : l->u.s.disp;
+      if (dist && dist > l->u.s.steps/dist) {
+       l->pops->pfree(&l->qq, l);
+       l->qq.p = 0;
+      }
+    }
+
     /* --- Gather up some factors --- */
 
-    if (l->qq.p)
-      mpmul_add(&mm, l->qq.p);
+    if (l->qq.p) mpmul_add(&mm, l->qq.p);
     for (i = 0; i < l->poolsz; i++) {
       if (!l->c[i])
        continue;
       if (!l->v[i].p) {
-       DRESET(&l->d);
-       dstr_putf(&l->d, "%s_%lu", ev->name, l->seq++);
+       DRESET(&d);
+       dstr_putf(&d, "%s_%lu", ev->name, l->seq++);
+       l->u.s.name = d.buf;
        l->pops->pgen(&l->v[i], l->ql, l);
       }
       mpmul_add(&mm, l->v[i].p);
     }
 
-    /* --- Check it for small factors --- */
+    /* --- Check on the large factor --- */
 
     p = mpmul_done(&mm);
+    if (!l->qq.p) {
+      DRESET(&d);
+      dstr_putf(&d, "%s*_%lu", ev->name, l->seq++);
+      l->u.s.name = d.buf;
+      l->pops->pgen(&l->qq, l->pl - mp_bits(p), l);
+      l->u.s.steps = l->u.s.disp = 0;
+      p = mp_mul(p, p, l->qq.p);
+    }
     p = mp_lsl(p, p, 1);
     p->v[0] |= 1;
+
+    nb = mp_bits(p);
+    l->u.s.steps++;
+    if (nb < l->pl) {
+      l->u.s.disp--;
+      continue;
+    } else if (nb > l->pl) {
+      l->u.s.disp++;
+      continue;
+    }
+
+    /* --- Check it for small factors --- */
+
     if ((rc = pfilt_smallfactor(p)) != PGEN_FAIL)
       break;
     mp_drop(p);
   }
 
   ev->m = p;
+  DDESTROY(&d);
   return (rc);
 }
 
@@ -305,7 +328,6 @@ static int done(pgen_event *ev, limlee_stepctx *l)
   /* --- Free other resources --- */
 
   xfree(l->c);
-  dstr_destroy(&l->d);
 
   /* --- Done --- */
 
index d22b5a7ec2b6215243c0f78c2b66d0d5f8e01f42..b46ae5c9c3d7806cc568ad6d6ade8b77cd30721c 100644 (file)
@@ -77,7 +77,13 @@ typedef struct limlee_stepctx {
   octet *c;                            /* Combination byte-flag vector */
   unsigned long seq;                   /* Sequence number for primes */
   size_t poolsz;                       /* Size of the small-prime pool */
-  dstr d;                              /* String for subprime name */
+  union {
+    dstr d;                            /* Obsolete; for ABI compat */
+    struct {
+      char *name;                      /* Name, for @primeops@ */
+      int steps, disp;                 /* Track how good @qq@ is */
+    } s;
+  } u;
   limlee_factor qq;                    /* Big prime to pick up slack */
 
 } limlee_stepctx;
index a2dd43858cb1d6d1ef1c624a76d97e425f335c5f..60a7cc75310b9c926aadb4afc6c3a206a04689bb 100644 (file)
@@ -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 --- */
index 730673e33388f0be37377dd73c07d2a1ca554f27..0653d1c3b2dafaf1bd2ca1b9eaadf0282e518790 100644 (file)
@@ -77,15 +77,25 @@ again:
   if ((rp->p = strongprime("p", MP_NEWSEC, nbits/2, r, n, event, ectx)) == 0)
     goto fail_p;
 
-  /* --- Do painful fiddling with GCD steppers --- */
+  /* --- Do painful fiddling with GCD steppers --- *
+   *
+   * Also, arrange that %$q \ge \lceil 2^{N-1}/p \rceil$%, so that %$p q$%
+   * has the right length.
+   */
 
   {
     mp *q;
+    mp *t = MP_NEW, *u = MP_NEW;
     rabin rb;
 
     if ((q = strongprime_setup("q", MP_NEWSEC, &g.jp, nbits / 2,
                               r, n, event, ectx)) == 0)
       goto fail_q;
+    t = mp_lsl(t, MP_ONE, nbits - 1);
+    mp_div(&t, &u, t, rp->p);
+    if (!MP_ZEROP(u)) t = mp_add(t, t, MP_ONE);
+    if (MP_CMP(q, <, t)) q = mp_leastcongruent(q, t, q, g.jp.m);
+
     g.r = mp_lsr(MP_NEW, rp->p, 1);
     g.g = MP_NEW;
     g.max = MP_256;