From: Mark Wooding Date: Tue, 20 Feb 2007 17:22:50 +0000 (+0000) Subject: primeiter: New functions for iterating over small primes. X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/1d6d3b01cb40e0cfaeb816c87c513695aea0816a primeiter: New functions for iterating over small primes. 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. --- diff --git a/Makefile.m4 b/Makefile.m4 index d593d8c..5d5420a 100644 --- a/Makefile.m4 +++ b/Makefile.m4 @@ -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 index 0000000..94d0572 --- /dev/null +++ b/genwheel.c @@ -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 +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +/*----- 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 index 0000000..d7126ba --- /dev/null +++ b/primeiter.c @@ -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 +#include + +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 index 0000000..95cc1b5 --- /dev/null +++ b/primeiter.h @@ -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 diff --git a/tests/pgen b/tests/pgen index de04dbf..6e78dee 100644 --- a/tests/pgen +++ b/tests/pgen @@ -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; +}