examples/eratosthenes.c

   1 /* eratosthenes.c
   2  *
   3  * An implementation of the sieve of Eratosthenes, to generate a list of primes.
   4  *
   5  */
   6
   7 /* nettle, low-level cryptographics library
   8  *
   9  * Copyright (C) 2007 Niels Möller
  10  *
  11  * The nettle library is free software; you can redistribute it and/or modify
  12  * it under the terms of the GNU Lesser General Public License as published by
  13  * the Free Software Foundation; either version 2.1 of the License, or (at your
  14  * option) any later version.
  15  *
  16  * The nettle library is distributed in the hope that it will be useful, but
  17  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  19  * License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public License
  22  * along with the nettle library; see the file COPYING.LIB.  If not, write to
  23  * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  24  * MA 02111-1301, USA.
  25  */
  26
  27 #if HAVE_CONFIG_H
  28 # include "config.h"
  29 #endif
  30
  31 #include <assert.h>
  32 #include <limits.h>
  33 #include <stdio.h>
  34 #include <stdlib.h>
  35
  36 #include "getopt.h"
  37
  38 #ifdef SIZEOF_LONG
  39 # define BITS_PER_LONG (CHAR_BIT * SIZEOF_LONG)
  40 # if BITS_PER_LONG > 32
  41 #  define NEED_HANDLE_LARGE_LONG 1
  42 # else
  43 #  define NEED_HANDLE_LARGE_LONG 0
  44 # endif
  45 #else
  46 # define BITS_PER_LONG (CHAR_BIT * sizeof(unsigned long))
  47 # define NEED_HANDLE_LARGE_LONG 1
  48 #endif
  49
  50
  51 static void
  52 usage(void)
  53 {
  54   fprintf(stderr, "Usage: erathostenes [OPTIONS] [LIMIT]\n\n"
  55           "Options:\n"
  56           "      -?         Display this message.\n"
  57           "      -b SIZE    Block size.\n"
  58           "      -v         Verbose output.\n"
  59           "      -s         No output.\n");
  60 }
  61
  62 static unsigned
  63 isqrt(unsigned long n)
  64 {
  65   unsigned long x;
  66
  67   /* FIXME: Better initialization. */
  68   if (n < ULONG_MAX)
  69     x = n;
  70   else
  71     /* Must avoid overflow in the first step. */
  72     x = n-1;
  73
  74   for (;;)
  75     {
  76       unsigned long y = (x + n/x) / 2;
  77       if (y >= x)
  78         return x;
  79
  80       x = y;
  81     }
  82 }
  83
  84 /* Size is in bits */
  85 static unsigned long *
  86 vector_alloc(unsigned long size)
  87 {
  88   unsigned long end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG;
  89   unsigned long *vector = malloc (end * sizeof(long));
  90
  91   if (!vector)
  92     {
  93       fprintf(stderr, "Insufficient memory.\n");
  94       exit(EXIT_FAILURE);
  95     }
  96   return vector;
  97 }
  98
  99 static void
 100 vector_init(unsigned long *vector, unsigned long size)
 101 {
 102   unsigned long end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG;
 103   unsigned long i;
 104
 105   for (i = 0; i < end; i++)
 106     vector[i] = ~0;
 107 }
 108
 109 static void
 110 vector_clear_bits (unsigned long *vector, unsigned long step,
 111                    unsigned long start, unsigned long size)
 112 {
 113   unsigned long bit;
 114
 115   for (bit = start; bit < size; bit += step)
 116     {
 117       unsigned long i = bit / BITS_PER_LONG;
 118       unsigned long mask = 1L << (bit % BITS_PER_LONG);
 119
 120       vector[i] &= ~mask;
 121     }
 122 }
 123
 124 static unsigned
 125 find_first_one (unsigned long x)
 126 {
 127   static const unsigned char table[0x101] =
 128     {
 129      15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 130       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 131       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 132       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 133      14, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 134       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 135      13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 136      12, 0, 0, 0, 0, 0, 0, 0,11, 0, 0, 0,10, 0, 9, 8,
 137       0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0,
 138       4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 139       5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 140       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 141       6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 142       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 143       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 144       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 145       7,
 146     };
 147
 148   unsigned i = 0;
 149
 150   /* Isolate least significant bit */
 151   x &= -x;
 152
 153 #if NEED_HANDLE_LARGE_LONG
 154 #ifndef SIZEOF_LONG
 155   /* Can not be tested by the preprocessor. May generate warnings
 156      when long is 32 bits. */
 157   if (BITS_PER_LONG > 32)
 158 #endif
 159     while (x >= 0x100000000L)
 160       {
 161         x >>= 32;
 162         i += 32;
 163       }
 164 #endif /* NEED_HANDLE_LARGE_LONG */
 165
 166   if (x >= 0x10000)
 167     {
 168       x >>= 16;
 169       i += 16;
 170     }
 171   return i + table[128 + (x & 0xff) - (x >> 8)];
 172 }
 173
 174 /* Returns size if there's no more bits set */
 175 static unsigned long
 176 vector_find_next (const unsigned long *vector, unsigned long bit, unsigned long size)
 177 {
 178   unsigned long end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG;
 179   unsigned long i = bit / BITS_PER_LONG;
 180   unsigned long mask = 1L << (bit % BITS_PER_LONG);
 181   unsigned long word;
 182
 183   if (i >= end)
 184     return size;
 185
 186   for (word = vector[i] & ~(mask - 1); !word; word = vector[i])
 187     if (++i >= end)
 188       return size;
 189
 190   /* Next bit is the least significant bit of word */
 191   return i * BITS_PER_LONG + find_first_one(word);
 192 }
 193
 194 /* For benchmarking, define to do nothing (otherwise, most of the time
 195    will be spent converting the output to decimal). */
 196 #define OUTPUT(n) printf("%lu\n", (n))
 197
 198 static long
 199 atosize(const char *s)
 200 {
 201   char *end;
 202   long value = strtol(s, &end, 10);
 203
 204   if (value <= 0)
 205     return 0;
 206
 207   /* FIXME: Doesn't check for overflow. */
 208   switch(*end)
 209     {
 210     default:
 211       return 0;
 212     case '\0':
 213       break;
 214     case 'k': case 'K':
 215       value <<= 10;
 216       break;
 217     case 'M':
 218       value <<= 20;
 219       break;
 220     }
 221   return value;
 222 }
 223
 224 int
 225 main (int argc, char **argv)
 226 {
 227   /* Generate all primes p <= limit */
 228   unsigned long limit;
 229   unsigned long root;
 230
 231   unsigned long limit_nbits;
 232
 233   /* Represents numbers up to sqrt(limit) */
 234   unsigned long sieve_nbits;
 235   unsigned long *sieve;
 236   /* Block for the rest of the sieving. Size should match the cache,
 237      the default value corresponds to 64 KB. */
 238   unsigned long block_nbits = 64L << 13;
 239   unsigned long block_start_bit;
 240   unsigned long *block;
 241
 242   unsigned long bit;
 243   int silent = 0;
 244   int verbose = 0;
 245   int c;
 246
 247   enum { OPT_HELP = 300 };
 248   static const struct option options[] =
 249     {
 250       /* Name, args, flag, val */
 251       { "help", no_argument, NULL, OPT_HELP },
 252       { "verbose", no_argument, NULL, 'v' },
 253       { "block-size", required_argument, NULL, 'b' },
 254       { "quiet", required_argument, NULL, 'q' },
 255       { NULL, 0, NULL, 0}
 256     };
 257
 258   while ( (c = getopt_long(argc, argv, "svb:", options, NULL)) != -1)
 259     switch (c)
 260       {
 261       case OPT_HELP:
 262         usage();
 263         return EXIT_SUCCESS;
 264       case 'b':
 265         block_nbits = CHAR_BIT * atosize(optarg);
 266         if (!block_nbits)
 267           {
 268             usage();
 269             return EXIT_FAILURE;
 270           }
 271         break;
 272
 273       case 'q':
 274         silent = 1;
 275         break;
 276
 277       case 'v':
 278         verbose++;
 279         break;
 280
 281       case '?':
 282         return EXIT_FAILURE;
 283
 284       default:
 285         abort();
 286       }
 287
 288   argc -= optind;
 289   argv += optind;
 290
 291   if (argc == 0)
 292     limit = 1000;
 293   else if (argc == 1)
 294     {
 295       limit = atol(argv[0]);
 296       if (limit < 2)
 297         return EXIT_SUCCESS;
 298     }
 299   else
 300     {
 301       usage();
 302       return EXIT_FAILURE;
 303     }
 304
 305   root = isqrt(limit);
 306   /* Round down to odd */
 307   root = (root - 1) | 1;
 308   /* Represents odd numbers from 3 up. */
 309   sieve_nbits = (root - 1) / 2;
 310   sieve = vector_alloc(sieve_nbits );
 311   vector_init(sieve, sieve_nbits);
 312
 313   if (verbose)
 314     fprintf(stderr, "Initial sieve using %lu bits.\n", sieve_nbits);
 315
 316   if (!silent)
 317     printf("2\n");
 318
 319   if (limit == 2)
 320     return EXIT_SUCCESS;
 321
 322   for (bit = 0;
 323        bit < sieve_nbits;
 324        bit = vector_find_next(sieve, bit + 1, sieve_nbits))
 325     {
 326       unsigned long n = 3 + 2 * bit;
 327       /* First bit to clear corresponds to n^2, which is bit
 328
 329          (n^2 - 3) / 2 = (n + 3) * bit + 3
 330       */
 331       unsigned long n2_bit = (n+3)*bit + 3;
 332
 333       if (!silent)
 334         printf("%lu\n", n);
 335
 336       vector_clear_bits (sieve, n, n2_bit, sieve_nbits);
 337     }
 338
 339   limit_nbits = (limit - 1) / 2;
 340
 341   if (sieve_nbits + block_nbits > limit_nbits)
 342     block_nbits = limit_nbits - sieve_nbits;
 343
 344   if (verbose)
 345     {
 346       double storage = block_nbits / 8.0;
 347       unsigned shift = 0;
 348       const char prefix[] = " KMG";
 349
 350       while (storage > 1024 && shift < 3)
 351         {
 352           storage /= 1024;
 353           shift++;
 354         }
 355       fprintf(stderr, "Blockwise sieving using blocks of %lu bits (%.3g %cByte)\n",
 356               block_nbits, storage, prefix[shift]);
 357     }
 358
 359   block = vector_alloc(block_nbits);
 360
 361   for (block_start_bit = bit; block_start_bit < limit_nbits; block_start_bit += block_nbits)
 362     {
 363       unsigned long block_start;
 364
 365       if (block_start_bit + block_nbits > limit_nbits)
 366         block_nbits = limit_nbits - block_start_bit;
 367
 368       vector_init(block, block_nbits);
 369
 370       block_start = 3 + 2*block_start_bit;
 371
 372       if (verbose > 1)
 373         fprintf(stderr, "Next block, n = %lu\n", block_start);
 374
 375       /* Sieve */
 376       for (bit = 0; bit < sieve_nbits;
 377            bit = vector_find_next(sieve, bit + 1, sieve_nbits))
 378         {
 379           unsigned long n = 3 + 2 * bit;
 380           unsigned long sieve_start_bit = (n + 3) * bit + 3;
 381
 382           if (sieve_start_bit < block_start_bit)
 383             {
 384               unsigned long k = (block_start + n - 1) / (2*n);
 385               sieve_start_bit = n * k + bit;
 386
 387               assert(sieve_start_bit < block_start_bit + n);
 388             }
 389           assert(sieve_start_bit >= block_start_bit);
 390
 391           vector_clear_bits(block, n, sieve_start_bit - block_start_bit, block_nbits);
 392         }
 393       for (bit = vector_find_next(block, 0, block_nbits);
 394            bit < block_nbits;
 395            bit = vector_find_next(block, bit + 1, block_nbits))
 396         {
 397           unsigned long n = block_start + 2 * bit;
 398           if (!silent)
 399             printf("%lu\n", n);
 400         }
 401     }
 402
 403   free(sieve);
 404   free(block);
 405
 406   return EXIT_SUCCESS;
 407 }