+++ /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