examples/eratosthenes.c

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