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 <rusty@rustcorp.com.au>

diff --git a/ccan/eratosthenes/LICENSE b/ccan/eratosthenes/LICENSE
new file mode 120000 (symlink)
index 0000000..dc314ec
--- /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 (file)
index 0000000..a09e8d1
--- /dev/null
+++ b/ccan/eratosthenes/_info
@@ -0,0 +1,45 @@
+#include "config.h"
+#include <stdio.h>
+#include <string.h>
+
+/**
+ * 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 <ccan/eratosthenes/eratosthenes.h>
+ *
+ *     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 <david@gibsond.dropbear.id.au>
+ */
+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 (file)
index 0000000..7c39c3a
--- /dev/null
+++ b/ccan/eratosthenes/eratosthenes.c
@@ -0,0 +1,116 @@
+/* Licensed under LGPLv2+ - see LICENSE file for details */
+#include <ccan/eratosthenes/eratosthenes.h>
+#include <ccan/bitmap/bitmap.h>
+
+#include <assert.h>
+#include <stdlib.h>
+
+#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 (file)
index 0000000..0ec4692
--- /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 <stdbool.h>
+
+#include <ccan/bitmap/bitmap.h>
+
+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 (file)
index 0000000..ea6b2b9
--- /dev/null
+++ b/ccan/eratosthenes/test/run-incremental.c
@@ -0,0 +1,36 @@
+#include <ccan/eratosthenes/eratosthenes.h>
+#include <ccan/tap/tap.h>
+
+#include <ccan/eratosthenes/eratosthenes.c>
+
+#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 (file)
index 0000000..a40b32b
--- /dev/null
+++ b/ccan/eratosthenes/test/run.c
@@ -0,0 +1,57 @@
+#include <ccan/eratosthenes/eratosthenes.h>
+#include <ccan/tap/tap.h>
+
+#include <ccan/eratosthenes/eratosthenes.c>
+
+#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();
+}