math/gfreduce.[ch]: Fix out-of-bounds memory access.

[u/mdw/catacomb] / math / strongprime.c

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

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

#include <mLib/dstr.h>

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

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

/* --- @strongprime_setup@ --- *
 *
 * Arguments:   @const char *name@ = pointer to name root
 *              @mp *d@ = destination for search start point
 *              @pfilt *f@ = where to store filter jump context
 *              @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 starting point for a `strong' prime search, or zero.
 *
 * Use:         Sets up for a strong prime search, so that primes with
 *              particular properties can be found.  It's probably important
 *              to note that the number left in the filter context @f@ is
 *              congruent to 2 (mod 4).
 */

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

  mp *rr = d;
  pgen_filterctx c;
  pgen_jumpctx j;
  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 12

  /* --- 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;
  c.step = 1;

  rr = mprand(rr, nbits, 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)
    goto fail_s;

  rr = mprand(rr, nbits, 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)
    goto fail_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);
  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);
  pfilt_destroy(&c.f);
  if (!q)
    goto fail_r;

  /* --- 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(f, x);
    x = mp_lsl(x, x, BITSLOP - 1);
    rr = mp_add(rr, rr, x);
    mp_drop(x);
  }

  /* --- Return the result --- */

  mp_drop(q);
  mp_drop(t);
  mp_drop(s);
  dstr_destroy(&dn);
  return (rr);

  /* --- Tidy up if something failed --- */

fail_r:
  mp_drop(t);
fail_t:
  mp_drop(s);
fail_s:
  mp_drop(rr);
  dstr_destroy(&dn);
  return (0);

#undef BITSLOP
}

/* --- @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)
{
  pfilt f;
  pgen_jumpctx j;
  rabin rb;

  d = strongprime_setup(name, d, &f, nbits, r, n, event, ectx);
  j.j = &f;
  d = pgen(name, d, d, event, ectx, n, pgen_jump, &j,
           rabin_iters(nbits), pgen_test, &rb);
  pfilt_destroy(&f);
  return (d);
}

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