From: David Gibson Date: Thu, 2 Oct 2014 14:15:01 +0000 (+1000) Subject: eratosthenes: Implementation of the Sieve of Eratosthenes X-Git-Url: https://git.ozlabs.org/?p=ccan;a=commitdiff_plain;h=c5e84ef2b5687fc35cdffa4768da13674afb8977 eratosthenes: Implementation of the Sieve of Eratosthenes This adds a new "eratosthenes" module which implements the standard Sieve of Eratosthenes algorithm for locating small prime numbers. Signed-off-by: Rusty Russell --- diff --git a/ccan/eratosthenes/LICENSE b/ccan/eratosthenes/LICENSE new file mode 120000 index 00000000..dc314eca --- /dev/null +++ b/ccan/eratosthenes/LICENSE @@ -0,0 +1 @@ +../../licenses/LGPL-2.1 \ No newline at end of file diff --git a/ccan/eratosthenes/_info b/ccan/eratosthenes/_info new file mode 100644 index 00000000..a09e8d1f --- /dev/null +++ b/ccan/eratosthenes/_info @@ -0,0 +1,45 @@ +#include "config.h" +#include +#include + +/** + * eratosthenes - Sieve of Eratosthenes + * + * This code implements Eratosthenes' Sieve for efficiently finding + * small prime numbers (in this context anything less than several + * billion is "small"). + * + * Example: + * #include + * + * int main(int argc, char *argv[]) + * { + * struct eratosthenes s; + * unsigned long p; + * + * eratosthenes_init(&s); + * eratosthenes_sieve(&s, atol(argv[1])); + * + * while ((p = eratosthenes_nextprime(&s, p)) != 0) { + * printf("%ld\n", p); + * } + * + * return 0; + * } + * + * License: LGPL (v2.1 or any later version) + * Author: David Gibson + */ +int main(int argc, char *argv[]) +{ + /* Expect exactly one argument */ + if (argc != 2) + return 1; + + if (strcmp(argv[1], "depends") == 0) { + printf("ccan/bitmap\n"); + return 0; + } + + return 1; +} diff --git a/ccan/eratosthenes/eratosthenes.c b/ccan/eratosthenes/eratosthenes.c new file mode 100644 index 00000000..7c39c3ad --- /dev/null +++ b/ccan/eratosthenes/eratosthenes.c @@ -0,0 +1,116 @@ +/* Licensed under LGPLv2+ - see LICENSE file for details */ +#include +#include + +#include +#include + +#define VAL_TO_BIT(v) (((v) - 3) / 2) +#define LIMIT_TO_NBITS(l) ((l > 2) ? ((l) - 2) / 2 : 0) + +#define BIT_TO_VAL(b) (((b) * 2) + 3) + +void eratosthenes_init(struct eratosthenes *s) +{ + s->limit = 0; + s->b = NULL; +} + +void eratosthenes_reset(struct eratosthenes *s) +{ + if (s->b) + free(s->b); + eratosthenes_init(s); +} + +static void eratosthenes_once(struct eratosthenes *s, unsigned long limit, unsigned long p) +{ + unsigned long n = VAL_TO_BIT(3*p); + unsigned long obits = LIMIT_TO_NBITS(s->limit); + + if (obits > n) { + n = obits + p - 1 - ((obits - n - 1) % p); + } + + assert((BIT_TO_VAL(n) % p) == 0); + assert((BIT_TO_VAL(n) / p) > 1); + + while (n < LIMIT_TO_NBITS(limit)) { + bitmap_clear_bit(s->b, n); + n += p; + } +} + +static void eratosthenes_sieve_(struct eratosthenes *s, unsigned long limit) +{ + unsigned long p = 3; + + while ((p * p) < limit) { + unsigned long n; + + eratosthenes_once(s, limit, p); + + n = bitmap_ffs(s->b, VAL_TO_BIT(p) + 1, LIMIT_TO_NBITS(limit)); + + /* We should never run out of primes */ + assert(n < LIMIT_TO_NBITS(limit)); + + p = BIT_TO_VAL(n); + } +} + +void eratosthenes_sieve(struct eratosthenes *s, unsigned long limit) +{ + if ((limit < 3) || (limit <= s->limit)) + /* Nothing to do */ + return; + + if (s->limit < 3) + s->b = bitmap_alloc1(LIMIT_TO_NBITS(limit)); + else + s->b = bitmap_realloc1(s->b, LIMIT_TO_NBITS(s->limit), + LIMIT_TO_NBITS(limit)); + + if (!s->b) + abort(); + + eratosthenes_sieve_(s, limit); + + s->limit = limit; +} + +bool eratosthenes_isprime(const struct eratosthenes *s, unsigned long n) +{ + assert(n < s->limit); + + if ((n % 2) == 0) + return (n == 2); + + if (n < 3) { + assert(n == 1); + return false; + } + + return bitmap_test_bit(s->b, VAL_TO_BIT(n)); +} + +unsigned long eratosthenes_nextprime(const struct eratosthenes *s, unsigned long n) +{ + unsigned long i; + + if ((n + 1) >= s->limit) + return 0; + + if (n < 2) + return 2; + + if (n == 2) + return 3; + + i = bitmap_ffs(s->b, VAL_TO_BIT(n) + 1, LIMIT_TO_NBITS(s->limit)); + if (i == LIMIT_TO_NBITS(s->limit)) + /* Reached the end of the sieve */ + return 0; + + return BIT_TO_VAL(i); +} diff --git a/ccan/eratosthenes/eratosthenes.h b/ccan/eratosthenes/eratosthenes.h new file mode 100644 index 00000000..0ec4692d --- /dev/null +++ b/ccan/eratosthenes/eratosthenes.h @@ -0,0 +1,27 @@ +/* Licensed under LGPLv2+ - see LICENSE file for details */ +#ifndef CCAN_ERATOSTHENES_H_ +#define CCAN_ERATOSTHENES_H_ + +#include "config.h" + +#include + +#include + +struct eratosthenes { + unsigned long limit; + bitmap *b; +}; + +void eratosthenes_init(struct eratosthenes *s); + +void eratosthenes_reset(struct eratosthenes *s); + +void eratosthenes_sieve(struct eratosthenes *s, unsigned long limit); + +bool eratosthenes_isprime(const struct eratosthenes *s, unsigned long n); + +unsigned long eratosthenes_nextprime(const struct eratosthenes *s, + unsigned long n); + +#endif /* CCAN_ERATOSTHENES_H_ */ diff --git a/ccan/eratosthenes/test/run-incremental.c b/ccan/eratosthenes/test/run-incremental.c new file mode 100644 index 00000000..ea6b2b9c --- /dev/null +++ b/ccan/eratosthenes/test/run-incremental.c @@ -0,0 +1,36 @@ +#include +#include + +#include + +#define LIMIT 500 + +#define ok_eq(a, b) \ + ok((a) == (b), "%s [%u] == %s [%u]", \ + #a, (unsigned)(a), #b, (unsigned)(b)) + +int main(void) +{ + struct eratosthenes s1, s2; + unsigned long n; + + /* This is how many tests you plan to run */ + plan_tests(LIMIT); + + eratosthenes_init(&s1); + eratosthenes_sieve(&s1, LIMIT); + + eratosthenes_init(&s2); + for (n = 1; n <= LIMIT; n++) + eratosthenes_sieve(&s2, n); + + for (n = 0; n < LIMIT; n++) + ok1(eratosthenes_isprime(&s1, n) + == eratosthenes_isprime(&s2, n)); + + eratosthenes_reset(&s1); + eratosthenes_reset(&s2); + + /* This exits depending on whether all tests passed */ + return exit_status(); +} diff --git a/ccan/eratosthenes/test/run.c b/ccan/eratosthenes/test/run.c new file mode 100644 index 00000000..a40b32ba --- /dev/null +++ b/ccan/eratosthenes/test/run.c @@ -0,0 +1,57 @@ +#include +#include + +#include + +#define LIMIT 500 + +#define ok_eq(a, b) \ + ok((a) == (b), "%s [%u] == %s [%u]", \ + #a, (unsigned)(a), #b, (unsigned)(b)) + +static bool test_isprime(unsigned long n) +{ + int i; + + if (n < 2) + return false; + + for (i = 2; i < n; i++) + if ((n % i) == 0) + return false; + + return true; +} + +static unsigned long test_nextprime(struct eratosthenes *s, unsigned long n) +{ + unsigned long i = n + 1; + + while ((i < LIMIT) && !eratosthenes_isprime(s, i)) + i++; + + return (i >= LIMIT) ? 0 : i; +} + +int main(void) +{ + struct eratosthenes s; + unsigned long n; + + /* This is how many tests you plan to run */ + plan_tests(2 * LIMIT); + + eratosthenes_init(&s); + + eratosthenes_sieve(&s, LIMIT); + + for (n = 0; n < LIMIT; n++) { + ok_eq(eratosthenes_isprime(&s, n), test_isprime(n)); + ok_eq(eratosthenes_nextprime(&s, n), test_nextprime(&s, n)); + } + + eratosthenes_reset(&s); + + /* This exits depending on whether all tests passed */ + return exit_status(); +}