tiger-mktab.c: Don't have printf swallow a kludge64 whole.
[u/mdw/catacomb] / pgen.c
diff --git a/pgen.c b/pgen.c
index a57c30a..1a2a8ac 100644 (file)
--- a/pgen.c
+++ b/pgen.c
@@ -1,13 +1,13 @@
 /* -*-c-*-
  *
- * $Id: pgen.c,v 1.2 1999/11/20 22:23:05 mdw Exp $
+ * $Id$
  *
- * Finding and testing prime numbers
+ * Prime generation glue
  *
  * (c) 1999 Straylight/Edgeware
  */
 
-/*----- Licensing notice --------------------------------------------------* 
+/*----- Licensing notice --------------------------------------------------*
  *
  * This file is part of Catacomb.
  *
  * 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: pgen.c,v $
- * Revision 1.2  1999/11/20 22:23:05  mdw
- * Add multiply-and-add function for Diffie-Hellman safe prime generation.
- *
- * Revision 1.1  1999/11/19 13:17:57  mdw
- * Prime number generator and tester.
- *
- */
-
 /*----- Header files ------------------------------------------------------*/
 
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "fibrand.h"
+#include "grand.h"
 #include "mp.h"
-#include "mpmont.h"
+#include "mprand.h"
 #include "pgen.h"
-#include "ptab.h"
+#include "pfilt.h"
 #include "rabin.h"
 
-/*----- Main code ---------------------------------------------------------*/
+/*----- Standard prime filter ---------------------------------------------*/
 
-/* --- @pgen_create@ --- *
- *
- * Arguments:  @pgen *p@ = pointer to prime generation context
- *             @mp *m@ = pointer to initial number to test
- *
- * Returns:    One of the @PGEN@ constants above.
- *
- * Use:                Tests an initial number for primality by computing its
- *             residue modulo various small prime numbers.  This is fairly
- *             quick, but not particularly certain.  If a @PGEN_MAYBE@
- *             result is returned, perform Rabin-Miller tests to confirm.
- */
+/* --- @pgen_filter@ --- */
 
-int pgen_create(pgen *p, mp *m)
+int pgen_filter(int rq, pgen_event *ev, void *p)
 {
-  int rc = PGEN_MAYBE;
-  int i;
-  mp *r = MP_NEW;
-  mpw qw;
-  mp q;
-
-  /* --- Take a copy of the number --- */
-
-  mp_shrink(m);
-  p->m = MP_COPY(m);
-
-  /* --- Fill in the residues --- */
-
-  mp_build(&q, &qw, &qw + 1);
-  for (i = 0; i < NPRIME; i++) {
-    qw = ptab[i];
-    mp_div(0, &r, m, &q);
-    p->r[i] = r->v[0];
-    if (!p->r[i] && rc == PGEN_MAYBE) {
-      if (MP_LEN(m) == 1 && m->v[0] == ptab[i])
-       rc = PGEN_PRIME;
-      else
-       rc = PGEN_COMPOSITE;
-    }
-  }
+  pgen_filterctx *f = p;
+  int rc = PGEN_FAIL;
 
-  /* --- Done --- */
+  switch (rq) {
+    case PGEN_BEGIN:
+      rc = pfilt_create(&f->f, ev->m);
+      mp_drop(ev->m);
+      break;
+    case PGEN_TRY:
+      mp_drop(ev->m);
+      break;
+    case PGEN_DONE:
+      pfilt_destroy(&f->f);
+      return (PGEN_DONE);
+    default:
+      rc = PGEN_ABORT;
+      break;
+  }
 
-  mp_drop(r);
+  if (rc == PGEN_FAIL && !((f->step | f->f.m->v[0]) & 1))
+    rc = pfilt_step(&f->f, 1);
+  while (rc == PGEN_FAIL)
+    rc = pfilt_step(&f->f, f->step);
+  ev->m = MP_COPY(f->f.m);
   return (rc);
 }
 
-/* --- @pgen_destroy@ --- *
- *
- * Arguments:  @pgen *p@ = pointer to prime generation context
- *
- * Returns:    ---
+/* --- @pgen_jump@ --- *
  *
- * Use:                Discards a context and all the resources it holds.
+ * Similar to the standard @pgen_filter@, but jumps in large steps rather
+ * than small ones.
  */
 
-void pgen_destroy(pgen *p)
+int pgen_jump(int rq, pgen_event *ev, void *p)
 {
-  mp_drop(p->m);
-}
-
-/* --- @pgen_step@ --- *
- *
- * Arguments:  @pgen *p@ = pointer to prime generation context
- *             @mpw step@ = how much to step the number
- *
- * Returns:    One of the @PGEN@ constants above.
- *
- * Use:                Steps a number by a small amount.  Stepping is much faster
- *             than initializing with a new number.  The test performed is
- *             the same simple one used by @ptab_create@, so @PGEN_MAYBE@
- *             results should be followed up by a Rabin-Miller test.
- */
+  pgen_jumpctx *f = p;
+  int rc = PGEN_ABORT;
+
+  switch (rq) {
+    case PGEN_BEGIN: {
+      mp *g = MP_NEW;
+      mp_gcd(&g, 0, 0, ev->m, f->j->m);
+      if (MP_CMP(g, >, MP_ONE)) {
+       mp_drop(g);
+       return (PGEN_ABORT);
+      }
+      mp_drop(g);
+      rc = pfilt_create(&f->f, ev->m);
+      mp_drop(ev->m);
+    } break;
+    case PGEN_TRY:
+      mp_drop(ev->m);
+      rc = pfilt_jump(&f->f, f->j);
+      break;
+    case PGEN_DONE:
+      pfilt_destroy(&f->f);
+      return (PGEN_DONE);
+  }
 
-int pgen_step(pgen *p, mpw step)
-{
-  int rc = PGEN_MAYBE;
-  int i;
+  while (rc == PGEN_FAIL)
+    rc = pfilt_jump(&f->f, f->j);
+  ev->m = MP_COPY(f->f.m);
+  return (rc);
+}
 
-  /* --- Add the step on to the number --- */
+/*----- Standard prime test -----------------------------------------------*/
 
-  p->m = mp_modify(p->m, MP_LEN(p->m) + 1);
-  mpx_uaddn(p->m->v, p->m->vl, step);
-  mp_shrink(p->m);
+/* --- @pgen_test@ --- */
 
-  /* --- Update the residue table --- */
+int pgen_test(int rq, pgen_event *ev, void *p)
+{
+  rabin *r = p;
+  int rc = PGEN_ABORT;
 
-  for (i = 0; i < NPRIME; i++) {
-    p->r[i] = (p->r[i] + step) % ptab[i];
-    if (!p->r[i] && rc == PGEN_MAYBE) {
-      if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
-       rc = PGEN_PRIME;
-      else
-       rc = PGEN_COMPOSITE;
-    }
+  switch (rq) {
+    case PGEN_BEGIN:
+      rabin_create(r, ev->m);
+      rc = PGEN_TRY;
+      break;
+    case PGEN_TRY:
+      if (!ev->tests)
+       rc = rabin_rtest(r, MP_TWO);
+      else {
+       mp *a = mprand_range(MP_NEW, ev->m, ev->r, 0);
+       rc = rabin_rtest(r, a);
+       mp_drop(a);
+      }
+      break;
+    case PGEN_DONE:
+      rabin_destroy(r);
+      rc = PGEN_DONE;
+      break;
   }
 
-  /* --- Small numbers must be prime --- */
-
-  if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
-      p->m->v[0] < MAXPRIME * MAXPRIME)
-    rc = PGEN_PRIME;
-
-  /* --- Done --- */
-
   return (rc);
 }
 
-/* --- @pgen_muladd@ --- *
- *
- * Arguments:  @pgen *p@ = destination prime generation context
- *             @const pgen *q@ = source prime generation context
- *             @mpw m@ = number to multiply by
- *             @mpw a@ = number to add
+/*----- The main driver ---------------------------------------------------*/
+
+/* --- @pgen@ --- *
  *
- * Returns:    One of the @PGEN@ constants above.
+ * Arguments:  @const char *name@ = name of the value being searched for
+ *             @mp *d@ = destination for the result integer
+ *             @mp *m@ = start value to pass to stepper
+ *             @pgen_proc *event@ = event handler function
+ *             @void *ectx@ = context argument for event andler
+ *             @unsigned steps@ = number of steps to take in search
+ *             @pgen_proc *step@ = stepper function to use
+ *             @void *sctx@ = context argument for stepper
+ *             @unsigned tests@ = number of tests to make
+ *             @pgen_proc *test@ = tester function to use
+ *             @void *tctx@ = context argument for tester
  *
- * Use:                Multiplies the number in a prime generation context by a
- *             small value and then adds a small value.  The destination
- *             should either be uninitialized or the same as the source.
+ * Returns:    Pointer to final result, or null.
  *
- *             Common things to do include multiplying by 2 and adding 0 to
- *             turn a prime into a jump for finding other primes with @q@ as
- *             a factor of @p - 1@, or multiplying by 2 and adding 1.
+ * Use:                A generalized prime-number search skeleton.  Yes, that's a
+ *             scary number of arguments.
  */
 
-int pgen_muladd(pgen *p, const pgen *q, mpw m, mpw a)
+mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
+        unsigned steps, pgen_proc *step, void *sctx,
+        unsigned tests, pgen_proc *test, void *tctx)
 {
-  int rc = PGEN_MAYBE;
-  int i;
-
-  /* --- Multiply the big number --- */
-
-  {
-    mp *d = mp_create(MP_LEN(q->m) + 2);
-    mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m);
-    mpx_uaddn(d->v, d->vl, a);
-    d->f = q->m->f;
-    if (p == q)
-      mp_drop(p->m);
-    mp_shrink(d);
-    p->m = d;
+  pgen_event ev;
+  int rq, rc;
+  pgen_proc *proc;
+  void *ctx;
+  int p;
+
+  enum { P_STEP, P_TEST };
+
+  /* --- Set up the initial event block --- */
+
+  ev.name = name;
+  if (m)
+    ev.m = MP_COPY(m);
+  else
+    ev.m = 0;
+  ev.steps = 0;
+  ev.tests = 0;
+  ev.r = fibrand_create(0);
+
+  /* --- Tell the event handler we're under way --- */
+
+  if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) {
+    ev.r->ops->destroy(ev.r);
+    return (0);
   }
 
-  /* --- Gallivant through the residue table --- */
-      
-  for (i = 0; i < NPRIME; i++) {
-    p->r[i] = (q->r[i] * m + a) % ptab[i];
-    if (!p->r[i] && rc == PGEN_MAYBE) {
-      if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
-       rc = PGEN_PRIME;
-      else
-       rc = PGEN_COMPOSITE;
+  /* --- Set up for the initial call --- */
+
+  proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN;
+
+  /* --- Enter the great maelstrom of state transitions --- */
+
+  for (;;) {
+    unsigned act = 0;
+
+#define A_STEP 1u
+#define A_TEST 2u
+#define A_EVENT 4u
+#define A_ENDTEST 8u
+#define A_ENDSTEP 16u
+#define A_DONE 32u
+
+    /* --- Call the procedure and decide what to do next --- */
+
+    rc = proc(rq, &ev, ctx);
+    switch (rc) {
+      case PGEN_TRY:
+       if (p == P_TEST)
+         rq = PGEN_TRY;
+       else {
+         act |= A_EVENT;
+         proc = test; ctx = tctx; p = P_TEST;
+         rq = PGEN_BEGIN;
+       }
+       break;
+      case PGEN_PASS:
+       act |= A_TEST | A_EVENT;
+       if (p == P_TEST)
+         rq = PGEN_TRY;
+       else {
+         proc = test; ctx = tctx; p = P_TEST;
+         rq = PGEN_BEGIN;
+       }
+       break;
+      case PGEN_FAIL:
+       act |= A_STEP;
+       if (p == P_TEST) {
+         act |= A_ENDTEST | A_EVENT;
+         proc = step; ctx = sctx; p = P_STEP;
+       }
+       rq = PGEN_TRY;
+       break;
+      case PGEN_DONE:
+       act |= A_EVENT | A_DONE | A_ENDSTEP;
+       if (p == P_TEST)
+         act |= A_ENDTEST;
+       break;
+      case PGEN_ABORT:
+       act |= A_EVENT | A_DONE;
+       if (p == P_TEST || rq == PGEN_TRY)
+         act |= A_ENDSTEP;
+       if (p == P_TEST && rq != PGEN_BEGIN)
+         act |= A_ENDTEST;
+       break;
+      default:
+       assert(((void)"Invalid response from function", 0));
+       break;
     }
-  }
 
-  /* --- Small numbers must be prime --- */
+    /* --- If decrementing counters is requested, do that --- */
 
-  if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
-      p->m->v[0] < MAXPRIME * MAXPRIME)
-    rc = PGEN_PRIME;
+    if ((act & A_STEP) && steps) {
+      ev.steps++;
+      if (ev.steps == steps) {
+       act |= A_EVENT | A_ENDSTEP | A_DONE;
+       rc = PGEN_ABORT;
+      }
+      ev.tests = 0;
+    }
 
-  /* --- Finished --- */
+    if ((act & A_TEST) && tests) {
+      ev.tests++;
+      if (ev.tests == tests) {
+       act |= A_ENDTEST | A_ENDSTEP | A_DONE;
+       rc = PGEN_DONE;
+      }
+    }
 
-  return (rc);
-}
+    /* --- Report an event if so directed --- */
 
+    if ((act & A_EVENT) && event && event(rc, &ev, ectx) == PGEN_ABORT) {
+      rc = PGEN_ABORT;
+      if (!(act & A_DONE)) {
+       act |= A_ENDSTEP | A_DONE;
+       if (p == P_TEST)
+         act |= A_ENDTEST;
+      }
+    }
 
-/* --- @pgen_jump@ --- *
- *
- * Arguments:  @pgen *p@ = pointer to prime generation context
- *             @pgen *j@ = pointer to another generation context
- *
- * Returns:    One of the @PGEN@ constants above.
- *
- * Use:                Steps a number by a large amount.  Even so, jumping is much
- *             faster than initializing a new number.  The test peformed is
- *             the same simple one used by @ptab_create@, so @PGEN_MAYBE@
- *             results should be followed up by a Rabin-Miller test.
- *
- *             Note that the number stored in the @j@ context is probably
- *             better off being even than prime.  The important thing is
- *             that all of the residues for the number have already been
- *             computed.
- */
+    /* --- Close down tester and stepper functions --- */
 
-int pgen_jump(pgen *p, pgen *j)
-{
-  int rc = PGEN_MAYBE;
-  int i;
+    if (act & A_ENDTEST)
+      test(PGEN_DONE, &ev, tctx);
+    if (act & A_ENDSTEP)
+      step(PGEN_DONE, &ev, sctx);
 
-  /* --- Add the step on --- */
+    /* --- Stop the entire test if necessary --- */
 
-  p->m = mp_add(p->m, p->m, j->m);
+    if (act & A_DONE)
+      break;
+  }
 
-  /* --- Update the residue table --- */
+  /* --- Tidy up and return --- */
 
-  for (i = 0; i < NPRIME; i++) {
-    p->r[i] = p->r[i] + j->r[i];
-    if (p->r[i] > ptab[i])
-      p->r[i] -= ptab[i];
-    if (!p->r[i] && rc == PGEN_MAYBE) {
-      if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
-       rc = PGEN_PRIME;
-      else
-       rc = PGEN_COMPOSITE;
-    }
+  if (rc == PGEN_ABORT) {
+    mp_drop(ev.m);
+    ev.m = 0;
   }
+  ev.r->ops->destroy(ev.r);
+  mp_drop(d);
 
-  /* --- Small numbers must be prime --- */
+  return (ev.m);
+}
 
-  if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
-      p->m->v[0] < MAXPRIME * MAXPRIME)
-    rc = PGEN_PRIME;
+/* --- @pgen_primep@ --- *
+ *
+ * Arguments:  @mp *p@ = a number to check
+ *             @grand *gr@ = a random number source
+ *
+ * Returns:    Nonzero if @p@ is really prime.
+ */
 
-  /* --- Done --- */
+int pgen_primep(mp *p, grand *gr)
+{
+  int i;
+  rabin r;
+  mp *x = MP_NEW;
 
-  return (rc);
+  if (MP_NEGP(p)) return (0);
+  switch (pfilt_smallfactor(p)) {
+    case PGEN_DONE: return (1);
+    case PGEN_FAIL: return (0);
+  }
+  rabin_create(&r, p);
+  for (i = 32; i; i--) {
+    x = mprand_range(x, p, gr, 0);
+    if (rabin_rtest(&r, x) == PGEN_FAIL)
+      break;
+  }
+  MP_DROP(x);
+  rabin_destroy(&r);
+  return (!i);
 }
 
-/*----- Test code ---------------------------------------------------------*/
+/*----- Test rig ----------------------------------------------------------*/
 
 #ifdef TEST_RIG
 
 #include <mLib/testrig.h>
 
-static int verify(dstr *v)
+static int t_primep(dstr *v)
 {
   mp *m = *(mp **)v[0].buf;
-  mp *r = *(mp **)v[1].buf;
-  pgen p;
-  mpw gw;
-  mp g;
-  int rc;
-  static char baton[] = "/-\\|";
-  char *b = baton;
+  int e = *(int *)v[1].buf;
+  int r;
+  grand *rng;
   int ok = 1;
 
-  mp_build(&g, &gw, &gw + 1);
-  rc = pgen_create(&p, m);
-  for (;;) {
-    if (rc == PGEN_PRIME)
-      break;
-    if (rc == PGEN_MAYBE) {
-      int i;
-      rabin r;
-      rabin_create(&r, p.m);
-      for (i = 0; i < 5; i++) {
-       gw = rand();
-       if (rabin_test(&r, &g) == PGEN_COMPOSITE)
-         break;
-      }
-      rabin_destroy(&r);
-      if (i == 5)
-       break;
-      putc(*b++, stderr);
-      putc('\b', stderr);
-      if (!*b)
-       b = baton;
-    }
-    rc = pgen_step(&p, 2);
+  rng = fibrand_create(0);
+  r = pgen_primep(m, rng);
+  GR_DESTROY(rng);
+  if (e != r) {
+    fputs("\n*** primep failed", stderr);
+    fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
+    fprintf(stderr, "\nexpected %d", e);
+    fprintf(stderr, "\nreported %d", r);
+    fputc('\n', stderr);
+    ok = 0;
   }
 
-  if (MP_CMP(p.m, !=, r)) {
-    fputs("\n*** failed.", stderr);
+  mp_drop(m);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
+  return (ok);
+}
+
+static int verify(dstr *v)
+{
+  mp *m = *(mp **)v[0].buf;
+  mp *q = *(mp **)v[1].buf;
+  mp *p;
+  int ok = 1;
+
+  pgen_filterctx pf;
+  rabin r;
+
+  pf.step = 2;
+  p = pgen("p", MP_NEW, m, pgen_evspin, 0, 0, pgen_filter, &pf,
+          rabin_iters(mp_bits(m)), pgen_test, &r);
+  if (!p || !MP_EQ(p, q)) {
+    fputs("\n*** pgen failed", stderr);
     fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
-    fputs("\np = ", stderr); mp_writefile(p.m, stderr, 10);
-    fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
+    fputs("\np = ", stderr); mp_writefile(p, stderr, 10);
+    fputs("\nq = ", stderr); mp_writefile(q, stderr, 10);
     fputc('\n', stderr);
     ok = 0;
   }
 
   mp_drop(m);
-  mp_drop(r);
-  pgen_destroy(&p);
+  mp_drop(q);
+  mp_drop(p);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 
 static test_chunk tests[] = {
   { "pgen", verify, { &type_mp, &type_mp, 0 } },
+  { "primep", t_primep, { &type_mp, &type_int, 0 } },
   { 0, 0, { 0 } }
 };
 
@@ -337,7 +417,6 @@ int main(int argc, char *argv[])
   test_run(argc, argv, tests, SRCDIR "/tests/pgen");
   return (0);
 }
-
 #endif
 
 /*----- That's all, folks -------------------------------------------------*/