eratosthenes.c 9.26 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
  static const unsigned char 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
  unsigned i = 0;

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

Niels Möller's avatar
Niels Möller committed
153 154
#if NEED_HANDLE_LARGE_LONG
#ifndef SIZEOF_LONG
155
  /* Can not be tested by the preprocessor. May generate warnings
Niels Möller's avatar
Niels Möller committed
156 157 158 159 160 161 162 163 164 165 166 167 168
     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
169
      i += 16;
Niels Möller's avatar
Niels Möller committed
170
    }
171
  return i + table[128 + (x & 0xff) - (x >> 8)];
Niels Möller's avatar
Niels Möller committed
172 173 174
}

/* Returns size if there's no more bits set */
Niels Möller's avatar
Niels Möller committed
175 176
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
177
{
Niels Möller's avatar
Niels Möller committed
178 179
  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
180 181 182 183 184 185 186 187 188 189 190 191 192 193
  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
194 195 196 197
/* 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))

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 223
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
224 225 226 227
int
main (int argc, char **argv)
{
  /* Generate all primes p <= limit */
Niels Möller's avatar
Niels Möller committed
228
  unsigned long limit;
229 230 231 232 233 234 235 236 237 238 239 240 241
  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
242
  unsigned long bit;
243 244
  int silent = 0;
  int verbose = 0;
Niels Möller's avatar
Niels Möller committed
245
  int c;
Niels Möller's avatar
Niels Möller committed
246

247 248 249 250 251 252 253 254 255 256 257 258
  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
259 260
    switch (c)
      {
261
      case OPT_HELP:
Niels Möller's avatar
Niels Möller committed
262
	usage();
263
	return EXIT_SUCCESS;
Niels Möller's avatar
Niels Möller committed
264
      case 'b':
265 266 267 268 269 270 271 272
	block_nbits = CHAR_BIT * atosize(optarg);
	if (!block_nbits)
	  {
	    usage();
	    return EXIT_FAILURE;
	  }
	break;

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

      case 'v':
	verbose++;
	break;

281 282 283
      case '?':
	return EXIT_FAILURE;

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

  argc -= optind;
  argv += optind;

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

305 306 307 308 309 310 311
  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);
312

313 314 315 316 317
  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
318 319 320 321

  if (limit == 2)
    return EXIT_SUCCESS;

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

329 330 331 332 333 334
	 (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
335

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

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

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

344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361
  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)
362
    {
Niels Möller's avatar
Niels Möller committed
363
      unsigned long block_start;
364 365 366
      
      if (block_start_bit + block_nbits > limit_nbits)
	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 373 374 375 376 377 378
      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
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 385 386 387
	      unsigned long k = (block_start + n - 1) / (2*n);
	      sieve_start_bit = n * k + bit;

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

391 392 393 394 395 396 397 398 399
	  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);
400
	}
Niels Möller's avatar
Niels Möller committed
401 402 403
    }
  return EXIT_SUCCESS;
}