From 48c2100a8291bf019cd71165de5a58a10c2f8311 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Niels=20M=C3=B6ller?= Date: Tue, 23 Oct 2007 23:50:19 +0200 Subject: [PATCH] Implemented block-wise sieving. Rev: nettle/examples/eratosthenes.c:1.2 --- examples/eratosthenes.c | 130 ++++++++++++++++++++++------------------ 1 file changed, 71 insertions(+), 59 deletions(-) diff --git a/examples/eratosthenes.c b/examples/eratosthenes.c index 82931372..1ab4d223 100644 --- a/examples/eratosthenes.c +++ b/examples/eratosthenes.c @@ -53,12 +53,13 @@ usage(void) fprintf(stderr, "Usage: erathostenes [OPTIONS] [LIMIT]\n\n" "Options:\n" " --help Display this message.\n" + " --quiet No summary line.\n" " --odd-only Omit the prime 2.\n" " --primes-only Suppress output of differences.\n" " --diff-only Supress output of primes.\n" " --tabular Tabular output (default is one prime per line).\n" " --binary Binary output.\n" - " --quiet No summary line.\n"); + " --block SIZE Block size.\n"); } static unsigned @@ -103,58 +104,14 @@ vector_alloc(unsigned size) static void vector_clear_bits (unsigned long *vector, unsigned step, unsigned start, unsigned size) { -#if 0 - if (step <= BITS_PER_LONG/2) - { - unsigned long buf[BITS_PER_LONG - 1]; - - unsigned end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG; - - unsigned i = start / BITS_PER_LONG; - unsigned bit = start % BITS_PER_LONG; - unsigned long mask; - unsigned j; - - if (bit >= step) - { - /* Must handle first word specially */ - for (mask = 0; bit < BITS_PER_LONG; bit += step) - mask |= 1L << (bit % BITS_PER_LONG); - - vector[i++] &= ~mask; - bit -= BITS_PER_LONG; - } - - /* Fill buf with bit pattern */ - for (j = 0; j < step; j++) - { - for (mask = 0; bit < BITS_PER_LONG; bit += step) - mask |= 1L << (bit % BITS_PER_LONG); - - buf[j] = ~mask; - bit -= BITS_PER_LONG; - } - - /* Apply bit pattern */ - while (i + step < end) - for (j = 0; j < step; ) - vector[i++] &= buf[j++]; + unsigned bit; - for (j = 0; i < end; ) - vector[i++] &= buf[j++]; - } - else -#endif + for (bit = start; bit < size; bit += step) { - unsigned bit; + unsigned i = bit / BITS_PER_LONG; + unsigned long mask = 1L << (bit % BITS_PER_LONG); - for (bit = start; bit < size; bit += step) - { - unsigned i = bit / BITS_PER_LONG; - unsigned long mask = 1L << (bit % BITS_PER_LONG); - - vector[i] &= ~mask; - } + vector[i] &= ~mask; } } @@ -369,28 +326,35 @@ main (int argc, char **argv) unsigned bit; unsigned prime_count; unsigned sieve_limit; + unsigned block_size; + unsigned block; struct output_info info; int quiet; int c; - enum { FLAG_ODD = -100, FLAG_PRIME, FLAG_DIFF, FLAG_TABULAR, FLAG_BINARY, FLAG_QUIET }; + enum { FLAG_ODD = -100, FLAG_PRIME, FLAG_DIFF, FLAG_TABULAR, FLAG_BINARY, + FLAG_QUIET, FLAG_BLOCK }; static const struct option options[] = { /* Name, args, flag, val */ { "help", no_argument, NULL, '?' }, + { "quiet", no_argument, NULL, FLAG_QUIET }, { "odd-only", no_argument, NULL, FLAG_ODD }, { "prime-only", no_argument, NULL, FLAG_PRIME }, { "diff-only", no_argument, NULL, FLAG_DIFF }, { "tabular", no_argument, NULL, FLAG_TABULAR }, { "binary", no_argument, NULL, FLAG_BINARY }, - { "quiet" , no_argument, NULL, FLAG_QUIET }, + { "block" , required_argument, NULL, FLAG_BLOCK }, { NULL, 0, NULL, 0} }; output_init(&info); + quiet = 0; + block_size = 0; + while ( (c = getopt_long(argc, argv, "?", options, NULL)) != -1) switch (c) { @@ -415,6 +379,17 @@ main (int argc, char **argv) case FLAG_QUIET: quiet = 1; break; + case FLAG_BLOCK: + { + int arg = atoi(optarg); + if (arg <= 10) + { + usage(); + return EXIT_FAILURE; + } + block_size = (arg - 3) / 2; + break; + } default: abort(); } @@ -438,6 +413,9 @@ main (int argc, char **argv) size = (limit - 1) / 2; + if (!block_size || block_size > size) + block_size = size; + vector = vector_alloc (size); if (!vector) { @@ -453,7 +431,7 @@ main (int argc, char **argv) if (limit == 2) return EXIT_SUCCESS; - sieve_limit = (isqrt(limit) - 1) / 2; + sieve_limit = (isqrt(2*block_size + 1) - 1) / 2; while (bit < sieve_limit) { @@ -466,22 +444,56 @@ main (int argc, char **argv) (n^2 - 3) / 2 = n * bit + 3 (bit + 1) */ - vector_clear_bits (vector, n, n*bit + 3*(bit + 1), size); + vector_clear_bits (vector, n, n*bit + 3*(bit + 1), block_size); - bit = vector_find_next (vector, bit + 1, size); + bit = vector_find_next (vector, bit + 1, block_size); } /* No more marking, just output the remaining primes. */ - while (bit < size) + while (bit < block_size) { - unsigned n = 3 + 2 * bit; - - output(&info, n); + output(&info, 3 + 2 * bit); prime_count++; bit = vector_find_next (vector, bit + 1, size); } + for (block = block_size; block < size; block += block_size) + { + unsigned block_start; + unsigned block_end; + + if (block + block_size > size) + /* For the final block */ + block_size = size - block; + + block_start = 2*block + 3; + block_end = 2*(block + block_size) + 3; + + sieve_limit = (isqrt(block_end) - 1) / 2; + for (bit = 0; bit < sieve_limit ;) + { + unsigned n = 3 + 2 * bit; + + unsigned start = n*bit + 3*(bit + 1); + if (start < block) + { + unsigned k = (block + 1) / n; + start = bit + k*n; + } + vector_clear_bits (vector, n, start, block + block_size); + + bit = vector_find_next (vector, bit + 1, block + block_size); + } + for (bit = vector_find_next (vector, block, block + block_size); + bit < block + block_size; + bit = vector_find_next (vector, bit + 1, block + block_size)) + { + output(&info, 3 + 2 * bit); + prime_count++; + } + } + if (!quiet) { printf("\n"); -- 2.47.3