Find `safe' primes (i.e., %$p = 2q + 1$%).

[u/mdw/catacomb] / strongprime.c

/* -*-c-*-
 *
 * $Id: strongprime.c,v 1.1 1999/12/22 15:51:22 mdw Exp $
 *
 * Generate `strong' prime numbers
 *
 * (c) 1999 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: strongprime.c,v $
 * Revision 1.1  1999/12/22 15:51:22  mdw
 * Find `strong' RSA primes using Gordon's algorithm.
 *
 */

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

#include <mLib/dstr.h>

#include "grand.h"
#include "rand.h"
#include "mp.h"
#include "mpmont.h"
#include "mprand.h"
#include "pgen.h"
#include "pfilt.h"
#include "rabin.h"

/*----- Main code ---------------------------------------------------------*/

/* --- @strongprime@ --- *
 *
 * Arguments:	@const char *name@ = pointer to name root
 *		@mp *d@ = destination integer
 *		@unsigned nbits@ = number of bits wanted
 *		@grand *r@ = random number source
 *		@unsigned n@ = number of attempts to make
 *		@pgen_proc *event@ = event handler function
 *		@void *ectx@ = argument for the event handler
 *
 * Returns:	A `strong' prime, or zero.
 *
 * Use:		Finds `strong' primes.  A strong prime %$p$% is such that
 *
 *		  * %$p - 1$% has a large prime factor %$r$%,
 *		  * %$p + 1$% has a large prime factor %$s$%, and
 *		  * %$r - 1$% has a large prime factor %$t$%.
 *
 *		The numbers produced may be slightly larger than requested,
 *		by a few bits.
 */

mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r,
		unsigned n, pgen_proc *event, void *ectx)
{
  mp *s, *t, *q, *p = 0;
  dstr dn = DSTR_INIT;

  mp *rr = MP_NEW;
  pgen_filterctx c;
  pgen_jumpctx cj;
  rabin rb;

  /* --- The bitslop parameter --- *
   *
   * There's quite a lot of prime searching to be done.  The constant
   * @BITSLOP@ is a (low) approximation to the base-2 log of the expected
   * number of steps to find a prime number.  Experimentation shows that
   * numbers around 10 seem to be good.
   */

#define BITSLOP 10

  /* --- Choose two primes %$s$% and %$t$% of half the required size --- */

  nbits = nbits/2 - BITSLOP;
  c.step = 1;

  rr = mprand(rr, nbits, r, 1);
  DRESET(&dn); dstr_putf(&dn, "%s [s]", name);
  if ((s = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_filter, &c,
		rabin_iters(nbits), pgen_test, &rb)) == 0)
    goto fail_s;
  mp_burn(s);

  rr = mprand(rr, nbits, r, 1);
  DRESET(&dn); dstr_putf(&dn, "%s [t]", name);
  if ((t = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_filter, &c,
		rabin_iters(nbits), pgen_test, &rb)) == 0)
    goto fail_t;
  mp_burn(t);

  /* --- Choose a suitable value for %$r = 2it + 1$% for some %$i$% --- */

  rr = mp_lsl(rr, t, 1);
  pfilt_create(&c.f, rr);
  rr = mp_lsl(rr, rr, BITSLOP - 1);
  rr = mp_add(rr, rr, MP_ONE);
  DRESET(&dn); dstr_putf(&dn, "%s [r]", name);
  cj.j = &c.f;
  nbits += BITSLOP;
  if ((q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &cj,
		rabin_iters(nbits), pgen_test, &rb)) == 0)
    goto fail_r;
  pfilt_destroy(&c.f);

  /* --- Select a suitable starting-point for finding %$p$% --- *
   *
   * This computes %$p_0 = 2(s^{r - 2} \bmod r)s - 1$%.
   */

  {
    mpmont mm;

    mpmont_create(&mm, q);
    rr = mp_sub(rr, q, MP_TWO);
    rr = mpmont_exp(&mm, rr, s, rr);
    mpmont_destroy(&mm);
    rr = mp_mul(rr, rr, s);
    rr = mp_lsl(rr, rr, 1);
    rr = mp_sub(rr, rr, MP_ONE);
  }

  /* --- Now find %$p = p_0 + 2jrs$% for some %$j$% --- */

  {
    mp *x;
    x = mp_mul(MP_NEW, q, s);
    x = mp_lsl(x, x, 1);
    pfilt_create(&c.f, x);
    x = mp_lsl(x, x, BITSLOP - 1);
    rr = mp_add(rr, rr, x);
    mp_drop(x);
  }

  if ((p = pgen(name, d, rr, event, ectx, n, pgen_jump, &cj,
		rabin_iters(nbits * 2), pgen_test, &rb)) == 0)
    goto fail_p;

  /* --- Tidy up because we've finished --- */

fail_p:
  mp_drop(q);
fail_r:
  pfilt_destroy(&c.f);
  mp_drop(t);
fail_t:
  mp_drop(s);
fail_s:
  mp_drop(rr);
  dstr_destroy(&dn);

  return (p);

#undef BITSLOP
}

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