eratosthenes.c 7.38 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
    }
}

static unsigned
find_first_one (unsigned long x)
115 116
{  
  unsigned table[0x101] =
Niels Möller's avatar
Niels Möller committed
117
    {
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
     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
135 136 137 138
    };

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

Niels Möller's avatar
Niels Möller committed
140 141 142
  unsigned i = 0;
#if NEED_HANDLE_LARGE_LONG
#ifndef SIZEOF_LONG
143
  /* Can not be tested by the preprocessor. May generate warnings
Niels Möller's avatar
Niels Möller committed
144 145 146 147 148 149 150 151 152 153 154 155 156
     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
157
      i += 16;
Niels Möller's avatar
Niels Möller committed
158
    }
159
  return i + table[128 + (x & 0xff) - (x >> 8)];
Niels Möller's avatar
Niels Möller committed
160 161 162
}

/* Returns size if there's no more bits set */
Niels Möller's avatar
Niels Möller committed
163 164
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
165
{
Niels Möller's avatar
Niels Möller committed
166 167
  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
168 169 170 171 172 173 174 175 176 177 178 179 180 181
  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
182 183 184 185
/* 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
186 187 188 189 190
int
main (int argc, char **argv)
{
  unsigned long *vector;
  /* Generate all primes p <= limit */
Niels Möller's avatar
Niels Möller committed
191 192 193 194 195 196
  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
197

Niels Möller's avatar
Niels Möller committed
198
  int c;
Niels Möller's avatar
Niels Möller committed
199

200
  block_size = 0;
Niels Möller's avatar
Niels Möller committed
201

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

  argc -= optind;
  argv += optind;

  if (argc == 0)
    limit = 1000;
  else if (argc == 1)
    {
230
      limit = atol(argv[0]);
Niels Möller's avatar
Niels Möller committed
231 232 233 234 235 236 237 238 239 240 241
      if (limit < 2)
	return EXIT_SUCCESS;
    }
  else
    {
      usage();
      return EXIT_FAILURE;
    }

  size = (limit - 1) / 2;

242 243 244
  if (!block_size || block_size > size)
    block_size = size;

245 246 247 248
  /* FIXME: Allocates a bit for all odd numbers up to limit. Could
     save a lot of memory by allocating space only for numbers up to
     the square root, plus one additional block. */
  vector = vector_alloc(size);
Niels Möller's avatar
Niels Möller committed
249 250 251 252 253 254
  if (!vector)
    {
      fprintf(stderr, "Insufficient memory.\n");
      return EXIT_FAILURE;
    }

Niels Möller's avatar
Niels Möller committed
255
  OUTPUT(2UL);
Niels Möller's avatar
Niels Möller committed
256 257 258 259 260 261

  bit = 0;

  if (limit == 2)
    return EXIT_SUCCESS;

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

Niels Möller's avatar
Niels Möller committed
264 265
  while (bit < sieve_limit)
    {
Niels Möller's avatar
Niels Möller committed
266
      unsigned long n = 3 + 2 * bit;
Niels Möller's avatar
Niels Möller committed
267

Niels Möller's avatar
Niels Möller committed
268
      OUTPUT(n);
Niels Möller's avatar
Niels Möller committed
269 270 271 272 273

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

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

276
      bit = vector_find_next (vector, bit + 1, block_size);
Niels Möller's avatar
Niels Möller committed
277 278 279
    }

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

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

285 286
  for (block = block_size; block < size; block += block_size)
    {
Niels Möller's avatar
Niels Möller committed
287 288
      unsigned long block_start;
      unsigned long block_end;
289 290 291 292 293 294 295 296 297 298 299

      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
300
	  unsigned long n = 3 + 2 * bit;
301

Niels Möller's avatar
Niels Möller committed
302
	  unsigned long start = n*bit + 3*(bit + 1);
303 304
	  if (start < block)
	    {
Niels Möller's avatar
Niels Möller committed
305
	      unsigned long k = (block + 1) / n;
306 307 308
	      start = bit + k*n;
	    }
	  vector_clear_bits (vector, n, start, block + block_size);
Niels Möller's avatar
Niels Möller committed
309

310 311 312 313 314
	  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
315
	OUTPUT(3 + 2 * bit);
Niels Möller's avatar
Niels Möller committed
316 317 318 319
    }

  return EXIT_SUCCESS;
}