"""Build the segmented sieve used for validation (C source embedded below).

Usage:
    python3 sieve_build.py            # writes sieve.c next to this file and compiles ./sieve with gcc
    ./sieve 10000000000 401 1e9 1e10 > cp_1e9_1e10.txt

Output columns at each of the 401 log-spaced checkpoints x in [1e9, 1e10]:
    x, theta(x), sum_{p<=x} log(1-1/p), sum_{p<=x} log(1-1/p^2), pi(x)
(accumulators in binary128). Requires gcc and libquadmath.
"""
import os, subprocess
SIEVE_C = r'''// Segmented sieve to LIMIT; at checkpoints x_j (sorted reals) output
// theta(x), M(x)=sum_{p<=x} log(1-1/p), Q2(x)=sum_{p<=x} log(1-1/p^2), pi(x)
// accumulators in __float128; per-term values in long double.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stdint.h>
#include <quadmath.h>
int main(int argc,char**argv){
  uint64_t LIMIT=strtoull(argv[1],0,10);
  int NC=atoi(argv[2]); double lo=atof(argv[3]), hi=atof(argv[4]);
  uint64_t *cp=malloc(sizeof(uint64_t)*NC);
  for(int j=0;j<NC;j++){ double x=exp(log(lo)+(log(hi)-log(lo))*j/(NC-1)); cp[j]=(uint64_t)floor(x);}
  uint64_t r=(uint64_t)sqrtl((long double)LIMIT)+2;
  char*small=calloc(r+1,1); uint64_t*sp=malloc(sizeof(uint64_t)*r); int ns=0;
  for(uint64_t i=2;i<=r;i++) if(!small[i]){ sp[ns++]=i; for(uint64_t j=i*i;j<=r;j+=i) small[j]=1; }
  const uint64_t SEG=1<<22; char*seg=malloc(SEG);
  __float128 th=0,M=0,Q2=0; uint64_t pi=0; int j=0;
  for(uint64_t L=2; L<=LIMIT; L+=SEG){
    uint64_t H=L+SEG-1; if(H>LIMIT) H=LIMIT;
    memset(seg,0,H-L+1);
    for(int k=0;k<ns;k++){ uint64_t p=sp[k]; if(p*p>H) break;
      uint64_t s=(L+p-1)/p*p; if(s<p*p) s=p*p;
      for(uint64_t m=s;m<=H;m+=p) seg[m-L]=1; }
    for(uint64_t n=L;n<=H;n++){
      while(j<NC && cp[j]<n){ {char a[64],b[64],c[64]; quadmath_snprintf(a,64,"%.32Qe",th); quadmath_snprintf(b,64,"%.32Qe",M); quadmath_snprintf(c,64,"%.32Qe",Q2); printf("%llu %s %s %s %llu\n",(unsigned long long)cp[j],a,b,c,(unsigned long long)pi);} j++; }
      if(!seg[n-L]){ long double p=(long double)n; th+=logl(p); M+=log1pl(-1.0L/p); Q2+=log1pl(-1.0L/(p*p)); pi++; }
    }
  }
  while(j<NC && cp[j]<=LIMIT){ {char a[64],b[64],c[64]; quadmath_snprintf(a,64,"%.32Qe",th); quadmath_snprintf(b,64,"%.32Qe",M); quadmath_snprintf(c,64,"%.32Qe",Q2); printf("%llu %s %s %s %llu\n",(unsigned long long)cp[j],a,b,c,(unsigned long long)pi);} j++; }
  return 0;
}
'''
if __name__ == "__main__":
    here = os.path.dirname(os.path.abspath(__file__))
    c = os.path.join(here, "sieve.c")
    with open(c, "w") as f:
        f.write(SIEVE_C)
    subprocess.run(["gcc", "-O2", "-o", os.path.join(here, "sieve"), c, "-lquadmath", "-lm"], check=True)
    print("built", os.path.join(here, "sieve"))