fn xor(row: &mut [u64], pivot: &[u64])Examples found in repository?
crates/competitive/src/math/bit_matrix.rs (line 198)
173 fn eliminate_impl(&mut self, cols: usize, full: bool, require_full_rank: bool) -> Vec<usize> {
174 let n = self.shape.0;
175 let mut pivots = Vec::with_capacity(n.min(cols));
176 if n < 32 {
177 let mut c = 0;
178 while c < cols {
179 let r = pivots.len();
180 if r == n {
181 break;
182 }
183 let Some(p) = (r..n).find(|&i| self[i].get(c)) else {
184 if require_full_rank {
185 return pivots;
186 }
187 c = self.next_column(r, c + 1, cols);
188 continue;
189 };
190 self.data.swap(r, p);
191 let (upper, lower) = self.data.split_at_mut(r);
192 let (pivot, lower) = lower.split_first_mut().unwrap();
193 for row in lower
194 .iter_mut()
195 .chain(upper.iter_mut().take(if full { r } else { 0 }))
196 {
197 if row.get(c) {
198 xor(&mut row.words_mut()[c / 64..], &pivot.words()[c / 64..]);
199 }
200 }
201 pivots.push(c);
202 c += 1;
203 }
204 return pivots;
205 }
206 // Eliminating a shared leading coefficient preserves the two-coefficient bound.
207 if self
208 .data
209 .iter()
210 .all(|row| row.iter_ones().take_while(|&c| c < cols).take(3).count() <= 2)
211 {
212 return self.eliminate_sparse(cols, full, require_full_rank);
213 }
214 let block: usize = if n < 512 {
215 4
216 } else if n < 1536 {
217 8
218 } else {
219 32
220 };
221 let mut table = vec![BitSet::new(self.shape.1); (1 << block.min(8)) * block.div_ceil(8)];
222 let mut reduced = vec![0; n];
223 let mut start = 0;
224 while start < cols {
225 let first = pivots.len();
226 let word = start / 64;
227 reduced.fill(first);
228 let end = cols.min(start + block);
229 let mut c = start;
230 while c < end {
231 let r = pivots.len();
232 let mut pivot = None;
233 for (i, reduced) in reduced.iter_mut().enumerate().skip(r) {
234 let (upper, lower) = self.data.split_at_mut(i);
235 let row = &mut lower[0];
236 for p in *reduced..r {
237 if row.get(pivots[p]) {
238 xor(&mut row.words_mut()[word..], &upper[p].words()[word..]);
239 }
240 }
241 *reduced = r;
242 if row.get(c) {
243 pivot = Some(i);
244 break;
245 }
246 }
247 if let Some(p) = pivot {
248 self.data.swap(r, p);
249 reduced.swap(r, p);
250 pivots.push(c);
251 c += 1;
252 } else if require_full_rank {
253 return pivots;
254 } else {
255 c = self.next_column(r, c + 1, end);
256 }
257 }
258 let rank = pivots.len();
259 if first == rank {
260 let next = self.next_column(rank, start + block, cols);
261 if next == cols {
262 break;
263 }
264 start = next / block * block;
265 continue;
266 }
267 // Make the panel's pivot columns an identity matrix before indexing its table.
268 for r in (first..rank).rev() {
269 let (upper, lower) = self.data.split_at_mut(r);
270 for row in &mut upper[first..] {
271 if row.get(pivots[r]) {
272 xor(&mut row.words_mut()[word..], &lower[0].words()[word..]);
273 }
274 }
275 }
276 let mut indices = [[0usize; 256]; 4];
277 let mut masks = [0usize; 4];
278 for group in 0..block.div_ceil(8) {
279 let mut keys = [0usize; 256];
280 let mut count = 0;
281 for (row, &c) in pivots[first..].iter().enumerate() {
282 if (c - start) / 8 != group {
283 continue;
284 }
285 let half = 1 << count;
286 count += 1;
287 for index in 0..half {
288 keys[index + half] = keys[index] | (1 << ((c - start) % 8));
289 let (lower, upper) = table.split_at_mut(group * 256 + index + half);
290 let target = &mut upper[0].words_mut()[word..];
291 let source = &lower[group * 256 + index].words()[word..];
292 let pivot = &self[first + row].words()[word..];
293 for ((x, y), z) in target.iter_mut().zip(source).zip(pivot) {
294 *x = y ^ z;
295 }
296 }
297 }
298 masks[group] = keys[(1 << count) - 1];
299 for (i, &key) in keys[..1 << count].iter().enumerate() {
300 indices[group][key] = i;
301 }
302 }
303 for i in (rank..n).chain(0..if full { first } else { 0 }) {
304 let key = (self[i].words()[word] >> (start & 63)) as usize;
305 let x = indices[0][key & masks[0]];
306 if block <= 8 {
307 if x != 0 {
308 xor(&mut self[i].words_mut()[word..], &table[x].words()[word..]);
309 }
310 continue;
311 }
312 let y = indices[1][(key >> 8) & masks[1]];
313 let z = indices[2][(key >> 16) & masks[2]];
314 let w = indices[3][(key >> 24) & masks[3]];
315 if x | y | z | w != 0 {
316 let p = &table[x].words()[word..];
317 let q = &table[256 + y].words()[word..];
318 let r = &table[512 + z].words()[word..];
319 let s = &table[768 + w].words()[word..];
320 for ((((x, y), z), r), s) in self[i].words_mut()[word..]
321 .iter_mut()
322 .zip(p)
323 .zip(q)
324 .zip(r)
325 .zip(s)
326 {
327 *x ^= y ^ z ^ r ^ s;
328 }
329 }
330 }
331 if rank == n {
332 break;
333 }
334 start += block;
335 }
336 pivots
337 }
338
339 fn next_column(&self, row: usize, mut start: usize, cols: usize) -> usize {
340 while start < cols {
341 let word = start / 64;
342 let bits = self.data[row..]
343 .iter()
344 .fold(0, |x, row| x | row.words()[word])
345 & (u64::MAX << (start & 63));
346 if bits != 0 {
347 return (word * 64 + bits.trailing_zeros() as usize).min(cols);
348 }
349 start = (word + 1) * 64;
350 }
351 cols
352 }
353
354 #[inline(always)]
355 fn eliminate_sparse(&mut self, cols: usize, full: bool, require_full_rank: bool) -> Vec<usize> {
356 let n = self.shape.0;
357 let mut basis = vec![n; cols];
358 let mut pivots = Vec::new();
359 for i in 0..n {
360 loop {
361 let Some(c) = self[i].iter_ones().next().filter(|&c| c < cols) else {
362 if require_full_rank {
363 return pivots;
364 }
365 break;
366 };
367 if basis[c] == n {
368 basis[c] = i;
369 pivots.push(c);
370 break;
371 }
372 let (upper, lower) = self.data.split_at_mut(i);
373 xor(
374 &mut lower[0].words_mut()[c / 64..],
375 &upper[basis[c]].words()[c / 64..],
376 );
377 }
378 }
379 pivots.sort_unstable();
380 self.data
381 .sort_by_cached_key(|row| row.iter_ones().next().filter(|&c| c < cols).unwrap_or(cols));
382 if full {
383 for (i, &c) in pivots.iter().enumerate() {
384 basis[c] = i;
385 }
386 for i in (0..pivots.len()).rev() {
387 let next = self[i].iter_ones().take_while(|&c| c < cols).nth(1);
388 if let Some(c) = next
389 && basis[c] != n
390 {
391 let (upper, lower) = self.data.split_at_mut(basis[c]);
392 xor(
393 &mut upper[i].words_mut()[c / 64..],
394 &lower[0].words()[c / 64..],
395 );
396 }
397 }
398 }
399 pivots
400 }
401
402 #[inline(always)]
403 fn mul_impl(&self, rhs: &Self) -> Self {
404 let mut result = Self::zeros((self.shape.0, rhs.shape.1));
405 let ones = self.data.iter().map(BitSet::count_ones).sum::<u64>();
406 let size = self.shape.0 as u64 * self.shape.1 as u64;
407 if self.shape.0 < 256 || self.shape.1 < 32 || ones <= size / 8 {
408 for (a, c) in self.data.iter().zip(&mut result.data) {
409 for j in a.iter_ones() {
410 xor(c.words_mut(), rhs[j].words());
411 }
412 }
413 return result;
414 }
415 if size - ones <= size / 8 {
416 let mut sum = BitSet::new(rhs.shape.1);
417 for row in &rhs.data {
418 sum ^= row;
419 }
420 for (a, c) in self.data.iter().zip(&mut result.data) {
421 c.words_mut().copy_from_slice(sum.words());
422 for j in (!a.clone()).iter_ones() {
423 xor(c.words_mut(), rhs[j].words());
424 }
425 }
426 return result;
427 }
428 let width = rhs.shape.1.div_ceil(64);
429 if width == 0 {
430 return result;
431 }
432
433 // Separate the table groups by a cache line to avoid mapping them to the same sets.
434 let group = 256 * width + 8;
435 let mut storage = BitSet::new(8 * group * 64);
436 let table = storage.words_mut();
437 for start in (0..self.shape.1).step_by(64) {
438 for (t, table) in table.chunks_exact_mut(group).enumerate() {
439 let col = start + t * 8;
440 for bit in 0..self.shape.1.saturating_sub(col).min(8) {
441 let row = rhs[col + bit].words();
442 let half = (1 << bit) * width;
443 let (lower, upper) = table.split_at_mut(half);
444 for (source, dest) in
445 lower.chunks_exact(width).zip(upper.chunks_exact_mut(width))
446 {
447 for ((x, y), z) in dest.iter_mut().zip(source).zip(row) {
448 *x = y ^ z;
449 }
450 }
451 }
452 }
453 for (a, c) in self.data.iter().zip(&mut result.data) {
454 let key = a.words()[start / 64];
455 let offset = (key & 255) as usize * width;
456 let p0 = &table[offset..offset + width];
457 let offset = group + (key >> 8 & 255) as usize * width;
458 let p1 = &table[offset..offset + width];
459 let offset = 2 * group + (key >> 16 & 255) as usize * width;
460 let p2 = &table[offset..offset + width];
461 let offset = 3 * group + (key >> 24 & 255) as usize * width;
462 let p3 = &table[offset..offset + width];
463 let offset = 4 * group + (key >> 32 & 255) as usize * width;
464 let p4 = &table[offset..offset + width];
465 let offset = 5 * group + (key >> 40 & 255) as usize * width;
466 let p5 = &table[offset..offset + width];
467 let offset = 6 * group + (key >> 48 & 255) as usize * width;
468 let p6 = &table[offset..offset + width];
469 let offset = 7 * group + (key >> 56 & 255) as usize * width;
470 let p7 = &table[offset..offset + width];
471 for ((((((((x, p0), p1), p2), p3), p4), p5), p6), p7) in c
472 .words_mut()
473 .iter_mut()
474 .zip(p0)
475 .zip(p1)
476 .zip(p2)
477 .zip(p3)
478 .zip(p4)
479 .zip(p5)
480 .zip(p6)
481 .zip(p7)
482 {
483 *x ^= p0 ^ p1 ^ p2 ^ p3 ^ p4 ^ p5 ^ p6 ^ p7;
484 }
485 }
486 }
487 result
488 }