primeiter: New functions for iterating over small primes.
authorMark Wooding <mdw@distorted.org.uk>
Tue, 20 Feb 2007 17:22:50 +0000 (17:22 +0000)
committerMark Wooding <mdw@distorted.org.uk>
Wed, 21 Feb 2007 12:38:03 +0000 (12:38 +0000)
The primeiter functions return consecutive prime numbers from a given
starting point.  To help do this efficiently we use a `wheel': a table
of steps to make which avoid integers with small factors.  The wheel is
generated by a new build-time utility genwheel.c.

Makefile.m4
genwheel.c [new file with mode: 0644]
primeiter.c [new file with mode: 0644]
primeiter.h [new file with mode: 0644]
tests/pgen

index d593d8c..5d5420a 100644 (file)
@@ -90,7 +90,8 @@ _(safer) _(mars) _(tiger) dnl
 _(gfshare) _(gfx-sqr)')
 
 autoheaders: \
-  addsuffix(`gen_tables', `-tab.h') primetab.h mptypes.h mplimits.h
+  addsuffix(`gen_tables', `-tab.h') \
+  primetab.h wheel.h mptypes.h mplimits.h
 define(`emit', `
 _item`'-tab.h: _item`'-mktab$(EXEEXT)
        ./_item`'-mktab >_item`'-tab.h.new
@@ -102,6 +103,13 @@ primetab.c: genprimes$(EXEEXT)
        ./genprimes -h primetab.h -c primetab.c \
                -s CATACOMB_PRIMETAB_H -n 256 \
                -t "unsigned short" -i primetab
+
+wheel.h: wheel.c
+wheel.c: genwheel$(EXEEXT)
+       ./genwheel -h wheel.h -c wheel.c \
+               -s CATACOMB_WHEEL_H -n 5
+               -t "unsigned char" -i wheel
+
 archinclude_HEADERS = mptypes.h mplimits.h
 mptypes.h: mptypes$(EXEEXT)
        ./mptypes >mptypes.h.new
@@ -168,8 +176,8 @@ pkginclude_HEADERS = \
        exp.h mpbarrett.h mpmont.h mpreduce.h mp-exp.h \
        mpcrt.h mprand.h mpmul.h \
        gfx.h gf.h gfreduce.h gfn.h gf-exp.h \
-       primetab.h pfilt.h rabin.h \
-       pgen.h prim.h strongprime.h limlee.h keycheck.h \
+       primetab.h wheel.h pfilt.h rabin.h \
+       pgen.h primeiter.h prim.h strongprime.h limlee.h keycheck.h \
        bbs.h rsa.h dh.h dsarand.h dsa.h gdsa.h gkcdsa.h \
        tlsprf.h sslprf.h \
        gfshare.h share.h \
@@ -214,7 +222,7 @@ define(`EC_SOURCES',
        ec-fetch.c ec-raw.c g-ec.c')
 
 define(`PGEN_SOURCES',
-       `pfilt.c rabin.c \
+       `pfilt.c primeiter.c rabin.c \
        pgen.c pgen-stdev.c pgen-gcd.c pgen-simul.c \
          prim.c strongprime.c limlee.c \
        keycheck.c keycheck-mp.c keycheck-report.c \
@@ -227,7 +235,7 @@ define(`PGEN_SOURCES',
        key-data.c key-flags.c key-text.c key-binary.c key-pass.c \
        key-pack.c key-misc.c key-file.c key-attr.c key-io.c key-moan.c \
        key-error.c key-fetch.c \
-       primetab.c share.c')
+       primetab.c wheel.c share.c')
 
 libcatacomb_la_SOURCES = \
        grand.c keysz.c keysz-conv.c \
@@ -271,7 +279,7 @@ bin_PROGRAMS = \
 noinst_LIBRARIES = libcatcrypt.a
 bin_SCRIPTS = catacomb-config xpixie
 noinst_PROGRAMS = \
-       genprimes mptypes genlimits serpent-check bittest mpdump \
+       genprimes genwheel mptypes genlimits serpent-check bittest mpdump \
        perftest \
        addsuffix(`gen_tables', `-mktab')
 LDADD = libcatcrypt.a libcatacomb.la
@@ -282,7 +290,7 @@ libcatcrypt_a_SOURCES = LIBCAT_SRC getdate.y
 
 patsubst(MP_BASE MP_SOURCES, `\.c\>', `.lo') dsig.o keyutil.o rspit.o \
        patsubst(LIBCAT_SRC, `\.c\>', `.o'): \
-       mptypes.h primetab.h
+       mptypes.h primetab.h wheel.h
 patsubst(MP_SOURCES, `\.c\>', `.lo'): mplimits.h
 
 dsig_SOURCES = dsig.c
@@ -312,10 +320,13 @@ serpent_check_LDADD =
 genprimes_SOURCES = genprimes.c
 genprimes_LDADD =
 
+genwheel_SOURCES = genwheel.c
+genwheel_LDADD =
+
 mptypes_SOURCES = mptypes.c
 mptypes_LDADD =
 
-genlimits_SOURCES = genlimits.c MP_BASE
+genlimits_SOURCES = genlimits.c MP_BASE mptypes.h
 genlimits_LDADD =
 genlimits_CFLAGS = $(AM_CFLAGS)
 
@@ -442,6 +453,7 @@ CTESTRIG(group-test)
 CTESTRIG(gdsa)
 CTESTRIG(gkcdsa)
 CTESTRIG(pgen)
+CTESTRIG(primeiter)
 CTESTRIG(dsa-gen)
 CTESTRIG(dsa-sign)
 CTESTRIG(dsa-verify)
@@ -455,7 +467,8 @@ TESTS = serpent-check bittest testprogs
 
 CLEANFILES = \
        *.t$(EXEEXT) *.to *.kr.old \
-       mptypes.h primetab.c primetab.h ectab.c ptab.c bintab.c \
+       mptypes.h primetab.c primetab.h wheel.c wheel.h \
+       ectab.c ptab.c bintab.c \
        addsuffix(`gen_tables', `-tab.h')
 
 ## --- Makefile building (haha!) ---
diff --git a/genwheel.c b/genwheel.c
new file mode 100644 (file)
index 0000000..94d0572
--- /dev/null
@@ -0,0 +1,194 @@
+/* -*-c-*-
+ *
+ * Generate a prime-iteration wheel
+ *
+ * (c) 2007 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 <ctype.h>
+#include <errno.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <mLib/darray.h>
+#include <mLib/dstr.h>
+#include <mLib/mdwopt.h>
+#include <mLib/quis.h>
+#include <mLib/report.h>
+
+/*----- Data structures ---------------------------------------------------*/
+
+DA_DECL(uintv, unsigned int);
+
+/*----- Main code ---------------------------------------------------------*/
+
+static unsigned long gcd(unsigned long a, unsigned long b)
+{
+  int t;
+  if (!a) return (b);
+  while (b) { t = a%b; a = b; b = t; }
+  return (a);
+}
+
+int main(int argc, char *argv[])
+{
+  int np = 5;
+  const char *type = "unsigned char";
+  const char *source = "wheel.c";
+  const char *header = "wheel.h";
+  const char *name = "wheel";
+  const char *sym = 0;
+  unsigned long i, n;
+  unsigned long mod;
+  int o;
+  uintv v = DA_INIT;
+
+  ego(argv[0]);
+  for (;;) {
+    o = getopt(argc, argv, "n:c:h:s:t:i:");
+    if (o < 0)
+      break;
+    switch (o) {
+      case 'n':
+       np = atoi(optarg);
+       break;
+      case 's':
+       sym = optarg;
+       break;
+      case 'c':
+       source = optarg;
+       break;
+      case 'h':
+       header = optarg;
+       break;
+      case 't':
+       type = optarg;
+       break;
+      case 'i':
+       name = optarg;
+       break;
+      default:
+       pquis(stderr, "Usage: $ [-n nprimes] [-s source] [-h header]\n");
+       exit(EXIT_FAILURE);
+    }
+  }
+
+  for (mod = 1, i = 2, n = 0;
+       n < np;
+       i++) {
+    if (gcd(i, mod) == 1) {
+      mod *= i;
+      n++;
+    }
+  }
+
+  n = 1;
+  for (i = 2; i < mod; i++) {
+    if (gcd(mod, i) == 1) {
+      DA_PUSH(&v, i - n);
+      n = i;
+    }
+  }
+  DA_PUSH(&v, mod + 1 - n);
+
+  {
+    FILE *fp = fopen(header, "w");
+    dstr d = DSTR_INIT;
+    const char *q;
+    if (!fp)
+      die(EXIT_FAILURE, "couldn't write `%s': %s", header, strerror(errno));
+    if (!sym) {
+      for (q = header; *q; q++) {
+       int ch = (unsigned char)*q;
+       if (isalnum(ch))
+         ch = toupper(ch);
+       else
+         ch = '_';
+       DPUTC(&d, ch);
+      }
+      DPUTZ(&d);
+      sym = d.buf;
+    }
+    fprintf(fp, "\
+/* -*-c-*-\n\
+ *\n\
+ * Wheel for small prime iteration [generated]\n\
+ */\n\
+\n\
+#ifndef %s\n\
+#define %s\n\
+\n\
+#define WHEELN %luu\n\
+#define WHEELMOD %luu\n\
+\n\
+extern const %s %s[];\n\
+\n\
+#endif\n\
+",
+           sym, sym,
+           (unsigned long)DA_LEN(&v),
+           mod,
+           type, name);
+    dstr_destroy(&d);
+    if (fclose(fp) == EOF) {
+      remove(header);
+      die(EXIT_FAILURE, "error writing `%s': %s", header, strerror(errno));
+    }
+  }
+
+  {
+    FILE *fp = fopen(source, "w");
+    int i;
+    if (!fp)
+      die(EXIT_FAILURE, "couldn't write `%s': %s", source, strerror(errno));
+    fprintf(fp, "\
+/* -*-c-*-\n\
+ *\n\
+ * Wheel for small prime iteration [generated]\n\
+ */\n\
+\n\
+#include \"%s\"\n\
+\n\
+const %s %s[] = {",
+           header, type, name);
+    for (i = 0; i < DA_LEN(&v); i++) {
+      if (i % 8 == 0)
+       fputs("\n  ", fp);
+      fprintf(fp, "%5u, ", DA(&v)[i]);
+    }
+    fputs("\n\
+};\n\
+", fp);
+    if (fclose(fp) == EOF) {
+      remove(source);
+      die(EXIT_FAILURE, "error writing `%s': %s", source, strerror(errno));
+    }
+  }
+
+  return (0);
+}
+
+/*----- That's all, folks -------------------------------------------------*/
diff --git a/primeiter.c b/primeiter.c
new file mode 100644 (file)
index 0000000..d7126ba
--- /dev/null
@@ -0,0 +1,254 @@
+/* -*-c-*-
+ *
+ * Iterate over small primes efficiently
+ *
+ * (c) 2007 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 "fibrand.h"
+#include "mp.h"
+#include "pgen.h"
+#include "primeiter.h"
+#include "primetab.h"
+#include "wheel.h"
+
+/*----- Theory ------------------------------------------------------------*
+ *
+ * For small primes, we can just pluck them out of the small primes table.
+ * For larger primes, we can test them individually, or build a sieve or
+ * something, but since we don't know when to stop, that could be tricky.
+ *
+ * We've built a `wheel', as follows.  Let %$m$% be the product of the first
+ * %$n$% primes.  There are %$\phi(m)$% integers %$n_i$%, with %$0 < n_i <
+ * m$% coprime to %$m$%, and any integer %$j > n$% must be congruent to some
+ * %$n_i$% modulo %$m$%.  The wheel itself doesn't list the %$n_i$%, but
+ * rather the differences %$\delta_i = n_i - n_{i-1}$% (wrapping around
+ * appropriately at the ends), so you can just add simple offsets to step
+ * onwards.  The wheel assumes you start at 1 and move on round.
+ */
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- @wheelsync@ --- *
+ *
+ * Arguments:  @primeiter *pi@ = iterator to synchronize
+ *             @mp *where@ = value to synchronize
+ *
+ * Returns:    ---
+ *
+ * Use:                Sets up the wheel index to match the given integer.  After
+ *             this, we can step along the wheel to find candidate primes.
+ */
+
+static void wheelsync(primeiter *pi, mp *where)
+{
+  mpw w;
+  mp t;
+  mpw rr;
+  mp *r = MP_NEW;
+  unsigned i, n;
+
+  w = WHEELMOD;
+  mp_build(&t, &w, &w + 1);
+  mp_div(0, &r, where, &t);
+  rr = MP_ZEROP(r) ? 0 : r->v[0];
+
+  for (i = 0, n = 1; rr > n; n += wheel[i], i++);
+  w = n - rr;
+  pi->p = mp_add(MP_NEW, where, &t);
+  pi->i = i;
+  pi->r = fibrand_create(0);
+  MP_DROP(r);
+}
+
+/* --- @primeiter_create@ --- *
+ *
+ * Arguments:  @primeiter *pi@ = pointer to an iterator structure
+ *             @mp *start@ = where to start
+ *
+ * Returns:    ---
+ *
+ * Use:                Initializes a prime iterator.  The first output will be the
+ *             smallest prime not less than @start@.
+ */
+
+void primeiter_create(primeiter *pi, mp *start)
+{
+  mpw n;
+  unsigned l, h, m;
+
+  if (!start || MP_CMP(start, <=, MP_TWO))
+    start = MP_TWO;
+
+  if (MP_LEN(start) <= 1) {
+    n = start->v[0];
+    if (n <= MAXPRIME) {
+      l = 0;
+      h = NPRIME;
+      for (;;) {
+       m = l + (h - l)/2;
+       if (primetab[m] == n) break;
+       else if (m == l) { m++; break; }
+       else if (primetab[m] < n) l = m;
+       else h = m;
+      }
+      pi->i = m;
+      pi->mode = PIMODE_PTAB;
+      mp_build(&pi->pp, &pi->w, &pi->w + 1);
+      pi->p = &pi->pp;
+      return;
+    }
+  }
+
+  wheelsync(pi, start);
+  pi->mode = PIMODE_STALL;
+}
+
+/* --- @primeiter_destroy@ --- *
+ *
+ * Arguments:  @primeiter *pi@ = pointer to iterator structure
+ *
+ * Returns:    ---
+ *
+ * Use:                Frees up an iterator structure when it's no longer wanted.
+ */
+
+void primeiter_destroy(primeiter *pi)
+{
+  switch (pi->mode) {
+    case PIMODE_PTAB:
+      break;
+    case PIMODE_STALL:
+    case PIMODE_WHEEL:
+      MP_DROP(pi->p);
+      GR_DESTROY(pi->r);
+      break;
+    default:
+      abort();
+  }
+}
+
+/* --- @primeiter_next@ --- *
+ *
+ * Arguments:  @primeiter *pi@ = pointer to an iterator structure
+ *             @mp *d@ = fake destination
+ *
+ * Returns:    The next prime number from the iterator.
+ *
+ * Use:                Returns a new prime number.
+ */
+
+mp *primeiter_next(primeiter *pi, mp *d)
+{
+  mp *p;
+
+  switch (pi->mode) {
+    case PIMODE_PTAB:
+      pi->w = primetab[pi->i++];
+      if (pi->i >= NPRIME) {
+       wheelsync(pi, pi->p);
+       pi->mode = PIMODE_WHEEL;
+      }
+      p = MP_COPY(pi->p);
+      MP_SPLIT(p);
+      break;
+    case PIMODE_STALL:
+      pi->mode = PIMODE_WHEEL;
+      goto loop;
+    case PIMODE_WHEEL:
+      do {
+       MP_DEST(pi->p, MP_LEN(pi->p) + 1, pi->p->f);
+       MPX_UADDN(pi->p->v, pi->p->vl, wheel[pi->i++]);
+       MP_SHRINK(pi->p);
+       if (pi->i >= WHEELN) pi->i = 0;
+      loop:;
+      } while (!pgen_primep(pi->p, pi->r));
+      p = MP_COPY(pi->p);
+      break;
+    default:
+      abort();
+  }
+  if (d) MP_DROP(d);
+  return (p);
+}
+
+/*----- Test rig ----------------------------------------------------------*/
+
+#ifdef TEST_RIG
+
+#include <mLib/macros.h>
+#include <mLib/testrig.h>
+
+static int test(dstr *v)
+{
+  mp *start = *(mp **)v[0].buf;
+  mp *pp[5], *ret[5];
+  int i;
+  primeiter pi;
+  int ok = 1;
+
+  for (i = 0; i < N(pp); i++)
+    pp[i] = *(mp **)v[i + 1].buf;
+  primeiter_create(&pi, start);
+  for (i = 0; i < N(pp); i++) {
+    ret[i] = primeiter_next(&pi, MP_NEW);
+    if (!MP_EQ(ret[i], pp[i])) ok = 0;
+  }
+  if (!ok) {
+    fprintf(stderr, "\n*** primeiter test failure:\n***   start = ");
+    mp_writefile(start, stderr, 10);
+    for (i = 0; i < N(pp); i++) {
+      fprintf(stderr, "\n***   p[%d] = ", i);
+      mp_writefile(ret[i], stderr, 10);
+      fprintf(stderr, " %s ", MP_EQ(ret[i], pp[i]) ? "==" : "!=");
+      mp_writefile(pp[i], stderr, 10);
+    }
+    fputc('\n', stderr);
+  }
+  for (i = 0; i < N(pp); i++) {
+    MP_DROP(pp[i]);
+    MP_DROP(ret[i]);
+  }
+  primeiter_destroy(&pi);
+  MP_DROP(start);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
+  return (ok);
+}
+
+static test_chunk tests[] = {
+  { "primeiter", test,
+    { &type_mp,  &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, } },
+  { 0 }
+};
+
+int main(int argc, char *argv[])
+{
+  test_run(argc, argv, tests, SRCDIR "/tests/pgen");
+  return (0);
+}
+
+#endif
+
+/*----- That's all, folks -------------------------------------------------*/
diff --git a/primeiter.h b/primeiter.h
new file mode 100644 (file)
index 0000000..95cc1b5
--- /dev/null
@@ -0,0 +1,106 @@
+/* -*-c-*-
+ *
+ * Iterate over small primes efficiently
+ *
+ * (c) 2007 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.
+ */
+
+#ifndef CATACOMB_PRIMEITER_H
+#define CATACOMB_PRIMEITER_H
+
+#ifdef __cplusplus
+  extern "C" {
+#endif
+
+/*----- Header files ------------------------------------------------------*/
+
+#ifndef CATACOMB_GRAND_H
+#  include "grand.h"
+#endif
+
+#ifndef CATACOMB_MP_H
+#  include "mp.h"
+#endif
+
+/*----- Data structures ---------------------------------------------------*/
+
+typedef struct primeiter {
+  mpw w;
+  mp pp;
+  mp *p;
+  grand *r;
+  unsigned mode;
+  int i;
+} primeiter;
+
+enum {
+  PIMODE_PTAB,
+  PIMODE_STALL,
+  PIMODE_WHEEL
+};
+
+/*----- Functions provided ------------------------------------------------*/
+
+/* --- @primeiter_create@ --- *
+ *
+ * Arguments:  @primeiter *pi@ = pointer to an iterator structure
+ *             @mp *start@ = where to start
+ *
+ * Returns:    ---
+ *
+ * Use:                Initializes a prime iterator.  The first output will be the
+ *             smallest prime not less than @start@.
+ */
+
+extern void primeiter_create(primeiter */*pi*/, mp */*start*/);
+
+/* --- @primeiter_destroy@ --- *
+ *
+ * Arguments:  @primeiter *pi@ = pointer to iterator structure
+ *
+ * Returns:    ---
+ *
+ * Use:                Frees up an iterator structure when it's no longer wanted.
+ */
+
+void primeiter_destroy(primeiter */*pi*/);
+
+/* --- @primeiter_next@ --- *
+ *
+ * Arguments:  @primeiter *pi@ = pointer to an iterator structure
+ *             @mp *d@ = fake destination
+ *
+ * Returns:    The next prime number from the iterator.
+ *
+ * Use:                Returns a new prime number.
+ */
+
+mp *primeiter_next(primeiter */*pi*/, mp */*d*/);
+
+/*----- That's all, folks -------------------------------------------------*/
+
+#ifdef __cplusplus
+  }
+#endif
+
+#endif
index de04dbf..6e78dee 100644 (file)
@@ -27,3 +27,14 @@ primep {
   40301809 1;
   40301811 0;
 }
+
+primeiter {
+  0 2 3 5 7 11;
+  2 2 3 5 7 11;
+  3 3 5 7 11 13;
+  4 5 7 11 13 17;
+
+  2309 2309 2311 2333 2339 2341;
+  7878 7879 7883 7901 7907 7919;
+  7879 7879 7883 7901 7907 7919;
+}