Release 2.1.2.
[u/mdw/catacomb] / pgen.c
diff --git a/pgen.c b/pgen.c
index 74a10bc..1a2a8ac 100644 (file)
--- a/pgen.c
+++ b/pgen.c
@@ -1,13 +1,13 @@
 /* -*-c-*-
  *
- * $Id: pgen.c,v 1.7 2001/02/03 16:05:32 mdw Exp $
+ * $Id$
  *
  * 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.7  2001/02/03 16:05:32  mdw
- * Now @mp_drop@ checks its argument is non-NULL before attempting to free
- * it.  Note that the macro version @MP_DROP@ doesn't do this.
- *
- * Revision 1.6  2000/10/08 12:11:22  mdw
- * Use @MP_EQ@ instead of @MP_CMP@.
- *
- * Revision 1.5  2000/06/17 11:52:36  mdw
- * Signal a pgen abort if the jump and base share a common factor.
- *
- * Revision 1.4  1999/12/22 16:01:11  mdw
- * Same file, completely different code.  Main interface for new prime-
- * search system.
- *
- */
-
 /*----- Header files ------------------------------------------------------*/
 
 #include <assert.h>
@@ -68,7 +49,7 @@
 int pgen_filter(int rq, pgen_event *ev, void *p)
 {
   pgen_filterctx *f = p;
-  int rc = PGEN_ABORT;
+  int rc = PGEN_FAIL;
 
   switch (rq) {
     case PGEN_BEGIN:
@@ -77,16 +58,17 @@ int pgen_filter(int rq, pgen_event *ev, void *p)
       break;
     case PGEN_TRY:
       mp_drop(ev->m);
-      if (!((f->step | f->f.m->v[0]) & 1))
-       rc = pfilt_step(&f->f, 1);
-      else
-       rc = pfilt_step(&f->f, f->step);
       break;
     case PGEN_DONE:
       pfilt_destroy(&f->f);
       return (PGEN_DONE);
+    default:
+      rc = PGEN_ABORT;
+      break;
   }
-       
+
+  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);
@@ -124,7 +106,7 @@ int pgen_jump(int rq, pgen_event *ev, void *p)
       pfilt_destroy(&f->f);
       return (PGEN_DONE);
   }
-       
+
   while (rc == PGEN_FAIL)
     rc = pfilt_jump(&f->f, f->j);
   ev->m = MP_COPY(f->f.m);
@@ -145,11 +127,15 @@ int pgen_test(int rq, pgen_event *ev, void *p)
       rabin_create(r, ev->m);
       rc = PGEN_TRY;
       break;
-    case PGEN_TRY: {
-      mp *a = mprand_range(MP_NEW, ev->m, ev->r, 0);
-      rc = rabin_test(r, a);
-      mp_drop(a);
-    } 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;
@@ -189,6 +175,9 @@ mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
   int rq, rc;
   pgen_proc *proc;
   void *ctx;
+  int p;
+
+  enum { P_STEP, P_TEST };
 
   /* --- Set up the initial event block --- */
 
@@ -197,73 +186,73 @@ mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
     ev.m = MP_COPY(m);
   else
     ev.m = 0;
-  ev.steps = steps;
-  ev.tests = tests;
+  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)
+  if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) {
+    ev.r->ops->destroy(ev.r);
     return (0);
+  }
 
   /* --- Set up for the initial call --- */
 
-  proc = step; ctx = sctx; rq = PGEN_BEGIN;
+  proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN;
 
   /* --- Enter the great maelstrom of state transitions --- */
 
   for (;;) {
     unsigned act = 0;
 
-    enum {
-      A_STEP = 1u,
-      A_TEST = 2u,
-      A_EVENT = 4u,
-      A_ENDTEST = 8u,
-      A_ENDSTEP = 16u,
-      A_DONE = 32u
-    };
+#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 (proc == test)
+       if (p == P_TEST)
          rq = PGEN_TRY;
        else {
          act |= A_EVENT;
-         proc = test; ctx = tctx;
+         proc = test; ctx = tctx; p = P_TEST;
          rq = PGEN_BEGIN;
        }
        break;
       case PGEN_PASS:
        act |= A_TEST | A_EVENT;
-       if (proc == test)
+       if (p == P_TEST)
          rq = PGEN_TRY;
        else {
-         proc = test; ctx = tctx;
+         proc = test; ctx = tctx; p = P_TEST;
          rq = PGEN_BEGIN;
        }
        break;
       case PGEN_FAIL:
        act |= A_STEP;
-       if (proc == test) {
+       if (p == P_TEST) {
          act |= A_ENDTEST | A_EVENT;
-         proc = step; ctx = sctx;
+         proc = step; ctx = sctx; p = P_STEP;
        }
        rq = PGEN_TRY;
        break;
       case PGEN_DONE:
        act |= A_EVENT | A_DONE | A_ENDSTEP;
-       if (proc == test)
+       if (p == P_TEST)
          act |= A_ENDTEST;
        break;
       case PGEN_ABORT:
        act |= A_EVENT | A_DONE;
-       if (proc == test || rq == PGEN_TRY)
+       if (p == P_TEST || rq == PGEN_TRY)
          act |= A_ENDSTEP;
-       if (proc == test && rq == PGEN_BEGIN)
+       if (p == P_TEST && rq != PGEN_BEGIN)
          act |= A_ENDTEST;
        break;
       default:
@@ -274,17 +263,17 @@ mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
     /* --- If decrementing counters is requested, do that --- */
 
     if ((act & A_STEP) && steps) {
-      ev.steps--;
-      if (!ev.steps) {
+      ev.steps++;
+      if (ev.steps == steps) {
        act |= A_EVENT | A_ENDSTEP | A_DONE;
        rc = PGEN_ABORT;
       }
-      ev.tests = tests;
+      ev.tests = 0;
     }
 
     if ((act & A_TEST) && tests) {
-      ev.tests--;
-      if (!ev.tests) {
+      ev.tests++;
+      if (ev.tests == tests) {
        act |= A_ENDTEST | A_ENDSTEP | A_DONE;
        rc = PGEN_DONE;
       }
@@ -296,10 +285,10 @@ mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
       rc = PGEN_ABORT;
       if (!(act & A_DONE)) {
        act |= A_ENDSTEP | A_DONE;
-       if (proc == test)
+       if (p == P_TEST)
          act |= A_ENDTEST;
       }
-    }  
+    }
 
     /* --- Close down tester and stepper functions --- */
 
@@ -326,12 +315,67 @@ mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
   return (ev.m);
 }
 
+/* --- @pgen_primep@ --- *
+ *
+ * Arguments:  @mp *p@ = a number to check
+ *             @grand *gr@ = a random number source
+ *
+ * Returns:    Nonzero if @p@ is really prime.
+ */
+
+int pgen_primep(mp *p, grand *gr)
+{
+  int i;
+  rabin r;
+  mp *x = MP_NEW;
+
+  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 rig ----------------------------------------------------------*/
 
 #ifdef TEST_RIG
 
 #include <mLib/testrig.h>
 
+static int t_primep(dstr *v)
+{
+  mp *m = *(mp **)v[0].buf;
+  int e = *(int *)v[1].buf;
+  int r;
+  grand *rng;
+  int ok = 1;
+
+  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;
+  }
+
+  mp_drop(m);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
+  return (ok);
+}
+
 static int verify(dstr *v)
 {
   mp *m = *(mp **)v[0].buf;
@@ -363,6 +407,7 @@ static int verify(dstr *v)
 
 static test_chunk tests[] = {
   { "pgen", verify, { &type_mp, &type_mp, 0 } },
+  { "primep", t_primep, { &type_mp, &type_int, 0 } },
   { 0, 0, { 0 } }
 };