pgen: Implement general simultaneous-primality searching.

[u/mdw/catacomb] / pgen-simul.c

diff --git a/pgen-simul.c b/pgen-simul.c
new file mode 100644 (file)
index 0000000..e27f7da
--- /dev/null
+++ b/pgen-simul.c
@@ -0,0 +1,169 @@
+/* -*-c-*-
+ *
+ * $Id$
+ *
+ * Simultaneous prime search
+ *
+ * (c) 2006 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 <stdlib.h>
+
+#include <mLib/macros.h>
+
+#include "mprand.h"
+#include "pgen.h"
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- @rcmerge@ --- *
+ *
+ * Arguments:  @int a, b@ = a pair of @PGEN@ return codes
+ *
+ * Returns:    The overall return code capturing both.
+ */
+
+static int rcmerge(int a, int b)
+{
+#define FROB(state)                                                    \
+    if (a == PGEN_##state || b == PGEN_##state) return (PGEN_##state)
+  FROB(FAIL);
+  FROB(TRY);
+  FROB(PASS);
+#undef FROB
+  return (PGEN_DONE);
+}
+
+/* --- @pgen_simulstep@ --- *
+ *
+ * Step a collection of numbers simultaneously.
+ */
+
+int pgen_simulstep(int rq, pgen_event *ev, void *p)
+{
+  pgen_simulctx *ss = p;
+  pgen_simulprime *sp;
+  int rc = PGEN_ABORT, nrc;
+  unsigned i;
+  mp *q;
+
+  assert(ss->n);
+  switch (rq) {
+    case PGEN_BEGIN:
+      rc = PGEN_DONE;
+      q = MP_NEW;
+      for (i = 0; i < ss->n; i++) {
+       sp = &ss->v[i];
+       q = mp_mul(q, ss->step, sp->mul);
+       if (MP_LEN(q) <= 1)
+         sp->u.step = q->v[0];
+       else {
+         sp->u.jump = xmalloc(sizeof(pfilt));
+         pfilt_create(sp->u.jump, q);
+         sp->f |= PGENF_JUMP;
+       }
+       q = mp_mul(q, ev->m, sp->mul);
+       q = mp_add(q, q, sp->add);
+       nrc = pfilt_create(&sp->p, q);
+       rc = rcmerge(rc, nrc);
+      }
+      MP_DROP(q);
+      if (rc != PGEN_FAIL)
+       goto done;
+    case PGEN_TRY:
+      for (;;) {
+       rc = PGEN_DONE;
+       for (i = 0; i < ss->n; i++) {
+         sp = &ss->v[i];
+         if (sp->f & PGENF_JUMP)
+           nrc = pfilt_jump(&sp->p, sp->u.jump);
+         else
+           nrc = pfilt_step(&sp->p, sp->u.step);
+         rc = rcmerge(rc, nrc);
+       }
+       if (rc != PGEN_FAIL)
+         goto done;
+      }
+    done:
+      mp_drop(ev->m);
+      ev->m = MP_COPY(ss->v[0].p.m);
+      break;
+    case PGEN_DONE:
+      for (i = 0; i < ss->n; i++) {
+       sp = &ss->v[i];
+       if (sp->f & PGENF_JUMP) {
+         pfilt_destroy(sp->u.jump);
+         xfree(sp->u.jump);
+       }
+       if (sp->f & PGENF_KEEP)
+         sp->u.x = MP_COPY(sp->p.m);
+       pfilt_destroy(&sp->p);
+      }
+      rc = PGEN_DONE;
+      break;
+  }
+  return (rc);
+}
+
+/* --- @pgen_simultest@ --- *
+ *
+ * Test a collection of numbers simultaneously.
+ */
+
+int pgen_simultest(int rq, pgen_event *ev, void *p)
+{
+  pgen_simulctx *ss = p;
+  pgen_simulprime *sp;
+  int rc;
+  unsigned i;
+  mp *m;
+
+  assert(ss->n);
+  switch (rq) {
+    case PGEN_BEGIN:
+      for (i = 0; i < ss->n; i++)
+       rabin_create(&ss->v[i].r, ss->v[i].p.m);
+      rc = PGEN_TRY;
+      break;
+    case PGEN_TRY:
+      m = MP_NEW;
+      for (i = 0; i < ss->n; i++) {
+       sp = &ss->v[i];
+       m = mprand_range(m, sp->p.m, ev->r, 0);
+       rc = rabin_test(&sp->r, m);
+       if (rc != PGEN_PASS)
+         break;
+      }
+      mp_drop(m);
+      break;
+    case PGEN_DONE:
+      for (i = 0; i < ss->n; i++)
+       rabin_destroy(&ss->v[i].r);
+      break;
+  }
+  return (rc);
+}
+
+/*----- That's all, folks -------------------------------------------------*/