M4RI 20140914
mzd.h
Go to the documentation of this file.
1
10#ifndef M4RI_MZD
11#define M4RI_MZD
12
13/*******************************************************************
14*
15* M4RI: Linear Algebra over GF(2)
16*
17* Copyright (C) 2007, 2008 Gregory Bard <bard@fordham.edu>
18* Copyright (C) 2008-2013 Martin Albrecht <M.R.Albrecht@rhul.ac.uk>
19* Copyright (C) 2011 Carlo Wood <carlo@alinoe.com>
20*
21* Distributed under the terms of the GNU General Public License (GPL)
22* version 2 or higher.
23*
24* This code is distributed in the hope that it will be useful,
25* but WITHOUT ANY WARRANTY; without even the implied warranty of
26* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27* General Public License for more details.
28*
29* The full text of the GPL is available at:
30*
31* http://www.gnu.org/licenses/
32*
33********************************************************************/
34
35#ifdef HAVE_CONFIG_H
36#include "config.h"
37#endif
38
39#include <m4ri/m4ri_config.h>
40
41#include <assert.h>
42#include <math.h>
43#include <stdio.h>
44
45#if __M4RI_HAVE_SSE2
46#include <emmintrin.h>
47#endif
48
49#include <m4ri/misc.h>
50#include <m4ri/debug_dump.h>
51
58#define __M4RI_MAX_MZD_BLOCKSIZE (((size_t)1) << 27)
59
68#define __M4RI_MUL_BLOCKSIZE MIN(((int)sqrt((double)(4 * __M4RI_CPU_L3_CACHE))) / 2, 2048)
69
74typedef struct {
75 size_t size;
79
86typedef struct mzd_t {
87
100
109
124 uint8_t flags;
125
132
133 /* ensures sizeof(mzd_t) == 64 */
134 uint8_t padding[62 - 2 * sizeof(rci_t) - 4 * sizeof(wi_t) - sizeof(word) - 2 * sizeof(void *)];
135
140
144static wi_t const mzd_paddingwidth = 1;
145
150static uint8_t const mzd_flag_nonzero_excess = 0x2;
151
156static uint8_t const mzd_flag_windowed_zerooffset = 0x4;
157
162static uint8_t const mzd_flag_windowed_zeroexcess = 0x8;
163
168static uint8_t const mzd_flag_windowed_ownsblocks = 0x10;
169
174static uint8_t const mzd_flag_multiple_blocks = 0x20;
175
183static inline int mzd_is_windowed(mzd_t const *M) { return M->flags & (mzd_flag_windowed_zerooffset); }
184
192static inline int mzd_owns_blocks(mzd_t const *M) {
193 return M->blocks && (!mzd_is_windowed(M) || ((M->flags & mzd_flag_windowed_ownsblocks)));
194}
195
204static inline word *mzd_first_row(mzd_t const *M) {
205 word *result = M->blocks[0].begin + M->offset_vector;
206 assert(M->nrows == 0 || result == M->rows[0]);
207 return result;
208}
209
220static inline word *mzd_first_row_next_block(mzd_t const *M, int n) {
221 assert(n > 0);
222 return M->blocks[n].begin + M->offset_vector - M->row_offset * M->rowstride;
223}
224
234static inline int mzd_row_to_block(mzd_t const *M, rci_t row) { return (M->row_offset + row) >> M->blockrows_log; }
235
249static inline wi_t mzd_rows_in_block(mzd_t const *M, int n) {
251 if (__M4RI_UNLIKELY(n == 0)) {
252 return (1 << M->blockrows_log) - M->row_offset;
253 } else {
254 int const last_block = mzd_row_to_block(M, M->nrows - 1);
255 if (n < last_block) {
256 return (1 << M->blockrows_log);
257 }
258 return M->nrows + M->row_offset - (n << M->blockrows_log);
259 }
260 }
261 return n ? 0 : M->nrows;
262}
263
273static inline wi_t mzd_remaining_rows_in_block(mzd_t const *M, rci_t r) {
274 const int n = mzd_row_to_block(M, r);
275 r = (r - (n << M->blockrows_log));
277 if (__M4RI_UNLIKELY(n == 0)) {
278 return (1 << M->blockrows_log) - M->row_offset - r;
279 } else {
280 int const last_block = mzd_row_to_block(M, M->nrows - 1);
281 if (n < last_block) {
282 return (1 << M->blockrows_log) - r;
283 }
284 return M->nrows + M->row_offset - (n << M->blockrows_log) - r;
285 }
286 }
287 return n ? 0 : M->nrows - r;
288}
289
299static inline word *mzd_row(mzd_t const *M, rci_t row) {
300 wi_t big_vector = M->offset_vector + row * M->rowstride;
301 word *result = M->blocks[0].begin + big_vector;
303 int const n = (M->row_offset + row) >> M->blockrows_log;
304 result = M->blocks[n].begin + big_vector - n * (M->blocks[0].size / sizeof(word));
305 }
306 assert(result == M->rows[row]);
307 return result;
308}
309
320mzd_t *mzd_init(rci_t const r, rci_t const c);
321
328void mzd_free(mzd_t *A);
329
351mzd_t *mzd_init_window(mzd_t *M, rci_t const lowr, rci_t const lowc, rci_t const highr, rci_t const highc);
352
359static inline mzd_t const *mzd_init_window_const(mzd_t const *M, rci_t const lowr, rci_t const lowc, rci_t const highr,
360 rci_t const highc) {
361 return mzd_init_window((mzd_t *)M, lowr, lowc, highr, highc);
362}
363
370#define mzd_free_window mzd_free
371
381static inline void _mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb, wi_t const startblock) {
382 if ((rowa == rowb) || (startblock >= M->width)) {
383 return;
384 }
385
386 wi_t width = M->width - startblock - 1;
387 word *a = M->rows[rowa] + startblock;
388 word *b = M->rows[rowb] + startblock;
389 word tmp;
390 word const mask_end = M->high_bitmask;
391
392 for (wi_t i = 0; i < width; ++i) {
393 tmp = a[i];
394 a[i] = b[i];
395 b[i] = tmp;
396 }
397 tmp = (a[width] ^ b[width]) & mask_end;
398 a[width] ^= tmp;
399 b[width] ^= tmp;
400
401 __M4RI_DD_ROW(M, rowa);
402 __M4RI_DD_ROW(M, rowb);
403}
404
413static inline void mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb) { _mzd_row_swap(M, rowa, rowb, 0); }
414
427void mzd_copy_row(mzd_t *B, rci_t i, mzd_t const *A, rci_t j);
428
437void mzd_col_swap(mzd_t *M, rci_t const cola, rci_t const colb);
438
449static inline void mzd_col_swap_in_rows(mzd_t *M, rci_t const cola, rci_t const colb, rci_t const start_row,
450 rci_t const stop_row) {
451 if (cola == colb) {
452 return;
453 }
454
455 rci_t const _cola = cola;
456 rci_t const _colb = colb;
457
458 wi_t const a_word = _cola / m4ri_radix;
459 wi_t const b_word = _colb / m4ri_radix;
460
461 int const a_bit = _cola % m4ri_radix;
462 int const b_bit = _colb % m4ri_radix;
463
464 word *RESTRICT ptr = mzd_row(M, start_row);
465 int max_bit = MAX(a_bit, b_bit);
466 int count_remaining = stop_row - start_row;
467 int min_bit = a_bit + b_bit - max_bit;
468 int block = mzd_row_to_block(M, start_row);
469 int offset = max_bit - min_bit;
470 word mask = m4ri_one << min_bit;
471 int count = MIN(mzd_remaining_rows_in_block(M, start_row), count_remaining);
472
473 // Apparently we're calling with start_row == stop_row sometimes (seems a bug to me).
474 if (count <= 0) {
475 return;
476 }
477
478 if (a_word == b_word) {
479 while (1) {
480 count_remaining -= count;
481 ptr += a_word;
482 int fast_count = count / 4;
483 int rest_count = count - 4 * fast_count;
484 word xor_v[4];
485 wi_t const rowstride = M->rowstride;
486 while (fast_count--) {
487 xor_v[0] = ptr[0];
488 xor_v[1] = ptr[rowstride];
489 xor_v[2] = ptr[2 * rowstride];
490 xor_v[3] = ptr[3 * rowstride];
491 xor_v[0] ^= xor_v[0] >> offset;
492 xor_v[1] ^= xor_v[1] >> offset;
493 xor_v[2] ^= xor_v[2] >> offset;
494 xor_v[3] ^= xor_v[3] >> offset;
495 xor_v[0] &= mask;
496 xor_v[1] &= mask;
497 xor_v[2] &= mask;
498 xor_v[3] &= mask;
499 xor_v[0] |= xor_v[0] << offset;
500 xor_v[1] |= xor_v[1] << offset;
501 xor_v[2] |= xor_v[2] << offset;
502 xor_v[3] |= xor_v[3] << offset;
503 ptr[0] ^= xor_v[0];
504 ptr[rowstride] ^= xor_v[1];
505 ptr[2 * rowstride] ^= xor_v[2];
506 ptr[3 * rowstride] ^= xor_v[3];
507 ptr += 4 * rowstride;
508 }
509 while (rest_count--) {
510 word xor_v = *ptr;
511 xor_v ^= xor_v >> offset;
512 xor_v &= mask;
513 *ptr ^= xor_v | (xor_v << offset);
514 ptr += rowstride;
515 }
516 block++;
517 if ((count = MIN(mzd_rows_in_block(M, block), count_remaining)) <= 0) {
518 break;
519 }
520 ptr = mzd_first_row_next_block(M, block);
521 }
522 } else {
523 word *RESTRICT min_ptr;
524 wi_t max_offset;
525 if (min_bit == a_bit) {
526 min_ptr = ptr + a_word;
527 max_offset = b_word - a_word;
528 } else {
529 min_ptr = ptr + b_word;
530 max_offset = a_word - b_word;
531 }
532 while (1) {
533 count_remaining -= count;
534 wi_t const rowstride = M->rowstride;
535 while (count--) {
536 word xor_v = (min_ptr[0] ^ (min_ptr[max_offset] >> offset)) & mask;
537 min_ptr[0] ^= xor_v;
538 min_ptr[max_offset] ^= xor_v << offset;
539 min_ptr += rowstride;
540 }
541 block++;
542 if ((count = MIN(mzd_rows_in_block(M, +block), count_remaining)) <= 0) {
543 break;
544 }
545 ptr = mzd_first_row_next_block(M, block);
546 if (min_bit == a_bit) {
547 min_ptr = ptr + a_word;
548 } else {
549 min_ptr = ptr + b_word;
550 }
551 }
552 }
553
554 __M4RI_DD_MZD(M);
555}
556
568static inline BIT mzd_read_bit(mzd_t const *M, rci_t const row, rci_t const col) {
569 return __M4RI_GET_BIT(M->rows[row][col / m4ri_radix], col % m4ri_radix);
570}
571
584static inline void mzd_write_bit(mzd_t *M, rci_t const row, rci_t const col, BIT const value) {
585 __M4RI_WRITE_BIT(M->rows[row][col / m4ri_radix], col % m4ri_radix, value);
586}
587
598static inline void mzd_xor_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values) {
599 int const spot = y % m4ri_radix;
600 wi_t const block = y / m4ri_radix;
601 M->rows[x][block] ^= values << spot;
602 int const space = m4ri_radix - spot;
603 if (n > space) {
604 M->rows[x][block + 1] ^= values >> space;
605 }
606}
607
618static inline void mzd_and_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values) {
619 /* This is the best way, since this will drop out once we inverse the bits in values: */
620 values >>= (m4ri_radix - n); /* Move the bits to the lowest columns */
621
622 int const spot = y % m4ri_radix;
623 wi_t const block = y / m4ri_radix;
624 M->rows[x][block] &= values << spot;
625 int const space = m4ri_radix - spot;
626 if (n > space) {
627 M->rows[x][block + 1] &= values >> space;
628 }
629}
630
640static inline void mzd_clear_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
641 assert(n > 0 && n <= m4ri_radix);
642 word values = m4ri_ffff >> (m4ri_radix - n);
643 int const spot = y % m4ri_radix;
644 wi_t const block = y / m4ri_radix;
645 M->rows[x][block] &= ~(values << spot);
646 int const space = m4ri_radix - spot;
647 if (n > space) {
648 M->rows[x][block + 1] &= ~(values >> space);
649 }
650}
651
664static inline void mzd_row_add_offset(mzd_t *M, rci_t dstrow, rci_t srcrow, rci_t coloffset) {
665 assert(dstrow < M->nrows && srcrow < M->nrows && coloffset < M->ncols);
666 wi_t const startblock = coloffset / m4ri_radix;
667 wi_t wide = M->width - startblock;
668 word *src = M->rows[srcrow] + startblock;
669 word *dst = M->rows[dstrow] + startblock;
670 word const mask_begin = __M4RI_RIGHT_BITMASK(m4ri_radix - coloffset % m4ri_radix);
671 word const mask_end = M->high_bitmask;
672
673 *dst++ ^= *src++ & mask_begin;
674 --wide;
675
676#if __M4RI_HAVE_SSE2
677 int not_aligned = __M4RI_ALIGNMENT(src, 16) != 0; /* 0: Aligned, 1: Not aligned */
678 if (wide > not_aligned + 1) /* Speed up for small matrices */
679 {
680 if (not_aligned) {
681 *dst++ ^= *src++;
682 --wide;
683 }
684 /* Now wide > 1 */
685 __m128i *__src = (__m128i *)src;
686 __m128i *__dst = (__m128i *)dst;
687 __m128i *const eof = (__m128i *)((unsigned long)(src + wide) & ~0xFUL);
688 do {
689 __m128i xmm1 = _mm_xor_si128(*__dst, *__src);
690 *__dst++ = xmm1;
691 } while (++__src < eof);
692 src = (word *)__src;
693 dst = (word *)__dst;
694 wide = ((sizeof(word) * wide) % 16) / sizeof(word);
695 }
696#endif
697 wi_t i = -1;
698 while (++i < wide) {
699 dst[i] ^= src[i];
700 }
701 /*
702 * Revert possibly non-zero excess bits.
703 * Note that i == wide here, and wide can be 0.
704 * But really, src[wide - 1] is M->rows[srcrow][M->width - 1] ;)
705 * We use i - 1 here to let the compiler know these are the same addresses
706 * that we last accessed, in the previous loop.
707 */
708 dst[i - 1] ^= src[i - 1] & ~mask_end;
709
710 __M4RI_DD_ROW(M, dstrow);
711}
712
724void mzd_row_add(mzd_t *M, rci_t const sourcerow, rci_t const destrow);
725
740mzd_t *mzd_transpose(mzd_t *DST, mzd_t const *A);
741
756mzd_t *mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);
757
772mzd_t *mzd_addmul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);
773
785mzd_t *_mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B, int const clear);
786
796mzd_t *_mzd_mul_va(mzd_t *C, mzd_t const *v, mzd_t const *A, int const clear);
797
804void mzd_randomize(mzd_t *M);
805
814typedef word (*m4ri_random_callback)(void* data);
815
823void mzd_randomize_custom(mzd_t *M, m4ri_random_callback rc, void* data);
824
839void mzd_set_ui(mzd_t *M, unsigned int const value);
840
854rci_t mzd_gauss_delayed(mzd_t *M, rci_t const startcol, int const full);
855
869rci_t mzd_echelonize_naive(mzd_t *M, int const full);
870
878int mzd_equal(mzd_t const *A, mzd_t const *B);
879
891int mzd_cmp(mzd_t const *A, mzd_t const *B);
892
900mzd_t *mzd_copy(mzd_t *DST, mzd_t const *A);
901
920mzd_t *mzd_concat(mzd_t *C, mzd_t const *A, mzd_t const *B);
921
939mzd_t *mzd_stack(mzd_t *C, mzd_t const *A, mzd_t const *B);
940
953mzd_t *mzd_submatrix(mzd_t *S, mzd_t const *M, rci_t const lowr, rci_t const lowc, rci_t const highr,
954 rci_t const highc);
955
967mzd_t *mzd_invert_naive(mzd_t *INV, mzd_t const *A, mzd_t const *I);
968
980mzd_t *mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);
981
990mzd_t *_mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);
991
1000#define mzd_sub mzd_add
1001
1010#define _mzd_sub _mzd_add
1011
1021static inline word mzd_read_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
1022 int const spot = y % m4ri_radix;
1023 wi_t const block = y / m4ri_radix;
1024 int const spill = spot + n - m4ri_radix;
1025 word temp = (spill <= 0) ? M->rows[x][block] << -spill
1026 : (M->rows[x][block + 1] << (m4ri_radix - spill)) | (M->rows[x][block] >> spill);
1027 return temp >> (m4ri_radix - n);
1028}
1029
1046static inline void mzd_combine_even_in_place(mzd_t *A, rci_t const a_row, wi_t const a_startblock,
1047 mzd_t const *B, rci_t const b_row, wi_t const b_startblock) {
1048
1049 wi_t wide = A->width - a_startblock - 1;
1050
1051 word *a = A->rows[a_row] + a_startblock;
1052 word *b = B->rows[b_row] + b_startblock;
1053
1054#if __M4RI_HAVE_SSE2
1055 if (wide > 2) {
1057 if (__M4RI_ALIGNMENT(a, 16)) {
1058 *a++ ^= *b++;
1059 wide--;
1060 }
1061
1062 if (__M4RI_ALIGNMENT(a, 16) == 0 && __M4RI_ALIGNMENT(b, 16) == 0) {
1063 __m128i *a128 = (__m128i *)a;
1064 __m128i *b128 = (__m128i *)b;
1065 const __m128i *eof = (__m128i *)((unsigned long)(a + wide) & ~0xFUL);
1066
1067 do {
1068 *a128 = _mm_xor_si128(*a128, *b128);
1069 ++b128;
1070 ++a128;
1071 } while (a128 < eof);
1072
1073 a = (word *)a128;
1074 b = (word *)b128;
1075 wide = ((sizeof(word) * wide) % 16) / sizeof(word);
1076 }
1077 }
1078#endif // __M4RI_HAVE_SSE2
1079
1080 if (wide > 0) {
1081 wi_t n = (wide + 7) / 8;
1082 switch (wide % 8) {
1083 case 0: do { *(a++) ^= *(b++);
1084 case 7: *(a++) ^= *(b++);
1085 case 6: *(a++) ^= *(b++);
1086 case 5: *(a++) ^= *(b++);
1087 case 4: *(a++) ^= *(b++);
1088 case 3: *(a++) ^= *(b++);
1089 case 2: *(a++) ^= *(b++);
1090 case 1: *(a++) ^= *(b++);
1091 } while (--n > 0);
1092 }
1093 }
1094
1095 *a ^= *b & A->high_bitmask;
1096
1097 __M4RI_DD_MZD(A);
1098}
1099
1119static inline void mzd_combine_even(mzd_t *C, rci_t const c_row, wi_t const c_startblock,
1120 mzd_t const *A, rci_t const a_row, wi_t const a_startblock,
1121 mzd_t const *B, rci_t const b_row, wi_t const b_startblock) {
1122
1123 wi_t wide = A->width - a_startblock - 1;
1124 word *a = A->rows[a_row] + a_startblock;
1125 word *b = B->rows[b_row] + b_startblock;
1126 word *c = C->rows[c_row] + c_startblock;
1127
1128#if __M4RI_HAVE_SSE2
1129 if (wide > 2) {
1131 if (__M4RI_ALIGNMENT(a, 16)) {
1132 *c++ = *b++ ^ *a++;
1133 wide--;
1134 }
1135
1136 if ((__M4RI_ALIGNMENT(b, 16) | __M4RI_ALIGNMENT(c, 16)) == 0) {
1137 __m128i *a128 = (__m128i *)a;
1138 __m128i *b128 = (__m128i *)b;
1139 __m128i *c128 = (__m128i *)c;
1140 const __m128i *eof = (__m128i *)((unsigned long)(a + wide) & ~0xFUL);
1141
1142 do {
1143 *c128 = _mm_xor_si128(*a128, *b128);
1144 ++c128;
1145 ++b128;
1146 ++a128;
1147 } while (a128 < eof);
1148
1149 a = (word *)a128;
1150 b = (word *)b128;
1151 c = (word *)c128;
1152 wide = ((sizeof(word) * wide) % 16) / sizeof(word);
1153 }
1154 }
1155#endif // __M4RI_HAVE_SSE2
1156
1157 if (wide > 0) {
1158 wi_t n = (wide + 7) / 8;
1159 switch (wide % 8) {
1160 case 0: do { *(c++) = *(a++) ^ *(b++);
1161 case 7: *(c++) = *(a++) ^ *(b++);
1162 case 6: *(c++) = *(a++) ^ *(b++);
1163 case 5: *(c++) = *(a++) ^ *(b++);
1164 case 4: *(c++) = *(a++) ^ *(b++);
1165 case 3: *(c++) = *(a++) ^ *(b++);
1166 case 2: *(c++) = *(a++) ^ *(b++);
1167 case 1: *(c++) = *(a++) ^ *(b++);
1168 } while (--n > 0);
1169 }
1170 }
1171 *c ^= ((*a ^ *b ^ *c) & C->high_bitmask);
1172
1173 __M4RI_DD_MZD(C);
1174}
1175
1194static inline void mzd_combine(mzd_t *C, rci_t const c_row, wi_t const c_startblock,
1195 mzd_t const *A, rci_t const a_row, wi_t const a_startblock,
1196 mzd_t const *B, rci_t const b_row, wi_t const b_startblock) {
1197
1198 if ((C == A) & (a_row == c_row) & (a_startblock == c_startblock)) {
1199 mzd_combine_even_in_place(C, c_row, c_startblock, B, b_row, b_startblock);
1200 } else {
1201 mzd_combine_even(C, c_row, c_startblock, A, a_row, a_startblock, B, b_row, b_startblock);
1202 }
1203 return;
1204}
1205
1214static inline int mzd_read_bits_int(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
1215 return __M4RI_CONVERT_TO_INT(mzd_read_bits(M, x, y, n));
1216}
1217
1224int mzd_is_zero(mzd_t const *A);
1225
1234void mzd_row_clear_offset(mzd_t *M, rci_t const row, rci_t const coloffset);
1235
1252int mzd_find_pivot(mzd_t const *M, rci_t start_row, rci_t start_col, rci_t *r, rci_t *c);
1253
1266double mzd_density(mzd_t const *A, wi_t res);
1267
1282double _mzd_density(mzd_t const *A, wi_t res, rci_t r, rci_t c);
1283
1293
1300static inline word mzd_hash(mzd_t const *A) {
1301 word hash = 0;
1302 for (rci_t r = 0; r < A->nrows; ++r) {
1303 hash ^= rotate_word(calculate_hash(A->rows[r], A->width), r % m4ri_radix);
1304 }
1305 return hash;
1306}
1307
1317mzd_t *mzd_extract_u(mzd_t *U, mzd_t const *A);
1318
1328mzd_t *mzd_extract_l(mzd_t *L, mzd_t const *A);
1329
1330#endif // M4RI_MZD
Helper functions.
int rci_t
Type of row and column indexes.
Definition: misc.h:72
#define __M4RI_ALIGNMENT(addr, n)
Return alignment of addr w.r.t. n. For example the address 17 would be 1 aligned w....
Definition: misc.h:421
#define __M4RI_UNLIKELY(cond)
Macro to help with branch prediction.
Definition: misc.h:449
#define __M4RI_CONVERT_TO_INT(w)
Explicit conversion macro.
Definition: misc.h:99
#define __M4RI_GET_BIT(w, spot)
Get the bit spot (counting from the left) in the word w.
Definition: misc.h:226
#define __M4RI_WRITE_BIT(w, spot, value)
Write the value to the bit spot in the word w.
Definition: misc.h:236
uint64_t word
A word is the typical packed data structure to represent packed bits.
Definition: misc.h:87
int BIT
Pretty for a boolean int.
Definition: misc.h:64
#define MIN(x, y)
Return the minimal element of x and y.
Definition: misc.h:174
static word const m4ri_ffff
A word with all bits set.
Definition: misc.h:153
static int const m4ri_radix
The number of bits in a word.
Definition: misc.h:141
#define MAX(x, y)
Return the maximal element of x and y.
Definition: misc.h:163
int wi_t
Type of word indexes.
Definition: misc.h:80
#define __M4RI_RIGHT_BITMASK(n)
create a bit mask to zero out all but the n rightmost bits.
Definition: misc.h:296
static word const m4ri_one
The number one as a word.
Definition: misc.h:147
mzd_t * mzd_transpose(mzd_t *DST, mzd_t const *A)
Transpose a matrix.
Definition: mzd.c:1389
static word * mzd_first_row_next_block(mzd_t const *M, int n)
Get a pointer to the first word in block n.
Definition: mzd.h:220
static int mzd_is_windowed(mzd_t const *M)
Test if a matrix is windowed.
Definition: mzd.h:183
static uint8_t const mzd_flag_nonzero_excess
flag when ncols%64 == 0
Definition: mzd.h:150
mzd_t * mzd_submatrix(mzd_t *S, mzd_t const *M, rci_t const lowr, rci_t const lowc, rci_t const highr, rci_t const highc)
Copy a submatrix.
Definition: mzd.c:1867
static uint8_t const mzd_flag_windowed_zerooffset
flag for windowed matrix
Definition: mzd.h:156
static uint8_t const mzd_flag_windowed_zeroexcess
flag for windowed matrix where ncols%64 == 0
Definition: mzd.h:162
rci_t mzd_echelonize_naive(mzd_t *M, int const full)
Gaussian elimination.
Definition: mzd.c:335
word(* m4ri_random_callback)(void *data)
Random callback that produces uniformly distributed random words on every call.
Definition: mzd.h:814
static uint8_t const mzd_flag_multiple_blocks
flag for multiply blocks
Definition: mzd.h:174
static void mzd_col_swap_in_rows(mzd_t *M, rci_t const cola, rci_t const colb, rci_t const start_row, rci_t const stop_row)
Swap the two columns cola and colb but only between start_row and stop_row.
Definition: mzd.h:449
void mzd_randomize(mzd_t *M)
Fill matrix M with uniformly distributed bits.
Definition: mzd.c:1558
void mzd_row_clear_offset(mzd_t *M, rci_t const row, rci_t const coloffset)
Clear the given row, but only begins at the column coloffset.
Definition: mzd.c:288
void mzd_randomize_custom(mzd_t *M, m4ri_random_callback rc, void *data)
Fill matrix M with uniformly distributed bits.
Definition: mzd.c:1570
mzd_t * mzd_stack(mzd_t *C, mzd_t const *A, mzd_t const *B)
Stack A on top of B and write the result to C.
Definition: mzd.c:1709
void mzd_row_add(mzd_t *M, rci_t const sourcerow, rci_t const destrow)
Add the rows sourcerow and destrow and stores the total in the row destrow.
Definition: mzd.c:284
void mzd_free(mzd_t *A)
Free a matrix created with mzd_init.
Definition: mzd.c:271
static word * mzd_row(mzd_t const *M, rci_t row)
Get pointer to first word of row.
Definition: mzd.h:299
rci_t mzd_first_zero_row(mzd_t const *A)
Return the first row with all zero entries.
Definition: mzd.c:2211
static wi_t mzd_remaining_rows_in_block(mzd_t const *M, rci_t r)
Number of rows in this block including r.
Definition: mzd.h:273
static void mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb)
Swap the two rows rowa and rowb.
Definition: mzd.h:413
mzd_t * mzd_init(rci_t const r, rci_t const c)
Create a new matrix of dimension r x c.
Definition: mzd.c:149
static wi_t mzd_rows_in_block(mzd_t const *M, int n)
Total number of rows in this block.
Definition: mzd.h:249
int mzd_cmp(mzd_t const *A, mzd_t const *B)
Return -1,0,1 if if A < B, A == B or A > B respectively.
Definition: mzd.c:1627
static void mzd_row_add_offset(mzd_t *M, rci_t dstrow, rci_t srcrow, rci_t coloffset)
Add the rows sourcerow and destrow and stores the total in the row destrow, but only begins at the co...
Definition: mzd.h:664
void mzd_copy_row(mzd_t *B, rci_t i, mzd_t const *A, rci_t j)
copy row j from A to row i from B.
Definition: mzd.c:2014
mzd_t * mzd_invert_naive(mzd_t *INV, mzd_t const *A, mzd_t const *I)
Invert the matrix target using Gaussian elimination.
Definition: mzd.c:1740
mzd_t * _mzd_mul_va(mzd_t *C, mzd_t const *v, mzd_t const *A, int const clear)
Matrix multiplication optimized for v*A where v is a vector.
Definition: mzd.c:1542
static word * mzd_first_row(mzd_t const *M)
Get a pointer the first word.
Definition: mzd.h:204
mzd_t * mzd_addmul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B)
Naive cubic matrix multiplication and addition.
Definition: mzd.c:1438
int mzd_is_zero(mzd_t const *A)
Zero test for matrix.
Definition: mzd.c:2001
double mzd_density(mzd_t const *A, wi_t res)
Return the number of nonzero entries divided by nrows * ncols.
Definition: mzd.c:2207
static int mzd_owns_blocks(mzd_t const *M)
Test if this mzd_t should free blocks.
Definition: mzd.h:192
static uint8_t const mzd_flag_windowed_ownsblocks
flag for windowed matrix wich owns its memory
Definition: mzd.h:168
static void mzd_combine_even(mzd_t *C, rci_t const c_row, wi_t const c_startblock, mzd_t const *A, rci_t const a_row, wi_t const a_startblock, mzd_t const *B, rci_t const b_row, wi_t const b_startblock)
c_row[c_startblock:] = a_row[a_startblock:] + b_row[b_startblock:] for offset 0
Definition: mzd.h:1119
static void _mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb, wi_t const startblock)
Swap the two rows rowa and rowb starting at startblock.
Definition: mzd.h:381
static void mzd_combine_even_in_place(mzd_t *A, rci_t const a_row, wi_t const a_startblock, mzd_t const *B, rci_t const b_row, wi_t const b_startblock)
a_row[a_startblock:] += b_row[b_startblock:] for offset 0
Definition: mzd.h:1046
static word mzd_hash(mzd_t const *A)
Return hash value for matrix.
Definition: mzd.h:1300
static void mzd_clear_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n)
Clear n bits in M starting a position (x,y).
Definition: mzd.h:640
mzd_t * mzd_concat(mzd_t *C, mzd_t const *A, mzd_t const *B)
Concatenate B to A and write the result to C.
Definition: mzd.c:1680
static int mzd_read_bits_int(mzd_t const *M, rci_t const x, rci_t const y, int const n)
Get n bits starting a position (x,y) from the matrix M.
Definition: mzd.h:1214
static int mzd_row_to_block(mzd_t const *M, rci_t row)
Convert row to blocks index.
Definition: mzd.h:234
static void mzd_xor_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values)
XOR n bits from values to M starting a position (x,y).
Definition: mzd.h:598
rci_t mzd_gauss_delayed(mzd_t *M, rci_t const startcol, int const full)
Gaussian elimination.
Definition: mzd.c:308
static void mzd_combine(mzd_t *C, rci_t const c_row, wi_t const c_startblock, mzd_t const *A, rci_t const a_row, wi_t const a_startblock, mzd_t const *B, rci_t const b_row, wi_t const b_startblock)
row3[col3:] = row1[col1:] + row2[col2:]
Definition: mzd.h:1194
double _mzd_density(mzd_t const *A, wi_t res, rci_t r, rci_t c)
Return the number of nonzero entries divided by nrows * ncols considering only the submatrix starting...
Definition: mzd.c:2170
static word mzd_read_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n)
Definition: mzd.h:1021
static void mzd_write_bit(mzd_t *M, rci_t const row, rci_t const col, BIT const value)
Write the bit value to position M[row,col].
Definition: mzd.h:584
mzd_t * mzd_extract_u(mzd_t *U, mzd_t const *A)
Definition: mzd.c:2231
void mzd_col_swap(mzd_t *M, rci_t const cola, rci_t const colb)
Swap the two columns cola and colb.
Definition: mzd.c:1907
mzd_t * _mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B, int const clear)
Naive cubic matrix multiplication with the pre-transposed B.
Definition: mzd.c:1453
mzd_t * mzd_init_window(mzd_t *M, rci_t const lowr, rci_t const lowc, rci_t const highr, rci_t const highc)
Create a window/view into the matrix M.
Definition: mzd.c:229
int mzd_equal(mzd_t const *A, mzd_t const *B)
Return TRUE if A == B.
Definition: mzd.c:1605
static wi_t const mzd_paddingwidth
The minimum width where padding occurs.
Definition: mzd.h:144
static BIT mzd_read_bit(mzd_t const *M, rci_t const row, rci_t const col)
Read the bit at position M[row,col].
Definition: mzd.h:568
mzd_t * mzd_copy(mzd_t *DST, mzd_t const *A)
Copy matrix A to DST.
Definition: mzd.c:1655
mzd_t * mzd_extract_l(mzd_t *L, mzd_t const *A)
Definition: mzd.c:2247
static void mzd_and_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values)
AND n bits from values to M starting a position (x,y).
Definition: mzd.h:618
mzd_t * mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B)
Set C = A+B.
Definition: mzd.c:1760
int mzd_find_pivot(mzd_t const *M, rci_t start_row, rci_t start_col, rci_t *r, rci_t *c)
Find the next nonzero entry in M starting at start_row and start_col.
Definition: mzd.c:2036
static mzd_t const * mzd_init_window_const(mzd_t const *M, rci_t const lowr, rci_t const lowc, rci_t const highr, rci_t const highc)
Create a const window/view into a const matrix M.
Definition: mzd.h:359
mzd_t * _mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B)
Same as mzd_add but without any checks on the input.
Definition: mzd.c:1774
void mzd_set_ui(mzd_t *M, unsigned int const value)
Set the matrix M to the value equivalent to the integer value provided.
Definition: mzd.c:1582
mzd_t * mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B)
Naive cubic matrix multiplication.
Definition: mzd.c:1420
Data containers containing the values packed into words.
Definition: mzd.h:74
word * begin
Definition: mzd.h:76
size_t size
Definition: mzd.h:75
word * end
Definition: mzd.h:77
Dense matrices over GF(2).
Definition: mzd.h:86
wi_t rowstride
Definition: mzd.h:99
wi_t width
Definition: mzd.h:90
wi_t row_offset
Definition: mzd.h:110
rci_t nrows
Definition: mzd.h:88
uint8_t blockrows_log
Definition: mzd.h:131
uint8_t flags
Definition: mzd.h:124
wi_t offset_vector
Definition: mzd.h:108
word high_bitmask
Definition: mzd.h:136
word ** rows
Definition: mzd.h:138
rci_t ncols
Definition: mzd.h:89
mzd_block_t * blocks
Definition: mzd.h:137