39 #include <m4ri/m4ri_config.h>
46 #include <emmintrin.h>
49 #include <m4ri/misc.h>
50 #include <m4ri/debug_dump.h>
58 #define __M4RI_MAX_MZD_BLOCKSIZE (((size_t)1) << 27)
68 #define __M4RI_MUL_BLOCKSIZE MIN(((int)sqrt((double)(4 * __M4RI_CPU_L3_CACHE))) / 2, 2048)
134 uint8_t padding[62 - 2*
sizeof(
rci_t) - 4*
sizeof(
wi_t) -
sizeof(
word) - 2*
sizeof(
void*)];
208 assert(M->
nrows == 0 || result == M->
rows[0]);
264 return n ? 0 : M->
nrows;
289 return n ? 0 : M->
nrows - r;
308 assert(result == M->
rows[row]);
373 #define mzd_free_window mzd_free
385 if ((rowa == rowb) || (startblock >= M->
width))
389 word *a = M->
rows[rowa] + startblock;
390 word *b = M->
rows[rowb] + startblock;
394 for(
wi_t i = 0; i < width; ++i) {
399 tmp = (a[width] ^ b[width]) & mask_end;
403 __M4RI_DD_ROW(M, rowa);
404 __M4RI_DD_ROW(M, rowb);
457 rci_t const _cola = cola;
458 rci_t const _colb = colb;
467 int max_bit =
MAX(a_bit, b_bit);
468 int count_remaining = stop_row - start_row;
469 int min_bit = a_bit + b_bit - max_bit;
471 int offset = max_bit - min_bit;
479 if (a_word == b_word) {
481 count_remaining -= count;
483 int fast_count = count / 4;
484 int rest_count = count - 4 * fast_count;
487 while (fast_count--) {
489 xor_v[1] = ptr[rowstride];
490 xor_v[2] = ptr[2 * rowstride];
491 xor_v[3] = ptr[3 * rowstride];
492 xor_v[0] ^= xor_v[0] >> offset;
493 xor_v[1] ^= xor_v[1] >> offset;
494 xor_v[2] ^= xor_v[2] >> offset;
495 xor_v[3] ^= xor_v[3] >> offset;
500 xor_v[0] |= xor_v[0] << offset;
501 xor_v[1] |= xor_v[1] << offset;
502 xor_v[2] |= xor_v[2] << offset;
503 xor_v[3] |= xor_v[3] << offset;
505 ptr[rowstride] ^= xor_v[1];
506 ptr[2 * rowstride] ^= xor_v[2];
507 ptr[3 * rowstride] ^= xor_v[3];
508 ptr += 4 * rowstride;
510 while (rest_count--) {
512 xor_v ^= xor_v >> offset;
514 *ptr ^= xor_v | (xor_v << offset);
523 word* RESTRICT min_ptr;
525 if (min_bit == a_bit) {
526 min_ptr = ptr + a_word;
527 max_offset = b_word - a_word;
529 min_ptr = ptr + b_word;
530 max_offset = a_word - b_word;
533 count_remaining -= count;
536 word xor_v = (min_ptr[0] ^ (min_ptr[max_offset] >> offset)) & mask;
538 min_ptr[max_offset] ^= xor_v << offset;
539 min_ptr += rowstride;
545 if (min_bit == a_bit)
546 min_ptr = ptr + a_word;
548 min_ptr = ptr + b_word;
600 M->
rows[x][block] ^= values << spot;
603 M->
rows[x][block + 1] ^= values >> space;
622 M->
rows[x][block] &= values << spot;
625 M->
rows[x][block + 1] &= values >> space;
642 M->
rows[x][block] &= ~(values << spot);
645 M->
rows[x][block + 1] &= ~(values >> space);
661 assert(dstrow < M->nrows && srcrow < M->nrows && coloffset < M->ncols);
664 word *src = M->
rows[srcrow] + startblock;
665 word *dst = M->
rows[dstrow] + startblock;
669 *dst++ ^= *src++ & mask_begin;
674 if (wide > not_aligned + 1)
681 __m128i* __src = (__m128i*)src;
682 __m128i* __dst = (__m128i*)dst;
683 __m128i*
const eof = (__m128i*)((
unsigned long)(src + wide) & ~0xFUL);
686 __m128i xmm1 = _mm_xor_si128(*__dst, *__src);
689 while(++__src < eof);
692 wide = ((
sizeof(
word)*wide)%16)/
sizeof(
word);
705 dst[i - 1] ^= src[i - 1] & ~mask_end;
707 __M4RI_DD_ROW(M, dstrow);
979 #define mzd_sub mzd_add
989 #define _mzd_sub _mzd_add
1006 word temp = (spill <= 0) ? M->
rows[x][block] << -spill : (M->
rows[x][block + 1] << (m4ri_radix - spill)) | (M->
rows[x][block] >> spill);
1007 return temp >> (m4ri_radix - n);
1030 wi_t wide = A->
width - a_startblock - 1;
1032 word *a = A->
rows[a_row] + a_startblock;
1033 word *b = B->
rows[b_row] + b_startblock;
1035 #if __M4RI_HAVE_SSE2
1044 __m128i *a128 = (__m128i*)a;
1045 __m128i *b128 = (__m128i*)b;
1046 const __m128i *eof = (__m128i*)((
unsigned long)(a + wide) & ~0xFUL);
1049 *a128 = _mm_xor_si128(*a128, *b128);
1052 }
while(a128 < eof);
1056 wide = ((
sizeof(
word) * wide) % 16) /
sizeof(
word);
1059 #endif // __M4RI_HAVE_SSE2
1062 wi_t n = (wide + 7) / 8;
1064 case 0:
do { *(a++) ^= *(b++);
1065 case 7: *(a++) ^= *(b++);
1066 case 6: *(a++) ^= *(b++);
1067 case 5: *(a++) ^= *(b++);
1068 case 4: *(a++) ^= *(b++);
1069 case 3: *(a++) ^= *(b++);
1070 case 2: *(a++) ^= *(b++);
1071 case 1: *(a++) ^= *(b++);
1105 wi_t wide = A->
width - a_startblock - 1;
1106 word *a = A->
rows[a_row] + a_startblock;
1107 word *b = B->
rows[b_row] + b_startblock;
1108 word *c = C->
rows[c_row] + c_startblock;
1110 #if __M4RI_HAVE_SSE2
1119 __m128i *a128 = (__m128i*)a;
1120 __m128i *b128 = (__m128i*)b;
1121 __m128i *c128 = (__m128i*)c;
1122 const __m128i *eof = (__m128i*)((
unsigned long)(a + wide) & ~0xFUL);
1125 *c128 = _mm_xor_si128(*a128, *b128);
1129 }
while(a128 < eof);
1134 wide = ((
sizeof(
word) * wide) % 16) /
sizeof(
word);
1137 #endif // __M4RI_HAVE_SSE2
1140 wi_t n = (wide + 7) / 8;
1142 case 0:
do { *(c++) = *(a++) ^ *(b++);
1143 case 7: *(c++) = *(a++) ^ *(b++);
1144 case 6: *(c++) = *(a++) ^ *(b++);
1145 case 5: *(c++) = *(a++) ^ *(b++);
1146 case 4: *(c++) = *(a++) ^ *(b++);
1147 case 3: *(c++) = *(a++) ^ *(b++);
1148 case 2: *(c++) = *(a++) ^ *(b++);
1149 case 1: *(c++) = *(a++) ^ *(b++);
1181 if( (C == A) & (a_row == c_row) & (a_startblock == c_startblock) )
1184 mzd_combine_even(C, c_row, c_startblock, A, a_row, a_startblock, B, b_row, b_startblock);
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:637
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:1998
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:1693
#define __M4RI_WRITE_BIT(w, spot, value)
Write the value to the bit spot in the word w.
Definition: misc.h:236
static word mzd_hash(mzd_t const *A)
Return hash value for matrix.
Definition: mzd.h:1285
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:1177
rci_t mzd_gauss_delayed(mzd_t *M, rci_t const startcol, int const full)
Gaussian elimination.
Definition: mzd.c:308
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:222
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:1724
static word mzd_read_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n)
Definition: mzd.h:1002
word high_bitmask
Definition: mzd.h:136
struct mzd_t mzd_t
Dense matrices over GF(2).
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:1434
mzd_block_t * blocks
Definition: mzd.h:137
int BIT
Pretty for a boolean int.
Definition: misc.h:64
#define __M4RI_RIGHT_BITMASK(n)
create a bit mask to zero out all but the n rightmost bits.
Definition: misc.h:296
static int const m4ri_radix
The number of bits in a word.
Definition: misc.h:141
mzd_t * mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B)
Naive cubic matrix multiplication.
Definition: mzd.c:1416
uint8_t flags
Definition: mzd.h:124
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:362
static uint8_t const mzd_flag_windowed_zeroexcess
flag for windowed matrix where ncols%64 == 0
Definition: mzd.h:162
Dense matrices over GF(2).
Definition: mzd.h:86
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:1027
int rci_t
Type of row and column indexes.
Definition: misc.h:72
#define MIN(x, y)
Return the minimal element of x and y.
Definition: misc.h:174
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:1851
int mzd_equal(mzd_t const *A, mzd_t const *B)
Return TRUE if A == B.
Definition: mzd.c:1589
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:1538
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:616
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:1664
void mzd_free(mzd_t *A)
Free a matrix created with mzd_init.
Definition: mzd.c:271
static word const m4ri_ffff
A word with all bits set.
Definition: misc.h:153
mzd_t * mzd_transpose(mzd_t *DST, mzd_t const *A)
Transpose a matrix.
Definition: mzd.c:1385
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:582
mzd_t * mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B)
Set C = A+B.
Definition: mzd.c:1744
mzd_t * mzd_extract_l(mzd_t *L, mzd_t const *A)
Definition: mzd.c:2231
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:1449
rci_t mzd_echelonize_naive(mzd_t *M, int const full)
Gaussian elimination.
Definition: mzd.c:335
rci_t nrows
Definition: mzd.h:88
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:453
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:384
size_t size
Definition: mzd.h:75
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:276
word ** rows
Definition: mzd.h:138
static uint8_t const mzd_flag_windowed_zerooffset
flag for windowed matrix
Definition: mzd.h:156
#define MAX(x, y)
Return the maximal element of x and y.
Definition: misc.h:163
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
void mzd_randomize(mzd_t *M)
Fill matrix M with uniformly distributed bits.
Definition: mzd.c:1554
#define __M4RI_ALIGNMENT(addr, n)
Return alignment of addr w.r.t. n. For example the address 17 would be 1 aligned w.r.t. 16.
Definition: misc.h:421
mzd_t * mzd_copy(mzd_t *DST, mzd_t const *A)
Copy matrix A to DST.
Definition: mzd.c:1639
Data containers containing the values packed into words.
Definition: mzd.h:74
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
static wi_t mzd_rows_in_block(mzd_t const *M, int n)
Total number of rows in this block.
Definition: mzd.h:253
static wi_t const mzd_paddingwidth
The minimum width where padding occurs.
Definition: mzd.h:144
wi_t offset_vector
Definition: mzd.h:108
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:1758
wi_t row_offset
Definition: mzd.h:110
static int mzd_row_to_block(mzd_t const *M, rci_t row)
Convert row to blocks index.
Definition: mzd.h:236
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:1196
mzd_t * mzd_extract_u(mzd_t *U, mzd_t const *A)
Definition: mzd.c:2215
static uint8_t const mzd_flag_multiple_blocks
flag for multiply blocks
Definition: mzd.h:174
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:415
static int mzd_is_windowed(mzd_t const *M)
Test if a matrix is windowed.
Definition: mzd.h:183
static int mzd_owns_blocks(mzd_t const *M)
Test if this mzd_t should free blocks.
Definition: mzd.h:194
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 word * mzd_first_row(mzd_t const *M)
Get a pointer the first word.
Definition: mzd.h:206
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:2154
uint8_t blockrows_log
Definition: mzd.h:131
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:1101
static word const m4ri_one
The number one as a word.
Definition: misc.h:147
int mzd_is_zero(mzd_t const *A)
Zero test for matrix.
Definition: mzd.c:1985
double mzd_density(mzd_t const *A, wi_t res)
Return the number of nonzero entries divided by nrows * ncols.
Definition: mzd.c:2191
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:1891
word * end
Definition: mzd.h:77
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:1566
wi_t rowstride
Definition: mzd.h:99
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:2020
static uint8_t const mzd_flag_nonzero_excess
flag when ncols%64 == 0
Definition: mzd.h:150
static word * mzd_row(mzd_t const *M, rci_t row)
Get pointer to first word of row.
Definition: mzd.h:301
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:660
rci_t ncols
Definition: mzd.h:89
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:566
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:1611
#define __M4RI_UNLIKELY(cond)
Macro to help with branch prediction.
Definition: misc.h:449
wi_t width
Definition: mzd.h:90
uint64_t word
A word is the typical packed data structure to represent packed bits.
Definition: misc.h:87
int wi_t
Type of word indexes.
Definition: misc.h:80
#define __M4RI_GET_BIT(w, spot)
Get the bit spot (counting from the left) in the word w.
Definition: misc.h:226
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
static uint8_t const mzd_flag_windowed_ownsblocks
flag for windowed matrix wich owns its memory
Definition: mzd.h:168
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:597
#define __M4RI_CONVERT_TO_INT(w)
Explicit conversion macro.
Definition: misc.h:99
word * begin
Definition: mzd.h:76
rci_t mzd_first_zero_row(mzd_t const *A)
Return the first row with all zero entries.
Definition: mzd.c:2195