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 
74 typedef struct {
75  size_t size;
77  word* end;
78 } mzd_block_t;
79 
86 typedef struct mzd_t {
87 
100 
109 
124  uint8_t flags;
125 
131  uint8_t blockrows_log;
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 
138  word **rows;
139 } mzd_t;
140 
144 static wi_t const mzd_paddingwidth = 1;
145 
150 static uint8_t const mzd_flag_nonzero_excess = 0x2;
151 
156 static uint8_t const mzd_flag_windowed_zerooffset = 0x4;
157 
162 static uint8_t const mzd_flag_windowed_zeroexcess = 0x8;
163 
168 static uint8_t const mzd_flag_windowed_ownsblocks = 0x10;
169 
174 static uint8_t const mzd_flag_multiple_blocks = 0x20;
175 
183 static inline int mzd_is_windowed(mzd_t const *M) {
184  return M->flags & (mzd_flag_windowed_zerooffset);
185 }
186 
194 static inline int mzd_owns_blocks(mzd_t const *M) {
195  return M->blocks && (!mzd_is_windowed(M) || ((M->flags & mzd_flag_windowed_ownsblocks)));
196 }
197 
206 static inline word* mzd_first_row(mzd_t const *M) {
207  word* result = M->blocks[0].begin + M->offset_vector;
208  assert(M->nrows == 0 || result == M->rows[0]);
209  return result;
210 }
211 
222 static inline word* mzd_first_row_next_block(mzd_t const* M, int n) {
223  assert(n > 0);
224  return M->blocks[n].begin + M->offset_vector - M->row_offset * M->rowstride;
225 }
226 
236 static inline int mzd_row_to_block(mzd_t const* M, rci_t row) {
237  return (M->row_offset + row) >> M->blockrows_log;
238 }
239 
253 static inline wi_t mzd_rows_in_block(mzd_t const* M, int n) {
254  if (__M4RI_UNLIKELY(M->flags & mzd_flag_multiple_blocks)) {
255  if (__M4RI_UNLIKELY(n == 0)) {
256  return (1 << M->blockrows_log) - M->row_offset;
257  } else {
258  int const last_block = mzd_row_to_block(M, M->nrows - 1);
259  if (n < last_block)
260  return (1 << M->blockrows_log);
261  return M->nrows + M->row_offset - (n << M->blockrows_log);
262  }
263  }
264  return n ? 0 : M->nrows;
265 }
266 
276 static inline wi_t mzd_remaining_rows_in_block(mzd_t const* M, rci_t r) {
277  const int n = mzd_row_to_block(M, r);
278  r = (r - (n << M->blockrows_log));
279  if (__M4RI_UNLIKELY(M->flags & mzd_flag_multiple_blocks)) {
280  if (__M4RI_UNLIKELY(n == 0)) {
281  return (1 << M->blockrows_log) - M->row_offset - r;
282  } else {
283  int const last_block = mzd_row_to_block(M, M->nrows - 1);
284  if (n < last_block)
285  return (1 << M->blockrows_log) - r;
286  return M->nrows + M->row_offset - (n << M->blockrows_log) - r;
287  }
288  }
289  return n ? 0 : M->nrows - r;
290 }
291 
301 static inline word* mzd_row(mzd_t const* M, rci_t row) {
302  wi_t big_vector = M->offset_vector + row * M->rowstride;
303  word* result = M->blocks[0].begin + big_vector;
304  if (__M4RI_UNLIKELY(M->flags & mzd_flag_multiple_blocks)) {
305  int const n = (M->row_offset + row) >> M->blockrows_log;
306  result = M->blocks[n].begin + big_vector - n * (M->blocks[0].size / sizeof(word));
307  }
308  assert(result == M->rows[row]);
309  return result;
310 }
311 
322 mzd_t *mzd_init(rci_t const r, rci_t const c);
323 
330 void mzd_free(mzd_t *A);
331 
332 
354 mzd_t *mzd_init_window(mzd_t *M, rci_t const lowr, rci_t const lowc, rci_t const highr, rci_t const highc);
355 
362 static inline 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)
363 {
364  return mzd_init_window((mzd_t*)M, lowr, lowc, highr, highc);
365 }
366 
373 #define mzd_free_window mzd_free
374 
384 static inline void _mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb, wi_t const startblock) {
385  if ((rowa == rowb) || (startblock >= M->width))
386  return;
387 
388  wi_t width = M->width - startblock - 1;
389  word *a = M->rows[rowa] + startblock;
390  word *b = M->rows[rowb] + startblock;
391  word tmp;
392  word const mask_end = M->high_bitmask;
393 
394  for(wi_t i = 0; i < width; ++i) {
395  tmp = a[i];
396  a[i] = b[i];
397  b[i] = tmp;
398  }
399  tmp = (a[width] ^ b[width]) & mask_end;
400  a[width] ^= tmp;
401  b[width] ^= tmp;
402 
403  __M4RI_DD_ROW(M, rowa);
404  __M4RI_DD_ROW(M, rowb);
405 }
406 
415 static inline void mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb) {
416  _mzd_row_swap(M, rowa, rowb, 0);
417 }
418 
431 void mzd_copy_row(mzd_t *B, rci_t i, mzd_t const *A, rci_t j);
432 
441 void mzd_col_swap(mzd_t *M, rci_t const cola, rci_t const colb);
442 
453 static inline 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) {
454  if (cola == colb)
455  return;
456 
457  rci_t const _cola = cola;
458  rci_t const _colb = colb;
459 
460  wi_t const a_word = _cola / m4ri_radix;
461  wi_t const b_word = _colb / m4ri_radix;
462 
463  int const a_bit = _cola % m4ri_radix;
464  int const b_bit = _colb % m4ri_radix;
465 
466  word* RESTRICT ptr = mzd_row(M, start_row);
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;
470  int block = mzd_row_to_block(M, start_row);
471  int offset = max_bit - min_bit;
472  word mask = m4ri_one << min_bit;
473  int count = MIN(mzd_remaining_rows_in_block(M, start_row), count_remaining);
474 
475  // Apparently we're calling with start_row == stop_row sometimes (seems a bug to me).
476  if (count <= 0)
477  return;
478 
479  if (a_word == b_word) {
480  while(1) {
481  count_remaining -= count;
482  ptr += a_word;
483  int fast_count = count / 4;
484  int rest_count = count - 4 * fast_count;
485  word xor_v[4];
486  wi_t const rowstride = M->rowstride;
487  while (fast_count--) {
488  xor_v[0] = ptr[0];
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;
496  xor_v[0] &= mask;
497  xor_v[1] &= mask;
498  xor_v[2] &= mask;
499  xor_v[3] &= mask;
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;
504  ptr[0] ^= xor_v[0];
505  ptr[rowstride] ^= xor_v[1];
506  ptr[2 * rowstride] ^= xor_v[2];
507  ptr[3 * rowstride] ^= xor_v[3];
508  ptr += 4 * rowstride;
509  }
510  while (rest_count--) {
511  word xor_v = *ptr;
512  xor_v ^= xor_v >> offset;
513  xor_v &= mask;
514  *ptr ^= xor_v | (xor_v << offset);
515  ptr += rowstride;
516  }
517  block++;
518  if ((count = MIN(mzd_rows_in_block(M, block), count_remaining)) <= 0)
519  break;
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  ptr = mzd_first_row_next_block(M, block);
545  if (min_bit == a_bit)
546  min_ptr = ptr + a_word;
547  else
548  min_ptr = ptr + b_word;
549  }
550  }
551 
552  __M4RI_DD_MZD(M);
553 }
554 
566 static inline BIT mzd_read_bit(mzd_t const *M, rci_t const row, rci_t const col ) {
567  return __M4RI_GET_BIT(M->rows[row][col/m4ri_radix], col%m4ri_radix);
568 }
569 
582 static inline void mzd_write_bit(mzd_t *M, rci_t const row, rci_t const col, BIT const value) {
583  __M4RI_WRITE_BIT(M->rows[row][col/m4ri_radix], col%m4ri_radix, value);
584 }
585 
586 
597 static inline void mzd_xor_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values) {
598  int const spot = y % m4ri_radix;
599  wi_t const block = y / m4ri_radix;
600  M->rows[x][block] ^= values << spot;
601  int const space = m4ri_radix - spot;
602  if (n > space)
603  M->rows[x][block + 1] ^= values >> space;
604 }
605 
616 static inline void mzd_and_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values) {
617  /* This is the best way, since this will drop out once we inverse the bits in values: */
618  values >>= (m4ri_radix - n); /* Move the bits to the lowest columns */
619 
620  int const spot = y % m4ri_radix;
621  wi_t const block = y / m4ri_radix;
622  M->rows[x][block] &= values << spot;
623  int const space = m4ri_radix - spot;
624  if (n > space)
625  M->rows[x][block + 1] &= values >> space;
626 }
627 
637 static inline void mzd_clear_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
638  assert(n>0 && n <= m4ri_radix);
639  word values = m4ri_ffff >> (m4ri_radix - n);
640  int const spot = y % m4ri_radix;
641  wi_t const block = y / m4ri_radix;
642  M->rows[x][block] &= ~(values << spot);
643  int const space = m4ri_radix - spot;
644  if (n > space)
645  M->rows[x][block + 1] &= ~(values >> space);
646 }
647 
660 static inline void mzd_row_add_offset(mzd_t *M, rci_t dstrow, rci_t srcrow, rci_t coloffset) {
661  assert(dstrow < M->nrows && srcrow < M->nrows && coloffset < M->ncols);
662  wi_t const startblock= coloffset/m4ri_radix;
663  wi_t wide = M->width - startblock;
664  word *src = M->rows[srcrow] + startblock;
665  word *dst = M->rows[dstrow] + startblock;
666  word const mask_begin = __M4RI_RIGHT_BITMASK(m4ri_radix - coloffset % m4ri_radix);
667  word const mask_end = M->high_bitmask;
668 
669  *dst++ ^= *src++ & mask_begin;
670  --wide;
671 
672 #if __M4RI_HAVE_SSE2
673  int not_aligned = __M4RI_ALIGNMENT(src,16) != 0; /* 0: Aligned, 1: Not aligned */
674  if (wide > not_aligned + 1) /* Speed up for small matrices */
675  {
676  if (not_aligned) {
677  *dst++ ^= *src++;
678  --wide;
679  }
680  /* Now wide > 1 */
681  __m128i* __src = (__m128i*)src;
682  __m128i* __dst = (__m128i*)dst;
683  __m128i* const eof = (__m128i*)((unsigned long)(src + wide) & ~0xFUL);
684  do
685  {
686  __m128i xmm1 = _mm_xor_si128(*__dst, *__src);
687  *__dst++ = xmm1;
688  }
689  while(++__src < eof);
690  src = (word*)__src;
691  dst = (word*)__dst;
692  wide = ((sizeof(word)*wide)%16)/sizeof(word);
693  }
694 #endif
695  wi_t i = -1;
696  while(++i < wide)
697  dst[i] ^= src[i];
698  /*
699  * Revert possibly non-zero excess bits.
700  * Note that i == wide here, and wide can be 0.
701  * But really, src[wide - 1] is M->rows[srcrow][M->width - 1] ;)
702  * We use i - 1 here to let the compiler know these are the same addresses
703  * that we last accessed, in the previous loop.
704  */
705  dst[i - 1] ^= src[i - 1] & ~mask_end;
706 
707  __M4RI_DD_ROW(M, dstrow);
708 }
709 
721 void mzd_row_add(mzd_t *M, rci_t const sourcerow, rci_t const destrow);
722 
737 mzd_t *mzd_transpose(mzd_t *DST, mzd_t const *A);
738 
753 mzd_t *mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);
754 
769 mzd_t *mzd_addmul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);
770 
782 mzd_t *_mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B, int const clear);
783 
793 mzd_t *_mzd_mul_va(mzd_t *C, mzd_t const *v, mzd_t const *A, int const clear);
794 
803 void mzd_randomize(mzd_t *M);
804 
819 void mzd_set_ui(mzd_t *M, unsigned int const value);
820 
834 rci_t mzd_gauss_delayed(mzd_t *M, rci_t const startcol, int const full);
835 
849 rci_t mzd_echelonize_naive(mzd_t *M, int const full);
850 
858 int mzd_equal(mzd_t const *A, mzd_t const *B);
859 
871 int mzd_cmp(mzd_t const *A, mzd_t const *B);
872 
880 mzd_t *mzd_copy(mzd_t *DST, mzd_t const *A);
881 
900 mzd_t *mzd_concat(mzd_t *C, mzd_t const *A, mzd_t const *B);
901 
919 mzd_t *mzd_stack(mzd_t *C, mzd_t const *A, mzd_t const *B);
920 
933 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);
934 
946 mzd_t *mzd_invert_naive(mzd_t *INV, mzd_t const *A, mzd_t const *I);
947 
959 mzd_t *mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);
960 
969 mzd_t *_mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);
970 
979 #define mzd_sub mzd_add
980 
989 #define _mzd_sub _mzd_add
990 
991 
992 
1002 static inline word mzd_read_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
1003  int const spot = y % m4ri_radix;
1004  wi_t const block = y / m4ri_radix;
1005  int const spill = spot + n - m4ri_radix;
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);
1008 }
1009 
1010 
1027 static inline void mzd_combine_even_in_place(mzd_t *A, rci_t const a_row, wi_t const a_startblock,
1028  mzd_t const *B, rci_t const b_row, wi_t const b_startblock) {
1029 
1030  wi_t wide = A->width - a_startblock - 1;
1031 
1032  word *a = A->rows[a_row] + a_startblock;
1033  word *b = B->rows[b_row] + b_startblock;
1034 
1035 #if __M4RI_HAVE_SSE2
1036  if(wide > 2) {
1038  if (__M4RI_ALIGNMENT(a,16)) {
1039  *a++ ^= *b++;
1040  wide--;
1041  }
1042 
1043  if (__M4RI_ALIGNMENT(a, 16) == 0 && __M4RI_ALIGNMENT(b, 16) == 0) {
1044  __m128i *a128 = (__m128i*)a;
1045  __m128i *b128 = (__m128i*)b;
1046  const __m128i *eof = (__m128i*)((unsigned long)(a + wide) & ~0xFUL);
1047 
1048  do {
1049  *a128 = _mm_xor_si128(*a128, *b128);
1050  ++b128;
1051  ++a128;
1052  } while(a128 < eof);
1053 
1054  a = (word*)a128;
1055  b = (word*)b128;
1056  wide = ((sizeof(word) * wide) % 16) / sizeof(word);
1057  }
1058  }
1059 #endif // __M4RI_HAVE_SSE2
1060 
1061  if (wide > 0) {
1062  wi_t n = (wide + 7) / 8;
1063  switch (wide % 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++);
1072  } while (--n > 0);
1073  }
1074  }
1075 
1076  *a ^= *b & A->high_bitmask;
1077 
1078  __M4RI_DD_MZD(A);
1079 }
1080 
1081 
1101 static inline void mzd_combine_even(mzd_t *C, rci_t const c_row, wi_t const c_startblock,
1102  mzd_t const *A, rci_t const a_row, wi_t const a_startblock,
1103  mzd_t const *B, rci_t const b_row, wi_t const b_startblock) {
1104 
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;
1109 
1110 #if __M4RI_HAVE_SSE2
1111  if(wide > 2) {
1113  if (__M4RI_ALIGNMENT(a,16)) {
1114  *c++ = *b++ ^ *a++;
1115  wide--;
1116  }
1117 
1118  if ( (__M4RI_ALIGNMENT(b, 16) | __M4RI_ALIGNMENT(c, 16)) == 0) {
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);
1123 
1124  do {
1125  *c128 = _mm_xor_si128(*a128, *b128);
1126  ++c128;
1127  ++b128;
1128  ++a128;
1129  } while(a128 < eof);
1130 
1131  a = (word*)a128;
1132  b = (word*)b128;
1133  c = (word*)c128;
1134  wide = ((sizeof(word) * wide) % 16) / sizeof(word);
1135  }
1136  }
1137 #endif // __M4RI_HAVE_SSE2
1138 
1139  if (wide > 0) {
1140  wi_t n = (wide + 7) / 8;
1141  switch (wide % 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++);
1150  } while (--n > 0);
1151  }
1152  }
1153  *c ^= ((*a ^ *b ^ *c) & C->high_bitmask);
1154 
1155  __M4RI_DD_MZD(C);
1156 }
1157 
1158 
1177 static inline void mzd_combine(mzd_t *C, rci_t const c_row, wi_t const c_startblock,
1178  mzd_t const *A, rci_t const a_row, wi_t const a_startblock,
1179  mzd_t const *B, rci_t const b_row, wi_t const b_startblock) {
1180 
1181  if( (C == A) & (a_row == c_row) & (a_startblock == c_startblock) )
1182  mzd_combine_even_in_place(C, c_row, c_startblock, B, b_row, b_startblock);
1183  else
1184  mzd_combine_even(C, c_row, c_startblock, A, a_row, a_startblock, B, b_row, b_startblock);
1185  return;
1186 }
1187 
1196 static inline int mzd_read_bits_int(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
1197  return __M4RI_CONVERT_TO_INT(mzd_read_bits(M, x, y, n));
1198 }
1199 
1200 
1207 int mzd_is_zero(mzd_t const *A);
1208 
1217 void mzd_row_clear_offset(mzd_t *M, rci_t const row, rci_t const coloffset);
1218 
1235 int mzd_find_pivot(mzd_t const *M, rci_t start_row, rci_t start_col, rci_t *r, rci_t *c);
1236 
1237 
1250 double mzd_density(mzd_t const *A, wi_t res);
1251 
1266 double _mzd_density(mzd_t const *A, wi_t res, rci_t r, rci_t c);
1267 
1268 
1277 rci_t mzd_first_zero_row(mzd_t const *A);
1278 
1285 static inline word mzd_hash(mzd_t const *A) {
1286  word hash = 0;
1287  for (rci_t r = 0; r < A->nrows; ++r)
1288  hash ^= rotate_word(calculate_hash(A->rows[r], A->width), r % m4ri_radix);
1289  return hash;
1290 }
1291 
1301 mzd_t *mzd_extract_u(mzd_t *U, mzd_t const *A);
1302 
1312 mzd_t *mzd_extract_l(mzd_t *L, mzd_t const *A);
1313 
1314 #endif // M4RI_MZD
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