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>

ccan/eratosthenes/LICENSE

@@ -0,0 +1 @@
+../../licenses/LGPL-2.1

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;
+}

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);
+}

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_ */

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();
+}

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();
+}