math/strongprime.c: Replace inexplicable exponentiation with extended-gcd.

[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); that the jump value is twice the
 *              product of two large primes; and that the starting point is
 *              at least %$3 \cdot 2^{N-2}$%.  (Hence, if you multiply two
 *              such numbers, the product is at least
 *
 *                      %$9 \cdot 2^{2N-4} > 2^{2N-1}$%
 *
 *              i.e., it will be (at least) a %$2 N$%-bit value.
 */

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;
  unsigned slop, nb, u, i;

  mp *rr = d;
  pgen_filterctx c;
  pgen_jumpctx j;
  rabin rb;

  /* --- Figure out how large the smaller primes should be --- *
   *
   * We want them to be `as large as possible', subject to the constraint
   * that we produce a number of the requested size at the end.  This is
   * tricky, because the final prime search is going to involve quite large
   * jumps from its starting point; the size of the jumps are basically
   * determined by our choice here, and if they're too big then we won't find
   * a prime in time.
   *
   * Let's suppose we're trying to make an %$N$%-bit prime.  The expected
   * number of steps tends to increase linearly with size, i.e., we need to
   * take about %2^k N$% steps for some %$k$%.  If we're jumping by a
   * %$J$%-bit quantity each time, from an %$N$%-bit starting point, then we
   * will only be able to find a match if %$2^k N 2^{J-1} \le 2^{N-1}$%,
   * i.e., if %$J \le N - (k + \log_2 N)$%.
   *
   * Experimentation shows that taking %$k + \log_2 N = 12$% works well for
   * %$N = 1024$%, so %$k = 2$%.  Add a few extra bits for luck.
   */

  for (i = 1; i && nbits >> i; i <<= 1); assert(i);
  for (slop = 6, nb = nbits; nb > 1; i >>= 1) {
    u = nb >> i;
    if (u) { slop += i; nb = u; }
  }
  if (nbits/2 <= slop) return (0);

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

  nb = nbits/2 - slop;
  c.step = 1;

  rr = mprand(rr, nb, 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(nb), pgen_test, &rb)) == 0)
    goto fail_s;

  rr = mprand(rr, nb, 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(nb), 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, slop - 1);
  rr = mp_add(rr, rr, MP_ONE);
  DRESET(&dn); dstr_putf(&dn, "%s [r]", name);
  j.j = &c.f;
  q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &j,
           rabin_iters(nb + slop), pgen_test, &rb);
  pfilt_destroy(&c.f);
  if (!q)
    goto fail_r;

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

  rr = mp_modinv(rr, s, q);
  rr = mp_mul(rr, rr, s);
  rr = mp_lsl(rr, rr, 1);
  rr = mp_sub(rr, rr, MP_ONE);

  /* --- Pick a starting point for the search --- *
   *
   * Select %$3 \cdot 2^{N-2} < p_1 < 2^N$% at random, only with
   * %$p_1 \equiv p_0 \pmod{2 r s}$.
   */

  {
    mp *x, *y;
    x = mp_mul(MP_NEW, q, s);
    x = mp_lsl(x, x, 1);
    pfilt_create(f, x); /* %$2 r s$% */
    y = mprand(MP_NEW, nbits, r, 0);
    y = mp_setbit(y, y, nbits - 2);
    rr = mp_leastcongruent(rr, y, rr, x);
    mp_drop(x); mp_drop(y);
  }

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

/* --- @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$%.
 */

mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r,
                unsigned n, pgen_proc *event, void *ectx)
{
  mp *p;
  pfilt f;
  pgen_jumpctx j;
  rabin rb;

  if (d) mp_copy(d);
  p = strongprime_setup(name, d, &f, nbits, r, n, event, ectx);
  if (!p) { mp_drop(d); return (0); }
  j.j = &f;
  p = pgen(name, p, p, event, ectx, n, pgen_jump, &j,
           rabin_iters(nbits), pgen_test, &rb);
  if (mp_bits(p) != nbits) { mp_drop(p); return (0); }
  pfilt_destroy(&f);
  mp_drop(d);
  return (p);
}

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