-
Notifications
You must be signed in to change notification settings - Fork 0
/
sieve-of-atkins.c
100 lines (86 loc) · 2.67 KB
/
sieve-of-atkins.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
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
void sieve_of_atkin(mpz_t limit) {
mpz_t x2, y2, n, limit_sqrt, remainder, temp, k;
mpz_init(x2);
mpz_init(y2);
mpz_init(n);
mpz_init(limit_sqrt);
mpz_init(remainder);
mpz_init(temp);
mpz_init(k);
int limit_int = mpz_get_ui(limit);
char* sieve = calloc(limit_int + 1, sizeof(char));
mpz_sqrt(limit_sqrt, limit);
for (mpz_set_ui(n, 1); mpz_cmp(n, limit_sqrt) <= 0; mpz_add_ui(n, n, 1)) {
mpz_mul(x2, n, n);
for (mpz_set_ui(y2, 1); mpz_cmp(y2, limit_sqrt) <= 0; mpz_add_ui(y2, y2, 1)) {
mpz_mul(y2, y2, y2);
mpz_set(temp, x2);
mpz_add(temp, temp, y2);
if (mpz_cmp(temp, limit) <= 0) {
mpz_mod_ui(remainder, temp, 12);
int r = mpz_get_ui(remainder);
if (r == 1 || r == 5) {
sieve[mpz_get_ui(temp)] ^= 1;
}
}
mpz_set(temp, x2);
mpz_mul_ui(temp, temp, 3);
mpz_sub(temp, temp, y2);
if (mpz_cmp_ui(temp, 0) > 0 && mpz_cmp(temp, limit) <= 0) {
mpz_mod_ui(remainder, temp, 12);
int r = mpz_get_ui(remainder);
if (r == 7) {
sieve[mpz_get_ui(temp)] ^= 1;
}
}
mpz_set(temp, x2);
mpz_mul_ui(temp, temp, 3);
mpz_add(temp, temp, y2);
if (mpz_cmp(temp, limit) <= 0) {
mpz_mod_ui(remainder, temp, 12);
int r = mpz_get_ui(remainder);
if (r == 11) {
sieve[mpz_get_ui(temp)] ^= 1;
}
}
}
}
for (mpz_set_ui(n, 5); mpz_cmp(n, limit_sqrt) <= 0; mpz_add_ui(n, n, 1)) {
if (sieve[mpz_get_ui(n)]) {
mpz_mul(temp, n, n);
for (mpz_set(k, temp); mpz_cmp(k, limit) <= 0; mpz_add(k, k, temp)) {
sieve[mpz_get_ui(k)] = 0;
}
}
}
if (limit_int >= 2) printf("2 ");
if (limit_int >= 3) printf("3 ");
for (mpz_set_ui(n, 5); mpz_cmp(n, limit) <= 0; mpz_add_ui(n, n, 1)) {
if (sieve[mpz_get_ui(n)]) {
gmp_printf("%Zd ", n);
}
}
printf("\n");
mpz_clear(x2);
mpz_clear(y2);
mpz_clear(n);
mpz_clear(limit_sqrt);
mpz_clear(remainder);
mpz_clear(temp);
mpz_clear(k);
free(sieve);
}
int main(int argc, char* argv[]) {
if (argc != 2) {
fprintf(stderr, "Usage: %s <limit>\n", argv[0]);
return 1;
}
mpz_t limit;
mpz_init_set_str(limit, argv[1], 10);
sieve_of_atkin(limit);
mpz_clear(limit);
return 0;
}