(pfilt_smallfactor): New function for doing trial division the hard

[u/mdw/catacomb] / limlee.c

/* -*-c-*-
 *
 * $Id: limlee.c,v 1.3 2000/07/29 09:58:32 mdw Exp $
 *
 * 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.
 */

/*----- Revision history --------------------------------------------------* 
 *
 * $Log: limlee.c,v $
 * Revision 1.3  2000/07/29 09:58:32  mdw
 * (limlee): Bug fix.  Old versions didn't set the filter step if @ql@ was
 * an exact divisor of @pl@.
 *
 * Revision 1.2  2000/07/26 18:00:00  mdw
 * No footer line!
 *
 * Revision 1.1  2000/07/09 21:30:58  mdw
 * Lim-Lee prime generation.
 *
 */

/*----- Header files ------------------------------------------------------*/

#include <mLib/alloc.h>
#include <mLib/dstr.h>

#include "limlee.h"
#include "mpmul.h"
#include "mprand.h"
#include "pgen.h"
#include "primorial.h"
#include "rabin.h"

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

static void comb_init(octet *c, unsigned n, unsigned r)
{
  memset(c, 0, n - r);
  memset(c + (n - r), 1, r);
}

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);
}

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)
{
  dstr dn = DSTR_INIT;
  unsigned qql;
  mp *qq = 0;
  unsigned nn;
  unsigned mm;
  mp **v;
  octet *c;
  unsigned i;
  unsigned long seq = 0;
  pgen_event ev;
  unsigned ntest;
  rabin rb;
  pgen_filterctx pf;

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

  nn = pl/ql;
  qql = pl%ql;
  if (!nn)
    return (0);
  else if (qql && nn > 1) {
    nn--;
    qql += 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.
   */

  mm = nn * 3 + 5;
  if (mm < 25)
    mm = 25;

  /* --- Now allocate the working memory --- */

  primorial_setup();
  v = xmalloc(mm * sizeof(mp *));
  c = xmalloc(mm);

  /* --- Initialize everything and try to find a prime --- */

  ev.name = name;
  ev.m = 0;
  ev.steps = on;
  ev.tests = ntest = rabin_iters(pl);
  ev.r = r;

  if (oev && oev(PGEN_BEGIN, &ev, oec) == PGEN_ABORT)
    goto fail;

  pf.step = 2;
  if (qql) {
    dstr_putf(&dn, "%s [+]", name);
    qq = mprand(d, qql, r, 1);
    qq = pgen(dn.buf, qq, qq, iev, iec,
	      0, pgen_filter, &pf, rabin_iters(qql), pgen_test, &rb);
  }

again:
  comb_init(c, mm, nn);
  for (i = 0; i < mm; i++)
    v[i] = 0;

  /* --- The main combinations loop --- */

  do {
    mpmul mmul = MPMUL_INIT;

    /* --- Multiply a bunch of primes together --- */

    if (qq) {
      mpmul_add(&mmul, qq);
    }
    for (i = 0; i < mm; i++) {
      if (!c[i])
	continue;
      if (!v[i]) {
	mp *z;

	DRESET(&dn);
	dstr_putf(&dn, "%s [%lu] = ", name, seq++);
	z = mprand(newp, ql, ev.r, 1);
	z = pgen(dn.buf, z, z, iev, iec,
		 0, pgen_filter, &pf, rabin_iters(ql), pgen_test, &rb);
	v[i] = z;
      }
      mpmul_add(&mmul, v[i]);
    }

    /* --- Now do some testing --- */

    {
      mp *p = mpmul_done(&mmul);
      mp *g = newp;
      int rc;

      /* --- Check for small factors --- */

      p = mp_lsl(p, p, 1);
      p = mp_add(p, p, MP_ONE);
      mp_gcd(&g, 0, 0, p, primorial);
      if (MP_CMP(g, !=, MP_ONE)) {
	mp_drop(g);
	mp_drop(p);
	continue;
      }
      mp_drop(g);

      /* --- Send an event out --- */

      ev.m = p;
      if (oev && oev(PGEN_TRY, &ev, oec) == PGEN_ABORT) {
	mp_drop(p);
	goto fail;
      }

      /* --- Do the Rabin testing --- */

      rabin_create(&rb, p);
      g = MP_NEW;
      do {
	g = mprand_range(g, p, ev.r, 1);
	rc = rabin_test(&rb, g);
	if (rc == PGEN_PASS) {
	  ev.tests--;
	  if (!ev.tests)
	    rc = PGEN_DONE;
	}
	if (oev &&oev(rc, &ev, oec) == PGEN_ABORT)
	  rc = PGEN_ABORT;
      } while (rc == PGEN_PASS);

      rabin_destroy(&rb);
      mp_drop(g);
      if (rc == PGEN_DONE)
	d = p;
      else
	mp_drop(p);
      if (rc == PGEN_ABORT)
	goto fail;
      if (rc == PGEN_DONE)
	goto done;
      ev.tests = ntest;
      ev.m = 0;
    }
  } while (comb_next(c, mm, nn));

  /* --- That failed --- */

  if (ev.steps) {
    ev.steps--;
    if (!ev.steps) {
      if (oev)
	oev(PGEN_ABORT, &ev, &oec);
      goto fail;
    }
  }

  for (i = 0; i < mm; i++)
    mp_drop(v[i]);
  goto again;

  /* --- We did it! --- */

done: {
    mp **vv = 0;
    if (f) {
      if (qq)
	nn++;
      *nf = nn;
      *f = vv = xmalloc(nn * sizeof(mp *));
    }

    for (i = 0; i < mm; i++) {
      if (c[i] && vv)
	*vv++ = v[i];
      else if (v[i])
	mp_drop(v[i]);
    }
    if (qq) {
      if (vv)
	*vv++ = qq;
      else
	mp_drop(qq);
    }
    xfree(v);
    xfree(c);
    dstr_destroy(&dn);
    return (d);
  }

  /* --- We blew it --- */

fail:
  for (i = 0; i < mm; i++)
    mp_drop(v[i]);
  if (qq)
    mp_drop(qq);
  xfree(v);
  xfree(c);
  dstr_destroy(&dn);
  return (0);
}

/*----- That's all, folks -------------------------------------------------*/