-
Notifications
You must be signed in to change notification settings - Fork 0
/
p10.c
155 lines (122 loc) · 4 KB
/
p10.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
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
/* Project Euler
*
* http://projecteuler.net/index.php?section=problems&id=10
*
* Problem 10
* 08 February 2002
*
* The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.
*
* Find the sum of all the primes below two million.
* Answer:
* 37550402023 (was: one million)
* 142913828922
*
* FIXME: worth to try: http://en.wikipedia.org/wiki/Sieve_of_Atkin
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define BYTES_PER_INT 4
#define BITS_PER_INT 32
#define ALL_ON 0xffffffff
static const unsigned int Nth_bit[] = {
0x01 << 0,
0x01 << 1,
0x01 << 2,
0x01 << 3,
0x01 << 4,
0x01 << 5,
0x01 << 6,
0x01 << 7,
0x01 << 8,
0x01 << 9,
0x01 << 10,
0x01 << 11,
0x01 << 12,
0x01 << 13,
0x01 << 14,
0x01 << 15,
0x01 << 16,
0x01 << 17,
0x01 << 18,
0x01 << 19,
0x01 << 20,
0x01 << 21,
0x01 << 22,
0x01 << 23,
0x01 << 24,
0x01 << 25,
0x01 << 26,
0x01 << 27,
0x01 << 28,
0x01 << 29,
0x01 << 30,
0x01 << 31
};
#define NTH_BIT(bit) Nth_bit[bit]
/* return an integer containing only Nth bit turend on
* (where N is in 0..31 range)
*/
#define BIT_TO_INDEX(bit) ( (bit) / BITS_PER_INT )
/* return index of array member containing given bit */
#define BIT_TO_BIT(bit) ( (bit) % BITS_PER_INT )
/* return number of the bit (in 0..31 range) inside of
* array member corresponding given bit
*/
#define SET_BIT(array, bit) array[ BIT_TO_INDEX(bit) ] |= NTH_BIT( BIT_TO_BIT(bit) )
/* turn on given bit in bitarray */
#define UNSET_BIT(array, bit) array[ BIT_TO_INDEX(bit) ] &= ~NTH_BIT( BIT_TO_BIT(bit) )
/* turn off given bit in bitarray */
#define IS_BIT_SET(array, bit) ( array[ BIT_TO_INDEX(bit) ] & NTH_BIT( BIT_TO_BIT(bit) ) )
/* return TRUE if given bit is turned on in bitarray */
#define PRIMES_UPPER_BOUND 2000000
/* max value to test for the primes */
int main(void)
{
unsigned int *bitarray = NULL;
unsigned int bitarray_size;
unsigned int i, j, n;
double s;
if (sizeof(unsigned int) != BYTES_PER_INT) {
fprintf(stderr, "size of 'unsigned int' is expected to be %d instead of %d\n",
BYTES_PER_INT,
sizeof(unsigned int));
return 1;
}
/* this is how many 'unsigned int's needed to store PRIMES_UPPER_BOUND bits */
bitarray_size = PRIMES_UPPER_BOUND / BITS_PER_INT
+ ((PRIMES_UPPER_BOUND % BITS_PER_INT) == 0 ? 0 : 1);
/* printf("PRIMES_UPPER_BOUND=%d\n", PRIMES_UPPER_BOUND); */
/* printf("BITS_PER_INT=%d\n", BITS_PER_INT); */
/* printf("bitarray_size=%d\n", bitarray_size); */
/* allocate bitarray and turn all the bits on */
bitarray = (unsigned int*)malloc(sizeof(unsigned int)*bitarray_size);
for (i = 0; i < bitarray_size; i++)
bitarray[i] = ALL_ON;
/* determine prime numbers by Sieve of Eratosthenes */
for (i = 2; i <= PRIMES_UPPER_BOUND; i++) {
if (! IS_BIT_SET(bitarray, i))
/* i is not a prime */
continue;
/* i is prime: remove from the list all the numbers
* divisible by i except itself
*/
/* printf("prime %6d\n", i); */
for (n = 2; ; n++) {
j = i*n;
if (j >= PRIMES_UPPER_BOUND)
break;
UNSET_BIT(bitarray, j);
}
}
/* compute sum of all the primes */
s = 0;
for (i = 2; i < PRIMES_UPPER_BOUND; i++)
if (IS_BIT_SET(bitarray, i))
s += i;
printf("%.0f\n", s);
free(bitarray);
return 0;
}
/* end of file */