--- /dev/null
+/* -*-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 -------------------------------------------------*/
projects / u / mdw / catacomb / blobdiff