eratosthenes.c 6.44 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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
  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;
      i =+ 16;
    }
  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
179
180
181
182
183
  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);
}

int
main (int argc, char **argv)
{
  unsigned long *vector;
  /* Generate all primes p <= limit */
Niels Möller's avatar
Niels Möller committed
184
185
186
187
188
189
  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
190

Niels Möller's avatar
Niels Möller committed
191
  int c;
Niels Möller's avatar
Niels Möller committed
192

193
  block_size = 0;
Niels Möller's avatar
Niels Möller committed
194

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

235
236
237
  if (!block_size || block_size > size)
    block_size = size;

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

Niels Möller's avatar
Niels Möller committed
245
  printf("2\n");
Niels Möller's avatar
Niels Möller committed
246
247
248
249
250
251

  bit = 0;

  if (limit == 2)
    return EXIT_SUCCESS;

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

Niels Möller's avatar
Niels Möller committed
254
255
  while (bit < sieve_limit)
    {
Niels Möller's avatar
Niels Möller committed
256
      unsigned long n = 3 + 2 * bit;
Niels Möller's avatar
Niels Möller committed
257

Niels Möller's avatar
Niels Möller committed
258
      printf("%lu\n", n);
Niels Möller's avatar
Niels Möller committed
259
260
261
262
263

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

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

266
      bit = vector_find_next (vector, bit + 1, block_size);
Niels Möller's avatar
Niels Möller committed
267
268
269
    }

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

Niels Möller's avatar
Niels Möller committed
273
    printf("%lu\n", 3 + 2 * bit);
Niels Möller's avatar
Niels Möller committed
274

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

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

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

300
301
302
303
304
	  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
305
	printf("%lu\n", 3 + 2 * bit);
Niels Möller's avatar
Niels Möller committed
306
307
308
309
    }

  return EXIT_SUCCESS;
}