eratosthenes.c 6.58 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 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
 * 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

#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
55 56
	  "      -?         Display this message.\n"
	  "      -b  SIZE   Block size.\n");
Niels Möller's avatar
Niels Möller committed
57 58 59
}

static unsigned
Niels Möller's avatar
Niels Möller committed
60
isqrt(unsigned long n)
Niels Möller's avatar
Niels Möller committed
61
{
Niels Möller's avatar
Niels Möller committed
62
  unsigned long x;
Niels Möller's avatar
Niels Möller committed
63 64

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

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

      x = y;
    }
}

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

Niels Möller's avatar
Niels Möller committed
89 90 91 92 93 94 95 96 97 98
  if (!vector)
    return NULL;

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

  return vector;
}

static void
Niels Möller's avatar
Niels Möller committed
99 100
vector_clear_bits (unsigned long *vector, unsigned long step,
		   unsigned long start, unsigned long size)
Niels Möller's avatar
Niels Möller committed
101
{
Niels Möller's avatar
Niels Möller committed
102
  unsigned long bit;
Niels Möller's avatar
Niels Möller committed
103

104
  for (bit = start; bit < size; bit += step)
Niels Möller's avatar
Niels Möller committed
105
    {
Niels Möller's avatar
Niels Möller committed
106
      unsigned long i = bit / BITS_PER_LONG;
107
      unsigned long mask = 1L << (bit % BITS_PER_LONG);
Niels Möller's avatar
Niels Möller committed
108

109
      vector[i] &= ~mask;
Niels Möller's avatar
Niels Möller committed
110 111 112 113 114 115 116 117 118
    }
}

static unsigned
find_first_one (unsigned long x)
{
  unsigned table[0x10] =
    {
      /* 0, 1,  2,  3,  4,  5,  6,  7 */
Niels Möller's avatar
Niels Möller committed
119
	-1, 0,  1,  0,  2,  0,  1 , 0,
Niels Möller's avatar
Niels Möller committed
120 121 122 123 124 125
      /* 8, 9, 10, 11, 12, 13, 14, 15 */
	 3, 0,  1,  0,  2,  0,  1,  0
    };

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

Niels Möller's avatar
Niels Möller committed
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
  unsigned i = 0;
#if NEED_HANDLE_LARGE_LONG
#ifndef SIZEOF_LONG
  /* Can't not be tested by the preprocessor. May generate warnings
     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
144
      i += 16;
Niels Möller's avatar
Niels Möller committed
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
    }
  if (x >= 0x100)
    {
      x >>= 8;
      i += 8;
    }
  if (x >= 0x10)
    {
      x >>= 4;
      i += 4;
    }
  return i + table[x & 0xf];
}

/* Returns size if there's no more bits set */
Niels Möller's avatar
Niels Möller committed
160 161
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
162
{
Niels Möller's avatar
Niels Möller committed
163 164
  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
165 166 167 168 169 170 171 172 173 174 175 176 177 178
  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
179 180 181 182
/* 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))

Niels Möller's avatar
Niels Möller committed
183 184 185 186 187
int
main (int argc, char **argv)
{
  unsigned long *vector;
  /* Generate all primes p <= limit */
Niels Möller's avatar
Niels Möller committed
188 189 190 191 192 193
  unsigned long limit;
  unsigned long size;
  unsigned long bit;
  unsigned long sieve_limit;
  unsigned long block_size;
  unsigned long block;
Niels Möller's avatar
Niels Möller committed
194

Niels Möller's avatar
Niels Möller committed
195
  int c;
Niels Möller's avatar
Niels Möller committed
196

197
  block_size = 0;
Niels Möller's avatar
Niels Möller committed
198

Niels Möller's avatar
Niels Möller committed
199
  while ( (c = getopt(argc, argv, "?b:")) != -1)
Niels Möller's avatar
Niels Möller committed
200 201 202 203 204
    switch (c)
      {
      case '?':
	usage();
	return EXIT_FAILURE;
Niels Möller's avatar
Niels Möller committed
205
      case 'b':
206
	{
Niels Möller's avatar
Niels Möller committed
207
	  long arg = atoi(optarg);
208 209 210 211 212 213 214 215
	  if (arg <= 10)
	    {
	      usage();
	      return EXIT_FAILURE;
	    }
	  block_size = (arg - 3) / 2;
	  break;
	}
Niels Möller's avatar
Niels Möller committed
216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238
      default:
	abort();
      }

  argc -= optind;
  argv += optind;

  if (argc == 0)
    limit = 1000;
  else if (argc == 1)
    {
      limit = atoi(argv[0]);
      if (limit < 2)
	return EXIT_SUCCESS;
    }
  else
    {
      usage();
      return EXIT_FAILURE;
    }

  size = (limit - 1) / 2;

239 240 241
  if (!block_size || block_size > size)
    block_size = size;

Niels Möller's avatar
Niels Möller committed
242 243 244 245 246 247 248
  vector = vector_alloc (size);
  if (!vector)
    {
      fprintf(stderr, "Insufficient memory.\n");
      return EXIT_FAILURE;
    }

Niels Möller's avatar
Niels Möller committed
249
  OUTPUT(2UL);
Niels Möller's avatar
Niels Möller committed
250 251 252 253 254 255

  bit = 0;

  if (limit == 2)
    return EXIT_SUCCESS;

256
  sieve_limit = (isqrt(2*block_size + 1) - 1) / 2;
Niels Möller's avatar
Niels Möller committed
257

Niels Möller's avatar
Niels Möller committed
258 259
  while (bit < sieve_limit)
    {
Niels Möller's avatar
Niels Möller committed
260
      unsigned long n = 3 + 2 * bit;
Niels Möller's avatar
Niels Möller committed
261

Niels Möller's avatar
Niels Möller committed
262
      OUTPUT(n);
Niels Möller's avatar
Niels Möller committed
263 264 265 266 267

      /* First bit to clear corresponds to n^2, which is bit

	 (n^2 - 3) / 2 = n * bit + 3 (bit + 1)
      */
268
      vector_clear_bits (vector, n, n*bit + 3*(bit + 1), block_size);
Niels Möller's avatar
Niels Möller committed
269

270
      bit = vector_find_next (vector, bit + 1, block_size);
Niels Möller's avatar
Niels Möller committed
271 272 273
    }

  /* No more marking, just output the remaining primes. */
Niels Möller's avatar
Niels Möller committed
274
  for (; bit < block_size ;
Niels Möller's avatar
Niels Möller committed
275
       bit = vector_find_next (vector, bit + 1, block_size))
Niels Möller's avatar
Niels Möller committed
276

Niels Möller's avatar
Niels Möller committed
277
    OUTPUT(3 + 2 * bit);
Niels Möller's avatar
Niels Möller committed
278

279 280
  for (block = block_size; block < size; block += block_size)
    {
Niels Möller's avatar
Niels Möller committed
281 282
      unsigned long block_start;
      unsigned long block_end;
283 284 285 286 287 288 289 290 291 292 293

      if (block + block_size > size)
	/* For the final block */
	block_size = size - block;

      block_start = 2*block + 3;
      block_end = 2*(block + block_size) + 3;

      sieve_limit = (isqrt(block_end) - 1) / 2;
      for (bit = 0; bit < sieve_limit ;)
	{
Niels Möller's avatar
Niels Möller committed
294
	  unsigned long n = 3 + 2 * bit;
295

Niels Möller's avatar
Niels Möller committed
296
	  unsigned long start = n*bit + 3*(bit + 1);
297 298
	  if (start < block)
	    {
Niels Möller's avatar
Niels Möller committed
299
	      unsigned long k = (block + 1) / n;
300 301 302
	      start = bit + k*n;
	    }
	  vector_clear_bits (vector, n, start, block + block_size);
Niels Möller's avatar
Niels Möller committed
303

304 305 306 307 308
	  bit = vector_find_next (vector, bit + 1, block + block_size);
	}
      for (bit = vector_find_next (vector, block, block + block_size);
	   bit < block + block_size;
	   bit = vector_find_next (vector, bit + 1, block + block_size))
Niels Möller's avatar
Niels Möller committed
309
	OUTPUT(3 + 2 * bit);
Niels Möller's avatar
Niels Möller committed
310 311 312 313
    }

  return EXIT_SUCCESS;
}