Rearrange the file tree.
[u/mdw/catacomb] / math / limlee.c
diff --git a/math/limlee.c b/math/limlee.c
new file mode 100644 (file)
index 0000000..b20965a
--- /dev/null
@@ -0,0 +1,396 @@
+/* -*-c-*-
+ *
+ * Generate Lim-Lee primes
+ *
+ * (c) 2000 Straylight/Edgeware
+ */
+
+/*----- Licensing notice --------------------------------------------------*
+ *
+ * This file is part of Catacomb.
+ *
+ * Catacomb is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Library General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * Catacomb is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU Library General Public
+ * License along with Catacomb; if not, write to the Free
+ * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+ * MA 02111-1307, USA.
+ */
+
+/*----- Header files ------------------------------------------------------*/
+
+#include <mLib/alloc.h>
+#include <mLib/dstr.h>
+
+#include "limlee.h"
+#include "mpmul.h"
+#include "mprand.h"
+#include "pgen.h"
+#include "rabin.h"
+
+/*----- Stepping through combinations -------------------------------------*/
+
+/* --- @comb_init@ --- *
+ *
+ * Arguments:  @octet *c@ = pointer to byte-flag array
+ *             @unsigned n@ = number of items in the array
+ *             @unsigned r@ = number of desired items
+ *
+ * Returns:    ---
+ *
+ * Use:                Initializes a byte-flag array which, under the control of
+ *             @comb_next@, will step through all combinations of @r@ chosen
+ *             elements.
+ */
+
+static void comb_init(octet *c, unsigned n, unsigned r)
+{
+  memset(c, 0, n - r);
+  memset(c + (n - r), 1, r);
+}
+
+/* --- @comb_next@ --- *
+ *
+ * Arguments:  @octet *c@ = pointer to byte-flag array
+ *             @unsigned n@ = number of items in the array
+ *             @unsigned r@ = number of desired items
+ *
+ * Returns:    Nonzero if another combination was returned, zero if we've
+ *             reached the end.
+ *
+ * Use:                Steps on to the next combination in sequence.
+ */
+
+static int comb_next(octet *c, unsigned n, unsigned r)
+{
+  unsigned g = 0;
+
+  /* --- How the algorithm works --- *
+   *
+   * Set bits start at the end and work their way towards the start.
+   * Excepting bits already at the start, we scan for the lowest set bit, and
+   * move it one place nearer the start.  A group of bits at the start are
+   * counted and reset just below the `moved' bit.  If there is no moved bit
+   * then we're done.
+   */
+
+  /* --- Count the group at the start --- */
+
+  for (; *c; c++) {
+    g++;
+    *c = 0;
+  }
+  if (g == r)
+    return (0);
+
+  /* --- Move the next bit down one --- *
+   *
+   * There must be one, because otherwise we'd have counted %$r$% bits
+   * earlier.
+   */
+
+  for (; !*c; c++)
+    ;
+  *c = 0;
+  g++;
+  for (; g; g--)
+    *--c = 1;
+  return (1);
+}
+
+/*----- Default prime generator -------------------------------------------*/
+
+static void llgen(limlee_factor *f, unsigned pl, limlee_stepctx *l)
+{
+  pgen_filterctx pf;
+  rabin r;
+  mp *p;
+
+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,
+          rabin_iters(pl), pgen_test, &r);
+  if (!p)
+    goto again;
+  f->p = p;
+}
+
+static void llfree(limlee_factor *f, limlee_stepctx *l)
+{
+  mp_drop(f->p);
+}
+
+static const limlee_primeops primeops_simple = { llgen, llfree };
+
+/*----- Lim-Lee stepper ---------------------------------------------------*/
+
+/* --- @init@ --- *
+ *
+ * Arguments:  @pgen_event *ev@ = pointer to event block
+ *             @limlee_stepctx *l@ = pointer to Lim-Lee context
+ *
+ * Returns:    A @PGEN@ result code.
+ *
+ * Use:                Initializes the stepper.
+ */
+
+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;
+  }
+
+  /* --- Now decide on how many primes I'll actually generate --- *
+   *
+   * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
+   * library.
+   */
+
+  l->poolsz = l->nf * 3 + 5;
+  if (l->poolsz < 25)
+    l->poolsz = 25;
+
+  /* --- Allocate and initialize the various tables --- */
+
+  l->c = xmalloc(l->poolsz);
+  l->v = xmalloc(l->poolsz * sizeof(limlee_factor));
+  comb_init(l->c, l->poolsz, l->nf);
+  for (i = 0; i < l->poolsz; i++)
+    l->v[i].p = 0;
+
+  /* --- 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 --- */
+
+  if (!qql)
+    l->qq.p = 0;
+  else {
+    dstr_putf(&l->d, "%s*", ev->name);
+    l->pops->pgen(&l->qq, qql, l);
+  }
+
+  return (PGEN_TRY);
+}
+
+/* --- @next@ --- *
+ *
+ * Arguments:  @int rq@ = request which triggered this call
+ *             @pgen_event *ev@ = pointer to event block
+ *             @limlee_stepctx *l@ = pointer to Lim-Lee context
+ *
+ * Returns:    A @PGEN@ result code.
+ *
+ * Use:                Initializes the stepper.
+ */
+
+static int next(int rq, pgen_event *ev, limlee_stepctx *l)
+{
+  mp *p;
+  int rc;
+
+  mp_drop(ev->m);
+
+  for (;;) {
+    size_t i;
+    mpmul mm = MPMUL_INIT;
+
+    /* --- Step on to next combination --- */
+
+    if (rq == PGEN_TRY && !comb_next(l->c, l->poolsz, l->nf)) {
+      for (i = 0; i < l->poolsz; i++) {
+       l->pops->pfree(&l->v[i], l);
+       l->v[i].p = 0;
+      }
+    }
+    rq = PGEN_TRY; /* For next time through */
+
+    /* --- Gather up some factors --- */
+
+    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++);
+       l->pops->pgen(&l->v[i], l->ql, l);
+      }
+      mpmul_add(&mm, l->v[i].p);
+    }
+
+    /* --- Check it for small factors --- */
+
+    p = mpmul_done(&mm);
+    p = mp_lsl(p, p, 1);
+    p->v[0] |= 1;
+    if ((rc = pfilt_smallfactor(p)) != PGEN_FAIL)
+      break;
+    mp_drop(p);
+  }
+
+  ev->m = p;
+  return (rc);
+}
+
+/* --- @done@ --- *
+ *
+ * Arguments:  @pgen_event *ev@ = pointer to event block
+ *             @limlee_stepctx *l@ = pointer to Lim-Lee context
+ *
+ * Returns:    A @PGEN@ result code.
+ *
+ * Use:                Finalizes the stepper.  The output values in the context
+ *             take on their final results; other resources are discarded.
+ */
+
+static int done(pgen_event *ev, limlee_stepctx *l)
+{
+  size_t i, j;
+  limlee_factor *v;
+
+  /* --- If an output vector of factors is wanted, produce one --- */
+
+  if (!(l->f & LIMLEE_KEEPFACTORS))
+    v = 0;
+  else {
+    if (l->qq.p)
+      l->nf++;
+    v = xmalloc(l->nf * sizeof(limlee_factor));
+  }
+
+  for (i = 0, j = 0; i < l->poolsz; i++) {
+    if (v && l->c[i])
+      v[j++] = l->v[i];
+    else if (l->v[i].p)
+      l->pops->pfree(&l->v[i], l);
+  }
+
+  if (l->qq.p) {
+    if (v)
+      v[j++] = l->qq;
+    else
+      l->pops->pfree(&l->qq, l);
+  }
+
+  xfree(l->v);
+  l->v = v;
+
+  /* --- Free other resources --- */
+
+  xfree(l->c);
+  dstr_destroy(&l->d);
+
+  /* --- Done --- */
+
+  return (PGEN_DONE);
+}
+
+/* --- @limlee_step@ --- */
+
+int limlee_step(int rq, pgen_event *ev, void *p)
+{
+  limlee_stepctx *l = p;
+  int rc;
+
+  switch (rq) {
+    case PGEN_BEGIN:
+      if ((rc = init(ev, l)) != PGEN_TRY)
+       return (rc);
+    case PGEN_TRY:
+      return (next(rq, ev, l));
+    case PGEN_DONE:
+      return (done(ev, l));
+  }
+  return (PGEN_ABORT);
+}
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- @limlee@ --- *
+ *
+ * Arguments:  @const char *name@ = pointer to name root
+ *             @mp *d@ = pointer to destination integer
+ *             @mp *newp@ = how to generate factor primes
+ *             @unsigned ql@ = size of individual factors
+ *             @unsigned pl@ = size of large prime
+ *             @grand *r@ = a random number source
+ *             @unsigned on@ = number of outer attempts to make
+ *             @pgen_proc *oev@ = outer event handler function
+ *             @void *oec@ = argument for the outer event handler
+ *             @pgen_proc *iev@ = inner event handler function
+ *             @void *iec@ = argument for the inner event handler
+ *             @size_t *nf@, @mp ***f@ = output array for factors
+ *
+ * Returns:    A Lim-Lee prime, or null if generation failed.
+ *
+ * Use:                Generates Lim-Lee primes.  A Lim-Lee prime %$p$% is one which
+ *             satisfies %$p = 2 \prod_i q_i + 1$%, where all of the %$q_i$%
+ *             are large enough to resist square-root discrete log
+ *             algorithms.
+ *
+ *             If we succeed, and @f@ is non-null, we write the array of
+ *             factors chosen to @f@ for the benefit of the caller.
+ */
+
+mp *limlee(const char *name, mp *d, mp *newp,
+          unsigned ql, unsigned pl, grand *r,
+          unsigned on, pgen_proc *oev, void *oec,
+          pgen_proc *iev, void *iec,
+          size_t *nf, mp ***f)
+{
+  limlee_stepctx l;
+  rabin rr;
+
+  l.f = 0; if (f) l.f |= LIMLEE_KEEPFACTORS;
+  l.newp = newp;
+  l.pl = pl; l.ql = ql;
+  l.pops = 0;
+  l.iev = iev;
+  l.iec = iec;
+  l.r = r;
+
+  d = pgen(name, d, 0, oev, oec, on, limlee_step, &l,
+          rabin_iters(pl), pgen_test, &rr);
+
+  if (d && f) {
+    mp **v;
+    size_t i;
+    v = xmalloc(l.nf * sizeof(mp *));
+    for (i = 0; i < l.nf; i++)
+      v[i] = l.v[i].p;
+    xfree(l.v);
+    *f = v;
+    *nf = l.nf;
+  }
+
+  return (d);
+}
+
+/*----- That's all, folks -------------------------------------------------*/