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., 59 Temple Place - Suite 330, Boston,
  24  * MA 02111-1307, 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   unsigned 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   /* Isolate least significant bit */
 149   x &= -x;
 150
 151   unsigned i = 0;
 152 #if NEED_HANDLE_LARGE_LONG
 153 #ifndef SIZEOF_LONG
 154   /* Can not be tested by the preprocessor. May generate warnings
 155      when long is 32 bits. */
 156   if (BITS_PER_LONG > 32)
 157 #endif
 158     while (x >= 0x100000000L)
 159       {
 160         x >>= 32;
 161         i += 32;
 162       }
 163 #endif /* NEED_HANDLE_LARGE_LONG */
 164
 165   if (x >= 0x10000)
 166     {
 167       x >>= 16;
 168       i += 16;
 169     }
 170   return i + table[128 + (x & 0xff) - (x >> 8)];
 171 }
 172
 173 /* Returns size if there's no more bits set */
 174 static unsigned long
 175 vector_find_next (const unsigned long *vector, unsigned long bit, unsigned long size)
 176 {
 177   unsigned long end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG;
 178   unsigned long i = bit / BITS_PER_LONG;
 179   unsigned long mask = 1L << (bit % BITS_PER_LONG);
 180   unsigned long word;
 181
 182   if (i >= end)
 183     return size;
 184
 185   for (word = vector[i] & ~(mask - 1); !word; word = vector[i])
 186     if (++i >= end)
 187       return size;
 188
 189   /* Next bit is the least significant bit of word */
 190   return i * BITS_PER_LONG + find_first_one(word);
 191 }
 192
 193 /* For benchmarking, define to do nothing (otherwise, most of the time
 194    will be spent converting the output to decimal). */
 195 #define OUTPUT(n) printf("%lu\n", (n))
 196
 197 static long
 198 atosize(const char *s)
 199 {
 200   char *end;
 201   long value = strtol(s, &end, 10);
 202
 203   if (value <= 0)
 204     return 0;
 205
 206   /* FIXME: Doesn't check for overflow. */
 207   switch(*end)
 208     {
 209     default:
 210       return 0;
 211     case '\0':
 212       break;
 213     case 'k': case 'K':
 214       value <<= 10;
 215       break;
 216     case 'M':
 217       value <<= 20;
 218       break;
 219     }
 220   return value;
 221 }
 222
 223 int
 224 main (int argc, char **argv)
 225 {
 226   /* Generate all primes p <= limit */
 227   unsigned long limit;
 228   unsigned long root;
 229
 230   unsigned long limit_nbits;
 231
 232   /* Represents numbers up to sqrt(limit) */
 233   unsigned long sieve_nbits;
 234   unsigned long *sieve;
 235   /* Block for the rest of the sieving. Size should match the cache,
 236      the default value corresponds to 64 KB. */
 237   unsigned long block_nbits = 64L << 13;
 238   unsigned long block_start_bit;
 239   unsigned long *block;
 240
 241   unsigned long bit;
 242   int silent = 0;
 243   int verbose = 0;
 244   int c;
 245
 246   while ( (c = getopt(argc, argv, "?svb:")) != -1)
 247     switch (c)
 248       {
 249       case '?':
 250         usage();
 251         return EXIT_FAILURE;
 252       case 'b':
 253         block_nbits = CHAR_BIT * atosize(optarg);
 254         if (!block_nbits)
 255           {
 256             usage();
 257             return EXIT_FAILURE;
 258           }
 259         break;
 260
 261       case 's':
 262         silent = 1;
 263         break;
 264
 265       case 'v':
 266         verbose++;
 267         break;
 268
 269       default:
 270         abort();
 271       }
 272
 273   argc -= optind;
 274   argv += optind;
 275
 276   if (argc == 0)
 277     limit = 1000;
 278   else if (argc == 1)
 279     {
 280       limit = atol(argv[0]);
 281       if (limit < 2)
 282         return EXIT_SUCCESS;
 283     }
 284   else
 285     {
 286       usage();
 287       return EXIT_FAILURE;
 288     }
 289
 290   root = isqrt(limit);
 291   /* Round down to odd */
 292   root = (root - 1) | 1;
 293   /* Represents odd numbers from 3 up. */
 294   sieve_nbits = (root - 1) / 2;
 295   sieve = vector_alloc(sieve_nbits );
 296   vector_init(sieve, sieve_nbits);
 297
 298   if (verbose)
 299     fprintf(stderr, "Initial sieve using %lu bits.\n", sieve_nbits);
 300
 301   if (!silent)
 302     printf("2\n");
 303
 304   if (limit == 2)
 305     return EXIT_SUCCESS;
 306
 307   for (bit = 0;
 308        bit < sieve_nbits;
 309        bit = vector_find_next(sieve, bit + 1, sieve_nbits))
 310     {
 311       unsigned long n = 3 + 2 * bit;
 312       /* First bit to clear corresponds to n^2, which is bit
 313
 314          (n^2 - 3) / 2 = (n + 3) * bit + 3
 315       */
 316       unsigned long n2_bit = (n+3)*bit + 3;
 317
 318       if (!silent)
 319         printf("%lu\n", n);
 320
 321       vector_clear_bits (sieve, n, n2_bit, sieve_nbits);
 322     }
 323
 324   limit_nbits = (limit - 1) / 2;
 325
 326   if (sieve_nbits + block_nbits > limit_nbits)
 327     block_nbits = limit_nbits - sieve_nbits;
 328
 329   if (verbose)
 330     {
 331       double storage = block_nbits / 8.0;
 332       unsigned shift = 0;
 333       const char prefix[] = " KMG";
 334
 335       while (storage > 1024 && shift < 3)
 336         {
 337           storage /= 1024;
 338           shift++;
 339         }
 340       fprintf(stderr, "Blockwise sieving using blocks of %lu bits (%.3g %cByte)\n",
 341               block_nbits, storage, prefix[shift]);
 342     }
 343
 344   block = vector_alloc(block_nbits);
 345
 346   for (block_start_bit = bit; block_start_bit < limit_nbits; block_start_bit += block_nbits)
 347     {
 348       unsigned long block_start;
 349
 350       if (block_start_bit + block_nbits > limit_nbits)
 351         block_nbits = limit_nbits - block_start_bit;
 352
 353       vector_init(block, block_nbits);
 354
 355       block_start = 3 + 2*block_start_bit;
 356
 357       if (verbose > 1)
 358         fprintf(stderr, "Next block, n = %lu\n", block_start);
 359
 360       /* Sieve */
 361       for (bit = 0; bit < sieve_nbits;
 362            bit = vector_find_next(sieve, bit + 1, sieve_nbits))
 363         {
 364           unsigned long n = 3 + 2 * bit;
 365           unsigned long sieve_start_bit = (n + 3) * bit + 3;
 366
 367           if (sieve_start_bit < block_start_bit)
 368             {
 369               unsigned long k = (block_start + n - 1) / (2*n);
 370               sieve_start_bit = n * k + bit;
 371
 372               assert(sieve_start_bit < block_start_bit + n);
 373             }
 374           assert(sieve_start_bit >= block_start_bit);
 375
 376           vector_clear_bits(block, n, sieve_start_bit - block_start_bit, block_nbits);
 377         }
 378       for (bit = vector_find_next(block, 0, block_nbits);
 379            bit < block_nbits;
 380            bit = vector_find_next(block, bit + 1, block_nbits))
 381         {
 382           unsigned long n = block_start + 2 * bit;
 383           if (!silent)
 384             printf("%lu\n", n);
 385         }
 386     }
 387   return EXIT_SUCCESS;
 388 }