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
|
#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
#include <string.h>
#include <errno.h>
#include <stdlib.h>
#include <math.h>
#ifdef __GNUC__
#define GCC_ALWAYS_INLINE inline __attribute__((always_inline))
#else
#define GCC_ALWAYS_INLINE inline
#endif
/* finds the popcount of a number without using a fancy intrinsic */
GCC_ALWAYS_INLINE unsigned fake_pc(uint64_t pc) {
pc = (pc & UINT64_C(0x5555555555555555)) + ((pc & UINT64_C(0xAAAAAAAAAAAAAAAA)) >> 1);
pc = (pc & UINT64_C(0x3333333333333333)) + ((pc & UINT64_C(0xCCCCCCCCCCCCCCCC)) >> 2);
pc = (pc & UINT64_C(0x0F0F0F0F0F0F0F0F)) + ((pc & UINT64_C(0xF0F0F0F0F0F0F0F0)) >> 4);
pc = (pc & UINT64_C(0x00FF00FF00FF00FF)) + ((pc & UINT64_C(0xFF00FF00FF00FF00)) >> 8);
pc = (pc & UINT64_C(0x0000FFFF0000FFFF)) + ((pc & UINT64_C(0xFFFF0000FFFF0000)) >> 16);
pc = (pc & UINT64_C(0x00000000FFFFFFFF)) + ((pc & UINT64_C(0xFFFFFFFF00000000)) >> 32);
return (unsigned)pc; /* it will always be small enough to fit into an unsigned int (it won't be more than 63) */
}
GCC_ALWAYS_INLINE unsigned real_pc(uint64_t pc) {
return (unsigned)(__builtin_popcount((unsigned)(pc & 0xFFFFFFFF))
+ __builtin_popcount((unsigned)(pc >> 32)));
}
#ifndef POPCOUNT
#define POPCOUNT fake_pc
#endif
#define PRIME_LONG(v, _st) ((_st)[(v) >> 6])
#define PRIME_SHIFT(v) ((v) & 63)
#define PRIME_NUM(v, _st) ((PRIME_LONG(v, _st) >> PRIME_SHIFT(v)) & 1)
unsigned long calc_pi(unsigned long in) {
uint64_t *storage;
unsigned long anum = in-1;
unsigned long nlongs = (anum + 63) >> 6;
unsigned long sum = 0;
storage = calloc(nlongs, sizeof(uint64_t));
if (!storage) {
fprintf(stderr, "Could not allocate space for %zu longs! (about %zu bytes): %s\n",
nlongs, nlongs * sizeof(uint64_t), strerror(errno));
abort();
}
/* fill storage with 1 bits */
for (size_t i = 0; i < nlongs - 1; ++i) {
storage[i] = ~UINT64_C(0);
}
storage[nlongs-1] = (UINT64_C(1) << (anum & 63)) - 1;
/* the sieve */
for (uint64_t n = 2; n * n <= in; ++n) {
if (!PRIME_NUM(n - 2, storage)) continue;
for (uint64_t m = n * 2; m <= in; m += n) {
PRIME_LONG(m - 2, storage) &= ~(UINT64_C(1) << PRIME_SHIFT(m - 2));
}
}
/* count up the primes */
for (size_t n = 0; n < nlongs; ++n) {
sum += POPCOUNT(storage[n]);
}
free(storage);
return sum;
}
unsigned long get_num(const char *str) {
char *end;
unsigned long ret = strtoul(str, &end, 10);
if (!*str || *end) {
fputs("Invalid number entered :(\n", stderr);
exit(EXIT_FAILURE);
}
if (ret < 2) {
fputs("Number must not be less than 2\n", stderr);
exit(EXIT_FAILURE);
}
return ret;
}
int main(int argc, char **argv) {
if (argc < 2) {
fprintf(stderr, "Usage: %s <num>\nWill print out the number of primes less than or equal to the supplied number.\n", argv[0] ? argv[0] : "(null)");
return EXIT_FAILURE;
}
unsigned long num = get_num(argv[1]);
unsigned long sum = calc_pi(num);
double v = (sum * log(num)) / num;
printf("%lu primes\n", sum);
printf("%lf num\n", v);
return EXIT_SUCCESS;
}
|
