math/pgen.c, math/pgen-simul.c: Add Baillie--PSW testers.

[catacomb] / math / pgen.c

diff --git a/math/pgen.c b/math/pgen.c
index f10d585..b7163b9 100644 (file)
--- a/math/pgen.c
+++ b/math/pgen.c
@@ -118,6 +118,7 @@ int pgen_jump(int rq, pgen_event *ev, void *p)
 int pgen_test(int rq, pgen_event *ev, void *p)
 {
   rabin *r = p;
+  mp *a = MP_NEW;
   int rc = PGEN_ABORT;
 
   switch (rq) {
@@ -126,13 +127,8 @@ int pgen_test(int rq, pgen_event *ev, void *p)
       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);
-      }
+      a = mprand_range(a, ev->m, ev->r, 0);
+      rc = rabin_rtest(r, a);
       break;
     case PGEN_DONE:
       rabin_destroy(r);
@@ -140,6 +136,48 @@ int pgen_test(int rq, pgen_event *ev, void *p)
       break;
   }
 
+  mp_drop(a);
+  return (rc);
+}
+
+/* --- @pgen_bailliepswtest@ --- */
+
+int pgen_bailliepswtest(int rq, pgen_event *ev, void *p)
+{
+  rabin r;
+  int rc;
+
+  switch (rq) {
+    case PGEN_BEGIN:
+      if (ev->tests != 2) rc = PGEN_ABORT;
+      else rc = PGEN_TRY;
+      break;
+
+    case PGEN_DONE:
+      rc = PGEN_DONE;
+      break;
+
+    case PGEN_TRY:
+      switch (ev->tests) {
+       case 2:
+         rabin_create(&r, ev->m);
+         rc = rabin_test(&r, MP_TWO);
+         rabin_destroy(&r);
+         break;
+       case 1:
+         rc = pgen_granfrob(ev->m, 0, 0);
+         break;
+       default:
+         rc = PGEN_ABORT;
+         break;
+      }
+      break;
+
+    default:
+      rc = PGEN_ABORT;
+      break;
+  }
+
   return (rc);
 }
 
@@ -184,8 +222,8 @@ 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 = 0;
-  ev.tests = 0;
+  ev.steps = steps;
+  ev.tests = tests;
   ev.r = fibrand_create(0);
 
   /* --- Tell the event handler we're under way --- */
@@ -261,17 +299,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 == steps) {
+      ev.steps--;
+      if (!ev.steps) {
        act |= A_EVENT | A_ENDSTEP | A_DONE;
        rc = PGEN_ABORT;
       }
-      ev.tests = 0;
+      ev.tests = tests;
     }
 
     if ((act & A_TEST) && tests) {
-      ev.tests++;
-      if (ev.tests == tests) {
+      ev.tests--;
+      if (!ev.tests) {
        act |= A_ENDTEST | A_ENDSTEP | A_DONE;
        rc = PGEN_DONE;
       }
@@ -283,7 +321,7 @@ 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 (p == P_TEST)
+       if (p == P_TEST && rq != PGEN_BEGIN)
          act |= A_ENDTEST;
       }
     }