math/scaf.c: Add some debugging utilities I found handy.

[catacomb] / math / limlee.c

/* -*-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->u.s.name, 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;

  /* --- First of all, decide on a number of factors to make --- */

  l->nf = l->pl / l->ql;
  if (l->nf < 2) return (PGEN_ABORT);
  l->nf--;

  /* --- 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;
  if (!l->pops) {
    l->pops = &primeops_simple;
    l->pc = 0;
  }

  /* --- Find a big prime later --- */

  l->qq.p = 0;

  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)
{
  dstr d = DSTR_INIT;
  mp *p;
  int rc;
  int dist;
  unsigned nb;

  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 */

    /* --- 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);
    for (i = 0; i < l->poolsz; i++) {
      if (!l->c[i])
	continue;
      if (!l->v[i].p) {
	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 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);
}

/* --- @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);

  /* --- 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 -------------------------------------------------*/