eratosthenes.c 9.24 KB
Newer Older
Niels Möller's avatar
Niels Möller committed
1 2 3 4 5
/* eratosthenes.c
 *
 * An implementation of the sieve of Eratosthenes, to generate a list of primes.
 *
 */
Niels Möller's avatar
Niels Möller committed
6

Niels Möller's avatar
Niels Möller committed
7 8 9
/* nettle, low-level cryptographics library
 *
 * Copyright (C) 2007 Niels Mller
Niels Möller's avatar
Niels Möller committed
10
 *
Niels Möller's avatar
Niels Möller committed
11 12 13 14
 * The nettle library is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
Niels Möller's avatar
Niels Möller committed
15
 *
Niels Möller's avatar
Niels Möller committed
16 17 18 19
 * The nettle library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
Niels Möller's avatar
Niels Möller committed
20
 *
Niels Möller's avatar
Niels Möller committed
21 22 23 24 25 26 27 28 29 30
 * You should have received a copy of the GNU Lesser General Public License
 * along with the nettle library; see the file COPYING.LIB.  If not, write to
 * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 * MA 02111-1307, USA.
 */

#if HAVE_CONFIG_H
# include "config.h"
#endif

31
#include <assert.h>
Niels Möller's avatar
Niels Möller committed
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>

#include "getopt.h"

#ifdef SIZEOF_LONG
# define BITS_PER_LONG (CHAR_BIT * SIZEOF_LONG)
# if BITS_PER_LONG > 32
#  define NEED_HANDLE_LARGE_LONG 1
# else
#  define NEED_HANDLE_LARGE_LONG 0
# endif
#else
# define BITS_PER_LONG (CHAR_BIT * sizeof(unsigned long))
# define NEED_HANDLE_LARGE_LONG 1
#endif


static void
usage(void)
{
  fprintf(stderr, "Usage: erathostenes [OPTIONS] [LIMIT]\n\n"
	  "Options:\n"
Niels Möller's avatar
Niels Möller committed
56
	  "      -?         Display this message.\n"
57 58 59
	  "      -b SIZE    Block size.\n"
	  "      -v         Verbose output.\n"
	  "      -s         No output.\n");
Niels Möller's avatar
Niels Möller committed
60 61 62
}

static unsigned
Niels Möller's avatar
Niels Möller committed
63
isqrt(unsigned long n)
Niels Möller's avatar
Niels Möller committed
64
{
Niels Möller's avatar
Niels Möller committed
65
  unsigned long x;
Niels Möller's avatar
Niels Möller committed
66 67

  /* FIXME: Better initialization. */
Niels Möller's avatar
Niels Möller committed
68
  if (n < ULONG_MAX)
Niels Möller's avatar
Niels Möller committed
69 70 71
    x = n;
  else
    /* Must avoid overflow in the first step. */
Niels Möller's avatar
Niels Möller committed
72
    x = n-1;
Niels Möller's avatar
Niels Möller committed
73 74 75

  for (;;)
    {
Niels Möller's avatar
Niels Möller committed
76
      unsigned long y = (x + n/x) / 2;
Niels Möller's avatar
Niels Möller committed
77 78 79 80 81 82 83 84 85
      if (y >= x)
	return x;

      x = y;
    }
}

/* Size is in bits */
static unsigned long *
Niels Möller's avatar
Niels Möller committed
86
vector_alloc(unsigned long size)
Niels Möller's avatar
Niels Möller committed
87
{
Niels Möller's avatar
Niels Möller committed
88
  unsigned long end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG;
Niels Möller's avatar
Niels Möller committed
89
  unsigned long *vector = malloc (end * sizeof(long));
Niels Möller's avatar
Niels Möller committed
90

Niels Möller's avatar
Niels Möller committed
91
  if (!vector)
92 93 94 95 96 97 98 99 100 101 102 103
    {
      fprintf(stderr, "Insufficient memory.\n");
      exit(EXIT_FAILURE);
    }
  return vector;
}

static void
vector_init(unsigned long *vector, unsigned long size)
{
  unsigned long end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG;
  unsigned long i;
Niels Möller's avatar
Niels Möller committed
104 105 106 107 108 109

  for (i = 0; i < end; i++)
    vector[i] = ~0;
}

static void
Niels Möller's avatar
Niels Möller committed
110 111
vector_clear_bits (unsigned long *vector, unsigned long step,
		   unsigned long start, unsigned long size)
Niels Möller's avatar
Niels Möller committed
112
{
Niels Möller's avatar
Niels Möller committed
113
  unsigned long bit;
Niels Möller's avatar
Niels Möller committed
114

115
  for (bit = start; bit < size; bit += step)
Niels Möller's avatar
Niels Möller committed
116
    {
Niels Möller's avatar
Niels Möller committed
117
      unsigned long i = bit / BITS_PER_LONG;
118
      unsigned long mask = 1L << (bit % BITS_PER_LONG);
Niels Möller's avatar
Niels Möller committed
119

120
      vector[i] &= ~mask;
Niels Möller's avatar
Niels Möller committed
121 122 123 124 125
    }
}

static unsigned
find_first_one (unsigned long x)
126 127
{  
  unsigned table[0x101] =
Niels Möller's avatar
Niels Möller committed
128
    {
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
     15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     14, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     12, 0, 0, 0, 0, 0, 0, 0,11, 0, 0, 0,10, 0, 9, 8,
      0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0,
      4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      7,
Niels Möller's avatar
Niels Möller committed
146 147 148 149
    };

  /* Isolate least significant bit */
  x &= -x;
Niels Möller's avatar
Niels Möller committed
150

Niels Möller's avatar
Niels Möller committed
151 152 153
  unsigned i = 0;
#if NEED_HANDLE_LARGE_LONG
#ifndef SIZEOF_LONG
154
  /* Can not be tested by the preprocessor. May generate warnings
Niels Möller's avatar
Niels Möller committed
155 156 157 158 159 160 161 162 163 164 165 166 167
     when long is 32 bits. */
  if (BITS_PER_LONG > 32)
#endif
    while (x >= 0x100000000L)
      {
	x >>= 32;
	i += 32;
      }
#endif /* NEED_HANDLE_LARGE_LONG */

  if (x >= 0x10000)
    {
      x >>= 16;
Niels Möller's avatar
Niels Möller committed
168
      i += 16;
Niels Möller's avatar
Niels Möller committed
169
    }
170
  return i + table[128 + (x & 0xff) - (x >> 8)];
Niels Möller's avatar
Niels Möller committed
171 172 173
}

/* Returns size if there's no more bits set */
Niels Möller's avatar
Niels Möller committed
174 175
static unsigned long
vector_find_next (const unsigned long *vector, unsigned long bit, unsigned long size)
Niels Möller's avatar
Niels Möller committed
176
{
Niels Möller's avatar
Niels Möller committed
177 178
  unsigned long end = (size + BITS_PER_LONG - 1) / BITS_PER_LONG;
  unsigned long i = bit / BITS_PER_LONG;
Niels Möller's avatar
Niels Möller committed
179 180 181 182 183 184 185 186 187 188 189 190 191 192
  unsigned long mask = 1L << (bit % BITS_PER_LONG);
  unsigned long word;

  if (i >= end)
    return size;

  for (word = vector[i] & ~(mask - 1); !word; word = vector[i])
    if (++i >= end)
      return size;

  /* Next bit is the least significant bit of word */
  return i * BITS_PER_LONG + find_first_one(word);
}

Niels Möller's avatar
Niels Möller committed
193 194 195 196
/* For benchmarking, define to do nothing (otherwise, most of the time
   will be spent converting the output to decimal). */
#define OUTPUT(n) printf("%lu\n", (n))

197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
static long
atosize(const char *s)
{
  char *end;
  long value = strtol(s, &end, 10);

  if (value <= 0)
    return 0;

  /* FIXME: Doesn't check for overflow. */
  switch(*end)
    {
    default:
      return 0;
    case '\0':
      break;
    case 'k': case 'K':
      value <<= 10;
      break;
    case 'M':
      value <<= 20;
      break;
    }
  return value;
}

Niels Möller's avatar
Niels Möller committed
223 224 225 226
int
main (int argc, char **argv)
{
  /* Generate all primes p <= limit */
Niels Möller's avatar
Niels Möller committed
227
  unsigned long limit;
228 229 230 231 232 233 234 235 236 237 238 239 240
  unsigned long root;

  unsigned long limit_nbits;

  /* Represents numbers up to sqrt(limit) */
  unsigned long sieve_nbits;
  unsigned long *sieve;
  /* Block for the rest of the sieving. Size should match the cache,
     the default value corresponds to 64 KB. */
  unsigned long block_nbits = 64L << 13;
  unsigned long block_start_bit;
  unsigned long *block;
  
Niels Möller's avatar
Niels Möller committed
241
  unsigned long bit;
242 243
  int silent = 0;
  int verbose = 0;
Niels Möller's avatar
Niels Möller committed
244
  int c;
Niels Möller's avatar
Niels Möller committed
245

246 247 248 249 250 251 252 253 254 255 256 257
  enum { OPT_HELP = 300 };
  static const struct option options[] =
    {
      /* Name, args, flag, val */
      { "help", no_argument, NULL, OPT_HELP },
      { "verbose", no_argument, NULL, 'v' },
      { "block-size", required_argument, NULL, 'b' },
      { "quiet", required_argument, NULL, 'q' },
      { NULL, 0, NULL, 0}
    };

  while ( (c = getopt_long(argc, argv, "svb:", options, NULL)) != -1)
Niels Möller's avatar
Niels Möller committed
258 259
    switch (c)
      {
260
      case OPT_HELP:
Niels Möller's avatar
Niels Möller committed
261
	usage();
262
	return EXIT_SUCCESS;
Niels Möller's avatar
Niels Möller committed
263
      case 'b':
264 265 266 267 268 269 270 271
	block_nbits = CHAR_BIT * atosize(optarg);
	if (!block_nbits)
	  {
	    usage();
	    return EXIT_FAILURE;
	  }
	break;

272
      case 'q':
273 274 275 276 277 278 279
	silent = 1;
	break;

      case 'v':
	verbose++;
	break;

280 281 282
      case '?':
	return EXIT_FAILURE;

Niels Möller's avatar
Niels Möller committed
283 284 285 286 287 288 289 290 291 292 293
      default:
	abort();
      }

  argc -= optind;
  argv += optind;

  if (argc == 0)
    limit = 1000;
  else if (argc == 1)
    {
294
      limit = atol(argv[0]);
Niels Möller's avatar
Niels Möller committed
295 296 297 298 299 300 301 302 303
      if (limit < 2)
	return EXIT_SUCCESS;
    }
  else
    {
      usage();
      return EXIT_FAILURE;
    }

304 305 306 307 308 309 310
  root = isqrt(limit);
  /* Round down to odd */
  root = (root - 1) | 1;
  /* Represents odd numbers from 3 up. */
  sieve_nbits = (root - 1) / 2;
  sieve = vector_alloc(sieve_nbits );
  vector_init(sieve, sieve_nbits);
311

312 313 314 315 316
  if (verbose)
    fprintf(stderr, "Initial sieve using %lu bits.\n", sieve_nbits);
      
  if (!silent)
    printf("2\n");
Niels Möller's avatar
Niels Möller committed
317 318 319 320

  if (limit == 2)
    return EXIT_SUCCESS;

321 322 323
  for (bit = 0;
       bit < sieve_nbits;
       bit = vector_find_next(sieve, bit + 1, sieve_nbits))
Niels Möller's avatar
Niels Möller committed
324
    {
Niels Möller's avatar
Niels Möller committed
325
      unsigned long n = 3 + 2 * bit;
Niels Möller's avatar
Niels Möller committed
326 327
      /* First bit to clear corresponds to n^2, which is bit

328 329 330 331 332 333
	 (n^2 - 3) / 2 = (n + 3) * bit + 3
      */      
      unsigned long n2_bit = (n+3)*bit + 3;

      if (!silent)
	printf("%lu\n", n);
Niels Möller's avatar
Niels Möller committed
334

335
      vector_clear_bits (sieve, n, n2_bit, sieve_nbits);
Niels Möller's avatar
Niels Möller committed
336 337
    }

338
  limit_nbits = (limit - 1) / 2;
Niels Möller's avatar
Niels Möller committed
339

340 341
  if (sieve_nbits + block_nbits > limit_nbits)
    block_nbits = limit_nbits - sieve_nbits;
Niels Möller's avatar
Niels Möller committed
342

343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
  if (verbose)
    {
      double storage = block_nbits / 8.0;
      unsigned shift = 0;
      const char prefix[] = " KMG";

      while (storage > 1024 && shift < 3)
	{
	  storage /= 1024;
	  shift++;
	}
      fprintf(stderr, "Blockwise sieving using blocks of %lu bits (%.3g %cByte)\n",
	      block_nbits, storage, prefix[shift]);
    }

  block = vector_alloc(block_nbits);

  for (block_start_bit = bit; block_start_bit < limit_nbits; block_start_bit += block_nbits)
361
    {
Niels Möller's avatar
Niels Möller committed
362
      unsigned long block_start;
363 364 365
      
      if (block_start_bit + block_nbits > limit_nbits)
	block_nbits = limit_nbits - block_start_bit;
366

367
      vector_init(block, block_nbits);
368

369
      block_start = 3 + 2*block_start_bit;
370

371 372 373 374 375 376 377
      if (verbose > 1)
	fprintf(stderr, "Next block, n = %lu\n", block_start);

      /* Sieve */
      for (bit = 0; bit < sieve_nbits;
	   bit = vector_find_next(sieve, bit + 1, sieve_nbits))
	{	  
Niels Möller's avatar
Niels Möller committed
378
	  unsigned long n = 3 + 2 * bit;
379
	  unsigned long sieve_start_bit = (n + 3) * bit + 3;
380

381
	  if (sieve_start_bit < block_start_bit)
382
	    {
383 384 385 386
	      unsigned long k = (block_start + n - 1) / (2*n);
	      sieve_start_bit = n * k + bit;

	      assert(sieve_start_bit < block_start_bit + n);
387
	    }
388
	  assert(sieve_start_bit >= block_start_bit);
Niels Möller's avatar
Niels Möller committed
389

390 391 392 393 394 395 396 397 398
	  vector_clear_bits(block, n, sieve_start_bit - block_start_bit, block_nbits);
	}
      for (bit = vector_find_next(block, 0, block_nbits);
	   bit < block_nbits;
	   bit = vector_find_next(block, bit + 1, block_nbits))
	{
	  unsigned long n = block_start + 2 * bit;
	  if (!silent)
	    printf("%lu\n", n);
399
	}
Niels Möller's avatar
Niels Möller committed
400 401 402
    }
  return EXIT_SUCCESS;
}