path: root/nss/lib/freebl/mpi/utils/sieve.c

/*
 * sieve.c
 *
 * Finds prime numbers using the Sieve of Eratosthenes
 *
 * This implementation uses a bitmap to represent all odd integers in a
 * given range.  We iterate over this bitmap, crossing off the
 * multiples of each prime we find.  At the end, all the remaining set
 * bits correspond to prime integers.
 *
 * Here, we make two passes -- once we have generated a sieve-ful of
 * primes, we copy them out, reset the sieve using the highest
 * generated prime from the first pass as a base.  Then we cross out
 * all the multiples of all the primes we found the first time through,
 * and re-sieve.  In this way, we get double use of the memory we
 * allocated for the sieve the first time though.  Since we also
 * implicitly ignore multiples of 2, this amounts to 4 times the
 * values.
 *
 * This could (and probably will) be generalized to re-use the sieve a
 * few more times.
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

typedef unsigned char byte;

typedef struct {
    int size;
    byte *bits;
    long base;
    int next;
    int nbits;
} sieve;

void sieve_init(sieve *sp, long base, int nbits);
void sieve_grow(sieve *sp, int nbits);
long sieve_next(sieve *sp);
void sieve_reset(sieve *sp, long base);
void sieve_cross(sieve *sp, long val);
void sieve_clear(sieve *sp);

#define S_ISSET(S, B) (((S)->bits[(B) / CHAR_BIT] >> ((B) % CHAR_BIT)) & 1)
#define S_SET(S, B) ((S)->bits[(B) / CHAR_BIT] |= (1 << ((B) % CHAR_BIT)))
#define S_CLR(S, B) ((S)->bits[(B) / CHAR_BIT] &= ~(1 << ((B) % CHAR_BIT)))
#define S_VAL(S, B) ((S)->base + (2 * (B)))
#define S_BIT(S, V) (((V) - ((S)->base)) / 2)

int
main(int argc, char *argv[])
{
    sieve s;
    long pr, *p;
    int c, ix, cur = 0;

    if (argc < 2) {
        fprintf(stderr, "Usage: %s <width>\n", argv[0]);
        return 1;
    }

    c = atoi(argv[1]);
    if (c < 0)
        c = -c;

    fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c);

    sieve_init(&s, 3, c);

    c = 0;
    while ((pr = sieve_next(&s)) > 0) {
        ++c;
    }

    p = calloc(c, sizeof(long));
    if (!p) {
        fprintf(stderr, "%s: out of memory after first half\n", argv[0]);
        sieve_clear(&s);
        exit(1);
    }

    fprintf(stderr, "%s: half done ... \n", argv[0]);

    for (ix = 0; ix < s.nbits; ix++) {
        if (S_ISSET(&s, ix)) {
            p[cur] = S_VAL(&s, ix);
            printf("%ld\n", p[cur]);
            ++cur;
        }
    }

    sieve_reset(&s, p[cur - 1]);
    fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur);
    for (ix = 0; ix < cur; ix++) {
        sieve_cross(&s, p[ix]);
        if (!(ix % 1000))
            fputc('.', stderr);
    }
    fputc('\n', stderr);

    free(p);

    fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]);
    c = 0;
    while ((pr = sieve_next(&s)) > 0) {
        ++c;
    }

    fprintf(stderr, "%s: done!\n", argv[0]);
    for (ix = 0; ix < s.nbits; ix++) {
        if (S_ISSET(&s, ix)) {
            printf("%ld\n", S_VAL(&s, ix));
        }
    }

    sieve_clear(&s);

    return 0;
}

void
sieve_init(sieve *sp, long base, int nbits)
{
    sp->size = (nbits / CHAR_BIT);

    if (nbits % CHAR_BIT)
        ++sp->size;

    sp->bits = calloc(sp->size, sizeof(byte));
    memset(sp->bits, UCHAR_MAX, sp->size);
    if (!(base & 1))
        ++base;
    sp->base = base;

    sp->next = 0;
    sp->nbits = sp->size * CHAR_BIT;
}

void
sieve_grow(sieve *sp, int nbits)
{
    int ns = (nbits / CHAR_BIT);

    if (nbits % CHAR_BIT)
        ++ns;

    if (ns > sp->size) {
        byte *tmp;
        int ix;

        tmp = calloc(ns, sizeof(byte));
        if (tmp == NULL) {
            fprintf(stderr, "Error: out of memory in sieve_grow\n");
            return;
        }

        memcpy(tmp, sp->bits, sp->size);
        for (ix = sp->size; ix < ns; ix++) {
            tmp[ix] = UCHAR_MAX;
        }

        free(sp->bits);
        sp->bits = tmp;
        sp->size = ns;

        sp->nbits = sp->size * CHAR_BIT;
    }
}

long
sieve_next(sieve *sp)
{
    long out;
    int ix = 0;
    long val;

    if (sp->next > sp->nbits)
        return -1;

    out = S_VAL(sp, sp->next);
#ifdef DEBUG
    fprintf(stderr, "Sieving %ld\n", out);
#endif

    /* Sieve out all multiples of the current prime */
    val = out;
    while (ix < sp->nbits) {
        val += out;
        ix = S_BIT(sp, val);
        if ((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */
            S_CLR(sp, ix);
#ifdef DEBUG
            fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix);
#endif
        }
    }

    /* Scan ahead to the next prime */
    ++sp->next;
    while (sp->next < sp->nbits && !S_ISSET(sp, sp->next))
        ++sp->next;

    return out;
}

void
sieve_cross(sieve *sp, long val)
{
    int ix = 0;
    long cur = val;

    while (cur < sp->base)
        cur += val;

    ix = S_BIT(sp, cur);
    while (ix < sp->nbits) {
        if (cur & 1)
            S_CLR(sp, ix);
        cur += val;
        ix = S_BIT(sp, cur);
    }
}

void
sieve_reset(sieve *sp, long base)
{
    memset(sp->bits, UCHAR_MAX, sp->size);
    sp->base = base;
    sp->next = 0;
}

void
sieve_clear(sieve *sp)
{
    if (sp->bits)
        free(sp->bits);

    sp->bits = NULL;
}