Rearrange the file tree.

[u/mdw/catacomb] / math / pgen-gcd.c

diff --git a/math/pgen-gcd.c b/math/pgen-gcd.c
new file mode 100644 (file)
index 0000000..1d205b5
--- /dev/null
+++ b/math/pgen-gcd.c
@@ -0,0 +1,101 @@
+/* -*-c-*-
+ *
+ * Prime search stepper ensuring a low GCD for %$(p - 1)/2$%
+ *
+ * (c) 2000 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 "mp.h"
+#include "pgen.h"
+
+/*----- Main code ---------------------------------------------------------*/
+
+int pgen_gcdstep(int rq, pgen_event *ev, void *p)
+{
+  pgen_gcdstepctx *g = p;
+  int rc = PGEN_ABORT;
+
+  switch (rq) {
+
+    /* --- Set everything up --- *
+     *
+     * Call things off if @p@ and @jp@ have common factors, or if @q@, @r@
+     * and @jq@ have common factors greater than @max@.
+     */
+
+    case PGEN_BEGIN: {
+      mp *p = ev->m;
+      mp_gcd(&g->g, 0, 0, p, g->jp.m);
+      if (MP_CMP(g->g, >, MP_ONE))
+       return (PGEN_ABORT);
+      g->q = mp_lsr(MP_NEW, p, 1);
+      g->jq = mp_lsr(MP_NEW, g->jp.m, 1);
+      mp_gcd(&g->g, 0, 0, g->q, g->jq);
+      mp_gcd(&g->g, 0, 0, g->g, g->r);
+      if (MP_CMP(g->g, >, g->max)) {
+       mp_drop(g->q);
+       mp_drop(g->jq);
+       return (PGEN_ABORT);
+      }
+      rc = pfilt_create(&g->p, p);
+      mp_drop(p);
+    } break;
+
+    /* --- Grind through another iteration --- */
+
+    case PGEN_TRY:
+      mp_drop(ev->m);
+      rc = pfilt_jump(&g->p, &g->jp);
+      g->q = mp_add(g->q, g->q, g->jq);
+      break;
+
+    /* --- Finished --- */
+
+    case PGEN_DONE:
+      pfilt_destroy(&g->p);
+      mp_drop(g->q);
+      mp_drop(g->jq);
+      return (PGEN_DONE);
+  }
+
+  /* --- Step on until everything is OK --- */
+
+  for (;;) {
+    if (rc != PGEN_FAIL) {
+      mp_gcd(&g->g, 0, 0, g->r, g->q);
+      if (MP_CMP(g->g, >, g->max))
+       rc = PGEN_FAIL;
+    }
+    if (rc != PGEN_FAIL)
+      break;
+    rc = pfilt_jump(&g->p, &g->jp);
+    g->q = mp_add(g->q, g->q, g->jq);
+  }
+
+  ev->m = MP_COPY(g->p.m);
+  return (rc);
+}
+
+/*----- That's all, folks -------------------------------------------------*/