1314 lines
34 KiB
C
1314 lines
34 KiB
C
/**
|
|
* \file mzd.h
|
|
* \brief Dense matrices over GF(2) represented as a bit field.
|
|
*
|
|
* \author Gregory Bard <bard@fordham.edu>
|
|
* \author Martin Albrecht <martinralbrecht+m4ri@googlemail.com>
|
|
* \author Carlo Wood <carlo@alinoe.com>
|
|
*/
|
|
|
|
#ifndef M4RI_MZD
|
|
#define M4RI_MZD
|
|
|
|
/*******************************************************************
|
|
*
|
|
* M4RI: Linear Algebra over GF(2)
|
|
*
|
|
* Copyright (C) 2007, 2008 Gregory Bard <bard@fordham.edu>
|
|
* Copyright (C) 2008-2013 Martin Albrecht <M.R.Albrecht@rhul.ac.uk>
|
|
* Copyright (C) 2011 Carlo Wood <carlo@alinoe.com>
|
|
*
|
|
* Distributed under the terms of the GNU General Public License (GPL)
|
|
* version 2 or higher.
|
|
*
|
|
* This code is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
* General Public License for more details.
|
|
*
|
|
* The full text of the GPL is available at:
|
|
*
|
|
* http://www.gnu.org/licenses/
|
|
*
|
|
********************************************************************/
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include <m4ri/m4ri_config.h>
|
|
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
|
|
#if __M4RI_HAVE_SSE2
|
|
#include <emmintrin.h>
|
|
#endif
|
|
|
|
#include <m4ri/misc.h>
|
|
#include <m4ri/debug_dump.h>
|
|
|
|
/**
|
|
* Maximum number of words allocated for one mzd_t block.
|
|
*
|
|
* \note This value must fit in an int, even though it's type is size_t.
|
|
*/
|
|
|
|
#define __M4RI_MAX_MZD_BLOCKSIZE (((size_t)1) << 27)
|
|
|
|
/**
|
|
* \brief Matrix multiplication block-ing dimension.
|
|
*
|
|
* Defines the number of rows of the matrix A that are
|
|
* processed as one block during the execution of a multiplication
|
|
* algorithm.
|
|
*/
|
|
|
|
#define __M4RI_MUL_BLOCKSIZE MIN(((int)sqrt((double)(4 * __M4RI_CPU_L3_CACHE))) / 2, 2048)
|
|
|
|
/**
|
|
* \brief Data containers containing the values packed into words
|
|
*/
|
|
|
|
typedef struct {
|
|
size_t size; /*!< number of words */
|
|
word* begin; /*!< first word */
|
|
word* end; /*!< last word */
|
|
} mzd_block_t;
|
|
|
|
/**
|
|
* \brief Dense matrices over GF(2).
|
|
*
|
|
* The most fundamental data type in this library.
|
|
*/
|
|
|
|
typedef struct mzd_t {
|
|
|
|
rci_t nrows; /*!< Number of rows. */
|
|
rci_t ncols; /*!< Number of columns. */
|
|
wi_t width; /*!< Number of words with valid bits: width = ceil(ncols / m4ri_radix) */
|
|
|
|
/**
|
|
* Offset in words between rows.
|
|
*
|
|
* rowstride = (width < mzd_paddingwidth || (width & 1) == 0) ? width : width + 1;
|
|
* where width is the width of the underlying non-windowed matrix.
|
|
*/
|
|
|
|
wi_t rowstride;
|
|
|
|
/**
|
|
* Offset in words from start of block to first word.
|
|
*
|
|
* rows[0] = blocks[0].begin + offset_vector;
|
|
* This, together with rowstride, makes the rows array obsolete.
|
|
*/
|
|
|
|
wi_t offset_vector;
|
|
|
|
wi_t row_offset; /*!< Number of rows to the first row counting from the start of the first block. */
|
|
|
|
/**
|
|
* Booleans to speed up things.
|
|
*
|
|
* The bits have the following meaning:
|
|
*
|
|
* 1: Has non-zero excess.
|
|
* 2: Is windowed, but has zero offset.
|
|
* 3: Is windowed, but has zero excess.
|
|
* 4: Is windowed, but owns the blocks allocations.
|
|
* 5: Spans more than 1 block.
|
|
*/
|
|
|
|
uint8_t flags;
|
|
|
|
/**
|
|
* blockrows_log = log2(blockrows);
|
|
* where blockrows is the number of rows in one block, which is a power of 2.
|
|
*/
|
|
|
|
uint8_t blockrows_log;
|
|
|
|
word high_bitmask; /*!< Mask for valid bits in the word with the highest index (width - 1). */
|
|
mzd_block_t *blocks; /*!< Pointers to the actual blocks of memory containing the values packed into words. */
|
|
word **rows; /*!< Address of first word in each row, so the first word of row i is is m->rows[i] */
|
|
uint64_t dummy; /*!< ensures sizeof(mzd_t) == 64 */
|
|
|
|
} mzd_t;
|
|
|
|
/**
|
|
* \brief The minimum width where padding occurs.
|
|
*/
|
|
static wi_t const mzd_paddingwidth = 1;
|
|
|
|
/**
|
|
* \brief flag when ncols%64 == 0
|
|
*/
|
|
|
|
static uint8_t const mzd_flag_nonzero_excess = 0x2;
|
|
|
|
/**
|
|
* \brief flag for windowed matrix
|
|
*/
|
|
|
|
static uint8_t const mzd_flag_windowed_zerooffset = 0x4;
|
|
|
|
/**
|
|
* \brief flag for windowed matrix where ncols%64 == 0
|
|
*/
|
|
|
|
static uint8_t const mzd_flag_windowed_zeroexcess = 0x8;
|
|
|
|
/**
|
|
* \brief flag for windowed matrix wich owns its memory
|
|
*/
|
|
|
|
static uint8_t const mzd_flag_windowed_ownsblocks = 0x10;
|
|
|
|
/**
|
|
* \brief flag for multiply blocks
|
|
*/
|
|
|
|
static uint8_t const mzd_flag_multiple_blocks = 0x20;
|
|
|
|
/**
|
|
* \brief Test if a matrix is windowed.
|
|
*
|
|
* \param M Matrix
|
|
*
|
|
* \return a non-zero value if the matrix is windowed, otherwise return zero.
|
|
*/
|
|
static inline int mzd_is_windowed(mzd_t const *M) {
|
|
return M->flags & (mzd_flag_windowed_zerooffset);
|
|
}
|
|
|
|
/**
|
|
* \brief Test if this mzd_t should free blocks.
|
|
*
|
|
* \param M Matrix
|
|
*
|
|
* \return TRUE iff blocks is non-zero and should be freed upon a call to mzd_free.
|
|
*/
|
|
static inline int mzd_owns_blocks(mzd_t const *M) {
|
|
return M->blocks && (!mzd_is_windowed(M) || ((M->flags & mzd_flag_windowed_ownsblocks)));
|
|
}
|
|
|
|
/**
|
|
* \brief Get a pointer the first word.
|
|
*
|
|
* \param M Matrix
|
|
*
|
|
* \return a pointer to the first word of the first row.
|
|
*/
|
|
|
|
static inline word* mzd_first_row(mzd_t const *M) {
|
|
word* result = M->blocks[0].begin + M->offset_vector;
|
|
assert(M->nrows == 0 || result == M->rows[0]);
|
|
return result;
|
|
}
|
|
|
|
/**
|
|
* \brief Get a pointer to the first word in block n.
|
|
*
|
|
* Use mzd_first_row for block number 0.
|
|
*
|
|
* \param M Matrix
|
|
* \param n The block number. Must be larger than 0.
|
|
*
|
|
* \return a pointer to the first word of the first row in block n.
|
|
*/
|
|
static inline word* mzd_first_row_next_block(mzd_t const* M, int n) {
|
|
assert(n > 0);
|
|
return M->blocks[n].begin + M->offset_vector - M->row_offset * M->rowstride;
|
|
}
|
|
|
|
/**
|
|
* \brief Convert row to blocks index.
|
|
*
|
|
* \param M Matrix.
|
|
* \param row The row to convert.
|
|
*
|
|
* \return the block number that contains this row.
|
|
*/
|
|
|
|
static inline int mzd_row_to_block(mzd_t const* M, rci_t row) {
|
|
return (M->row_offset + row) >> M->blockrows_log;
|
|
}
|
|
|
|
/**
|
|
* \brief Total number of rows in this block.
|
|
*
|
|
* Should be called with a constant n=0, or with
|
|
* n > 0 when n is a variable, for optimization
|
|
* reasons.
|
|
*
|
|
* \param M Matrix
|
|
* \param n The block number.
|
|
*
|
|
* \return the total number of rows in this block.
|
|
*/
|
|
|
|
static inline wi_t mzd_rows_in_block(mzd_t const* M, int n) {
|
|
if (__M4RI_UNLIKELY(M->flags & mzd_flag_multiple_blocks)) {
|
|
if (__M4RI_UNLIKELY(n == 0)) {
|
|
return (1 << M->blockrows_log) - M->row_offset;
|
|
} else {
|
|
int const last_block = mzd_row_to_block(M, M->nrows - 1);
|
|
if (n < last_block)
|
|
return (1 << M->blockrows_log);
|
|
return M->nrows + M->row_offset - (n << M->blockrows_log);
|
|
}
|
|
}
|
|
return n ? 0 : M->nrows;
|
|
}
|
|
|
|
/**
|
|
* \brief Number of rows in this block including r
|
|
*
|
|
* \param M Matrix
|
|
* \param r row
|
|
*
|
|
* \return the number of rows with index >= r in this block
|
|
*/
|
|
|
|
static inline wi_t mzd_remaining_rows_in_block(mzd_t const* M, rci_t r) {
|
|
const int n = mzd_row_to_block(M, r);
|
|
r = (r - (n << M->blockrows_log));
|
|
if (__M4RI_UNLIKELY(M->flags & mzd_flag_multiple_blocks)) {
|
|
if (__M4RI_UNLIKELY(n == 0)) {
|
|
return (1 << M->blockrows_log) - M->row_offset - r;
|
|
} else {
|
|
int const last_block = mzd_row_to_block(M, M->nrows - 1);
|
|
if (n < last_block)
|
|
return (1 << M->blockrows_log) - r;
|
|
return M->nrows + M->row_offset - (n << M->blockrows_log) - r;
|
|
}
|
|
}
|
|
return n ? 0 : M->nrows - r;
|
|
}
|
|
|
|
/**
|
|
* \brief Get pointer to first word of row.
|
|
*
|
|
* \param M Matrix
|
|
* \param row The row index.
|
|
*
|
|
* \return pointer to first word of the row.
|
|
*/
|
|
|
|
static inline word* mzd_row(mzd_t const* M, rci_t row) {
|
|
wi_t big_vector = M->offset_vector + row * M->rowstride;
|
|
word* result = M->blocks[0].begin + big_vector;
|
|
if (__M4RI_UNLIKELY(M->flags & mzd_flag_multiple_blocks)) {
|
|
int const n = (M->row_offset + row) >> M->blockrows_log;
|
|
result = M->blocks[n].begin + big_vector - n * (M->blocks[0].size / sizeof(word));
|
|
}
|
|
assert(result == M->rows[row]);
|
|
return result;
|
|
}
|
|
|
|
/**
|
|
* \brief Create a new matrix of dimension r x c.
|
|
*
|
|
* Use mzd_free to kill it.
|
|
*
|
|
* \param r Number of rows
|
|
* \param c Number of columns
|
|
*
|
|
*/
|
|
|
|
mzd_t *mzd_init(rci_t const r, rci_t const c);
|
|
|
|
/**
|
|
* \brief Free a matrix created with mzd_init.
|
|
*
|
|
* \param A Matrix
|
|
*/
|
|
|
|
void mzd_free(mzd_t *A);
|
|
|
|
|
|
/**
|
|
* \brief Create a window/view into the matrix M.
|
|
*
|
|
* A matrix window for M is a meta structure on the matrix M. It is
|
|
* setup to point into the matrix so M \em must \em not be freed while the
|
|
* matrix window is used.
|
|
*
|
|
* This function puts the restriction on the provided parameters that
|
|
* all parameters must be within range for M which is not enforced
|
|
* currently .
|
|
*
|
|
* Use mzd_free_window to free the window.
|
|
*
|
|
* \param M Matrix
|
|
* \param lowr Starting row (inclusive)
|
|
* \param lowc Starting column (inclusive, must be multiple of m4ri_radix)
|
|
* \param highr End row (exclusive)
|
|
* \param highc End column (exclusive)
|
|
*
|
|
*/
|
|
|
|
mzd_t *mzd_init_window(mzd_t *M, rci_t const lowr, rci_t const lowc, rci_t const highr, rci_t const highc);
|
|
|
|
/**
|
|
* \brief Create a const window/view into a const matrix M.
|
|
*
|
|
* See mzd_init_window, but for constant M.
|
|
*/
|
|
|
|
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)
|
|
{
|
|
return mzd_init_window((mzd_t*)M, lowr, lowc, highr, highc);
|
|
}
|
|
|
|
/**
|
|
* \brief Free a matrix window created with mzd_init_window.
|
|
*
|
|
* \param A Matrix
|
|
*/
|
|
|
|
#define mzd_free_window mzd_free
|
|
|
|
/**
|
|
* \brief Swap the two rows rowa and rowb starting at startblock.
|
|
*
|
|
* \param M Matrix with a zero offset.
|
|
* \param rowa Row index.
|
|
* \param rowb Row index.
|
|
* \param startblock Start swapping only in this block.
|
|
*/
|
|
|
|
static inline void _mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb, wi_t const startblock) {
|
|
if ((rowa == rowb) || (startblock >= M->width))
|
|
return;
|
|
|
|
wi_t width = M->width - startblock - 1;
|
|
word *a = M->rows[rowa] + startblock;
|
|
word *b = M->rows[rowb] + startblock;
|
|
word tmp;
|
|
word const mask_end = M->high_bitmask;
|
|
|
|
for(wi_t i = 0; i < width; ++i) {
|
|
tmp = a[i];
|
|
a[i] = b[i];
|
|
b[i] = tmp;
|
|
}
|
|
tmp = (a[width] ^ b[width]) & mask_end;
|
|
a[width] ^= tmp;
|
|
b[width] ^= tmp;
|
|
|
|
__M4RI_DD_ROW(M, rowa);
|
|
__M4RI_DD_ROW(M, rowb);
|
|
}
|
|
|
|
/**
|
|
* \brief Swap the two rows rowa and rowb.
|
|
*
|
|
* \param M Matrix
|
|
* \param rowa Row index.
|
|
* \param rowb Row index.
|
|
*/
|
|
|
|
static inline void mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb) {
|
|
_mzd_row_swap(M, rowa, rowb, 0);
|
|
}
|
|
|
|
/**
|
|
* \brief copy row j from A to row i from B.
|
|
*
|
|
* The offsets of A and B must match and the number of columns of A
|
|
* must be less than or equal to the number of columns of B.
|
|
*
|
|
* \param B Target matrix.
|
|
* \param i Target row index.
|
|
* \param A Source matrix.
|
|
* \param j Source row index.
|
|
*/
|
|
|
|
void mzd_copy_row(mzd_t *B, rci_t i, mzd_t const *A, rci_t j);
|
|
|
|
/**
|
|
* \brief Swap the two columns cola and colb.
|
|
*
|
|
* \param M Matrix.
|
|
* \param cola Column index.
|
|
* \param colb Column index.
|
|
*/
|
|
|
|
void mzd_col_swap(mzd_t *M, rci_t const cola, rci_t const colb);
|
|
|
|
/**
|
|
* \brief Swap the two columns cola and colb but only between start_row and stop_row.
|
|
*
|
|
* \param M Matrix.
|
|
* \param cola Column index.
|
|
* \param colb Column index.
|
|
* \param start_row Row index.
|
|
* \param stop_row Row index (exclusive).
|
|
*/
|
|
|
|
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) {
|
|
if (cola == colb)
|
|
return;
|
|
|
|
rci_t const _cola = cola;
|
|
rci_t const _colb = colb;
|
|
|
|
wi_t const a_word = _cola / m4ri_radix;
|
|
wi_t const b_word = _colb / m4ri_radix;
|
|
|
|
int const a_bit = _cola % m4ri_radix;
|
|
int const b_bit = _colb % m4ri_radix;
|
|
|
|
word* RESTRICT ptr = mzd_row(M, start_row);
|
|
int max_bit = MAX(a_bit, b_bit);
|
|
int count_remaining = stop_row - start_row;
|
|
int min_bit = a_bit + b_bit - max_bit;
|
|
int block = mzd_row_to_block(M, start_row);
|
|
int offset = max_bit - min_bit;
|
|
word mask = m4ri_one << min_bit;
|
|
int count = MIN(mzd_remaining_rows_in_block(M, start_row), count_remaining);
|
|
|
|
// Apparently we're calling with start_row == stop_row sometimes (seems a bug to me).
|
|
if (count <= 0)
|
|
return;
|
|
|
|
if (a_word == b_word) {
|
|
while(1) {
|
|
count_remaining -= count;
|
|
ptr += a_word;
|
|
int fast_count = count / 4;
|
|
int rest_count = count - 4 * fast_count;
|
|
word xor_v[4];
|
|
wi_t const rowstride = M->rowstride;
|
|
while (fast_count--) {
|
|
xor_v[0] = ptr[0];
|
|
xor_v[1] = ptr[rowstride];
|
|
xor_v[2] = ptr[2 * rowstride];
|
|
xor_v[3] = ptr[3 * rowstride];
|
|
xor_v[0] ^= xor_v[0] >> offset;
|
|
xor_v[1] ^= xor_v[1] >> offset;
|
|
xor_v[2] ^= xor_v[2] >> offset;
|
|
xor_v[3] ^= xor_v[3] >> offset;
|
|
xor_v[0] &= mask;
|
|
xor_v[1] &= mask;
|
|
xor_v[2] &= mask;
|
|
xor_v[3] &= mask;
|
|
xor_v[0] |= xor_v[0] << offset;
|
|
xor_v[1] |= xor_v[1] << offset;
|
|
xor_v[2] |= xor_v[2] << offset;
|
|
xor_v[3] |= xor_v[3] << offset;
|
|
ptr[0] ^= xor_v[0];
|
|
ptr[rowstride] ^= xor_v[1];
|
|
ptr[2 * rowstride] ^= xor_v[2];
|
|
ptr[3 * rowstride] ^= xor_v[3];
|
|
ptr += 4 * rowstride;
|
|
}
|
|
while (rest_count--) {
|
|
word xor_v = *ptr;
|
|
xor_v ^= xor_v >> offset;
|
|
xor_v &= mask;
|
|
*ptr ^= xor_v | (xor_v << offset);
|
|
ptr += rowstride;
|
|
}
|
|
block++;
|
|
if ((count = MIN(mzd_rows_in_block(M, block), count_remaining)) <= 0)
|
|
break;
|
|
ptr = mzd_first_row_next_block(M, block);
|
|
}
|
|
} else {
|
|
word* RESTRICT min_ptr;
|
|
wi_t max_offset;
|
|
if (min_bit == a_bit) {
|
|
min_ptr = ptr + a_word;
|
|
max_offset = b_word - a_word;
|
|
} else {
|
|
min_ptr = ptr + b_word;
|
|
max_offset = a_word - b_word;
|
|
}
|
|
while(1) {
|
|
count_remaining -= count;
|
|
wi_t const rowstride = M->rowstride;
|
|
while(count--) {
|
|
word xor_v = (min_ptr[0] ^ (min_ptr[max_offset] >> offset)) & mask;
|
|
min_ptr[0] ^= xor_v;
|
|
min_ptr[max_offset] ^= xor_v << offset;
|
|
min_ptr += rowstride;
|
|
}
|
|
block++;
|
|
if ((count = MIN(mzd_rows_in_block(M,+block), count_remaining)) <= 0)
|
|
break;
|
|
ptr = mzd_first_row_next_block(M, block);
|
|
if (min_bit == a_bit)
|
|
min_ptr = ptr + a_word;
|
|
else
|
|
min_ptr = ptr + b_word;
|
|
}
|
|
}
|
|
|
|
__M4RI_DD_MZD(M);
|
|
}
|
|
|
|
/**
|
|
* \brief Read the bit at position M[row,col].
|
|
*
|
|
* \param M Matrix
|
|
* \param row Row index
|
|
* \param col Column index
|
|
*
|
|
* \note No bounds checks whatsoever are performed.
|
|
*
|
|
*/
|
|
|
|
static inline BIT mzd_read_bit(mzd_t const *M, rci_t const row, rci_t const col ) {
|
|
return __M4RI_GET_BIT(M->rows[row][col/m4ri_radix], col%m4ri_radix);
|
|
}
|
|
|
|
/**
|
|
* \brief Write the bit value to position M[row,col]
|
|
*
|
|
* \param M Matrix
|
|
* \param row Row index
|
|
* \param col Column index
|
|
* \param value Either 0 or 1
|
|
*
|
|
* \note No bounds checks whatsoever are performed.
|
|
*
|
|
*/
|
|
|
|
static inline void mzd_write_bit(mzd_t *M, rci_t const row, rci_t const col, BIT const value) {
|
|
__M4RI_WRITE_BIT(M->rows[row][col/m4ri_radix], col%m4ri_radix, value);
|
|
}
|
|
|
|
|
|
/**
|
|
* \brief XOR n bits from values to M starting a position (x,y).
|
|
*
|
|
* \param M Source matrix.
|
|
* \param x Starting row.
|
|
* \param y Starting column.
|
|
* \param n Number of bits (<= m4ri_radix);
|
|
* \param values Word with values;
|
|
*/
|
|
|
|
static inline void mzd_xor_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values) {
|
|
int const spot = y % m4ri_radix;
|
|
wi_t const block = y / m4ri_radix;
|
|
M->rows[x][block] ^= values << spot;
|
|
int const space = m4ri_radix - spot;
|
|
if (n > space)
|
|
M->rows[x][block + 1] ^= values >> space;
|
|
}
|
|
|
|
/**
|
|
* \brief AND n bits from values to M starting a position (x,y).
|
|
*
|
|
* \param M Source matrix.
|
|
* \param x Starting row.
|
|
* \param y Starting column.
|
|
* \param n Number of bits (<= m4ri_radix);
|
|
* \param values Word with values;
|
|
*/
|
|
|
|
static inline void mzd_and_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n, word values) {
|
|
/* This is the best way, since this will drop out once we inverse the bits in values: */
|
|
values >>= (m4ri_radix - n); /* Move the bits to the lowest columns */
|
|
|
|
int const spot = y % m4ri_radix;
|
|
wi_t const block = y / m4ri_radix;
|
|
M->rows[x][block] &= values << spot;
|
|
int const space = m4ri_radix - spot;
|
|
if (n > space)
|
|
M->rows[x][block + 1] &= values >> space;
|
|
}
|
|
|
|
/**
|
|
* \brief Clear n bits in M starting a position (x,y).
|
|
*
|
|
* \param M Source matrix.
|
|
* \param x Starting row.
|
|
* \param y Starting column.
|
|
* \param n Number of bits (0 < n <= m4ri_radix);
|
|
*/
|
|
|
|
static inline void mzd_clear_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
|
|
assert(n>0 && n <= m4ri_radix);
|
|
word values = m4ri_ffff >> (m4ri_radix - n);
|
|
int const spot = y % m4ri_radix;
|
|
wi_t const block = y / m4ri_radix;
|
|
M->rows[x][block] &= ~(values << spot);
|
|
int const space = m4ri_radix - spot;
|
|
if (n > space)
|
|
M->rows[x][block + 1] &= ~(values >> space);
|
|
}
|
|
|
|
/**
|
|
* \brief Add the rows sourcerow and destrow and stores the total in the row
|
|
* destrow, but only begins at the column coloffset.
|
|
*
|
|
* \param M Matrix
|
|
* \param dstrow Index of target row
|
|
* \param srcrow Index of source row
|
|
* \param coloffset Start column (0 <= coloffset < M->ncols)
|
|
*
|
|
* \warning This function expects that there is at least one word worth of work.
|
|
*/
|
|
|
|
static inline void mzd_row_add_offset(mzd_t *M, rci_t dstrow, rci_t srcrow, rci_t coloffset) {
|
|
assert(dstrow < M->nrows && srcrow < M->nrows && coloffset < M->ncols);
|
|
wi_t const startblock= coloffset/m4ri_radix;
|
|
wi_t wide = M->width - startblock;
|
|
word *src = M->rows[srcrow] + startblock;
|
|
word *dst = M->rows[dstrow] + startblock;
|
|
word const mask_begin = __M4RI_RIGHT_BITMASK(m4ri_radix - coloffset % m4ri_radix);
|
|
word const mask_end = M->high_bitmask;
|
|
|
|
*dst++ ^= *src++ & mask_begin;
|
|
--wide;
|
|
|
|
#if __M4RI_HAVE_SSE2
|
|
int not_aligned = __M4RI_ALIGNMENT(src,16) != 0; /* 0: Aligned, 1: Not aligned */
|
|
if (wide > not_aligned + 1) /* Speed up for small matrices */
|
|
{
|
|
if (not_aligned) {
|
|
*dst++ ^= *src++;
|
|
--wide;
|
|
}
|
|
/* Now wide > 1 */
|
|
__m128i* __src = (__m128i*)src;
|
|
__m128i* __dst = (__m128i*)dst;
|
|
__m128i* const eof = (__m128i*)((unsigned long)(src + wide) & ~0xFUL);
|
|
do
|
|
{
|
|
__m128i xmm1 = _mm_xor_si128(*__dst, *__src);
|
|
*__dst++ = xmm1;
|
|
}
|
|
while(++__src < eof);
|
|
src = (word*)__src;
|
|
dst = (word*)__dst;
|
|
wide = ((sizeof(word)*wide)%16)/sizeof(word);
|
|
}
|
|
#endif
|
|
wi_t i = -1;
|
|
while(++i < wide)
|
|
dst[i] ^= src[i];
|
|
/*
|
|
* Revert possibly non-zero excess bits.
|
|
* Note that i == wide here, and wide can be 0.
|
|
* But really, src[wide - 1] is M->rows[srcrow][M->width - 1] ;)
|
|
* We use i - 1 here to let the compiler know these are the same addresses
|
|
* that we last accessed, in the previous loop.
|
|
*/
|
|
dst[i - 1] ^= src[i - 1] & ~mask_end;
|
|
|
|
__M4RI_DD_ROW(M, dstrow);
|
|
}
|
|
|
|
/**
|
|
* \brief Add the rows sourcerow and destrow and stores the total in
|
|
* the row destrow.
|
|
*
|
|
* \param M Matrix
|
|
* \param sourcerow Index of source row
|
|
* \param destrow Index of target row
|
|
*
|
|
* \note this can be done much faster with mzd_combine.
|
|
*/
|
|
|
|
void mzd_row_add(mzd_t *M, rci_t const sourcerow, rci_t const destrow);
|
|
|
|
/**
|
|
* \brief Transpose a matrix.
|
|
*
|
|
* This function uses the fact that:
|
|
\verbatim
|
|
[ A B ]T [AT CT]
|
|
[ C D ] = [BT DT]
|
|
\endverbatim
|
|
* and thus rearranges the blocks recursively.
|
|
*
|
|
* \param DST Preallocated return matrix, may be NULL for automatic creation.
|
|
* \param A Matrix
|
|
*/
|
|
|
|
mzd_t *mzd_transpose(mzd_t *DST, mzd_t const *A);
|
|
|
|
/**
|
|
* \brief Naive cubic matrix multiplication.
|
|
*
|
|
* That is, compute C such that C == AB.
|
|
*
|
|
* \param C Preallocated product matrix, may be NULL for automatic creation.
|
|
* \param A Input matrix A.
|
|
* \param B Input matrix B.
|
|
*
|
|
* \note Normally, if you will multiply several times by b, it is
|
|
* smarter to calculate bT yourself, and keep it, and then use the
|
|
* function called _mzd_mul_naive
|
|
*
|
|
*/
|
|
mzd_t *mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Naive cubic matrix multiplication and addition
|
|
*
|
|
* That is, compute C such that C == C + AB.
|
|
*
|
|
* \param C Preallocated product matrix.
|
|
* \param A Input matrix A.
|
|
* \param B Input matrix B.
|
|
*
|
|
* \note Normally, if you will multiply several times by b, it is
|
|
* smarter to calculate bT yourself, and keep it, and then use the
|
|
* function called _mzd_mul_naive
|
|
*/
|
|
|
|
mzd_t *mzd_addmul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Naive cubic matrix multiplication with the pre-transposed B.
|
|
*
|
|
* That is, compute C such that C == AB^t.
|
|
*
|
|
* \param C Preallocated product matrix.
|
|
* \param A Input matrix A.
|
|
* \param B Pre-transposed input matrix B.
|
|
* \param clear Whether to clear C before accumulating AB
|
|
*/
|
|
|
|
mzd_t *_mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B, int const clear);
|
|
|
|
/**
|
|
* \brief Matrix multiplication optimized for v*A where v is a vector.
|
|
*
|
|
* \param C Preallocated product matrix.
|
|
* \param v Input matrix v.
|
|
* \param A Input matrix A.
|
|
* \param clear If set clear C first, otherwise add result to C.
|
|
*
|
|
*/
|
|
mzd_t *_mzd_mul_va(mzd_t *C, mzd_t const *v, mzd_t const *A, int const clear);
|
|
|
|
/**
|
|
* \brief Fill matrix M with uniformly distributed bits.
|
|
*
|
|
* \param M Matrix
|
|
*
|
|
* \todo Allow the user to provide a RNG callback.
|
|
*/
|
|
|
|
void mzd_randomize(mzd_t *M);
|
|
|
|
/**
|
|
* \brief Set the matrix M to the value equivalent to the integer
|
|
* value provided.
|
|
*
|
|
* Specifically, this function does nothing if value%2 == 0 and
|
|
* returns the identity matrix if value%2 == 1.
|
|
*
|
|
* If the matrix is not square then the largest possible square
|
|
* submatrix is set to the identity matrix.
|
|
*
|
|
* \param M Matrix
|
|
* \param value Either 0 or 1
|
|
*/
|
|
|
|
void mzd_set_ui(mzd_t *M, unsigned int const value);
|
|
|
|
/**
|
|
* \brief Gaussian elimination.
|
|
*
|
|
* This will do Gaussian elimination on the matrix m but will start
|
|
* not at column 0 necc but at column startcol. If full=FALSE, then it
|
|
* will do triangular style elimination, and if full=TRUE, it will do
|
|
* Gauss-Jordan style, or full elimination.
|
|
*
|
|
* \param M Matrix
|
|
* \param startcol First column to consider for reduction.
|
|
* \param full Gauss-Jordan style or upper triangular form only.
|
|
*/
|
|
|
|
rci_t mzd_gauss_delayed(mzd_t *M, rci_t const startcol, int const full);
|
|
|
|
/**
|
|
* \brief Gaussian elimination.
|
|
*
|
|
* This will do Gaussian elimination on the matrix m. If full=FALSE,
|
|
* then it will do triangular style elimination, and if full=TRUE,
|
|
* it will do Gauss-Jordan style, or full elimination.
|
|
*
|
|
* \param M Matrix
|
|
* \param full Gauss-Jordan style or upper triangular form only.
|
|
*
|
|
* \sa mzd_echelonize_m4ri(), mzd_echelonize_pluq()
|
|
*/
|
|
|
|
rci_t mzd_echelonize_naive(mzd_t *M, int const full);
|
|
|
|
/**
|
|
* \brief Return TRUE if A == B.
|
|
*
|
|
* \param A Matrix
|
|
* \param B Matrix
|
|
*/
|
|
|
|
int mzd_equal(mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Return -1,0,1 if if A < B, A == B or A > B respectively.
|
|
*
|
|
* \param A Matrix.
|
|
* \param B Matrix.
|
|
*
|
|
* \note This comparison is not well defined mathematically and
|
|
* relatively arbitrary since elements of GF(2) don't have an
|
|
* ordering.
|
|
*/
|
|
|
|
int mzd_cmp(mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Copy matrix A to DST.
|
|
*
|
|
* \param DST May be NULL for automatic creation.
|
|
* \param A Source matrix.
|
|
*/
|
|
|
|
mzd_t *mzd_copy(mzd_t *DST, mzd_t const *A);
|
|
|
|
/**
|
|
* \brief Concatenate B to A and write the result to C.
|
|
*
|
|
* That is,
|
|
*
|
|
\verbatim
|
|
[ A ], [ B ] -> [ A B ] = C
|
|
\endverbatim
|
|
*
|
|
* The inputs are not modified but a new matrix is created.
|
|
*
|
|
* \param C Matrix, may be NULL for automatic creation
|
|
* \param A Matrix
|
|
* \param B Matrix
|
|
*
|
|
* \note This is sometimes called augment.
|
|
*/
|
|
|
|
mzd_t *mzd_concat(mzd_t *C, mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Stack A on top of B and write the result to C.
|
|
*
|
|
* That is,
|
|
*
|
|
\verbatim
|
|
[ A ], [ B ] -> [ A ] = C
|
|
[ B ]
|
|
\endverbatim
|
|
*
|
|
* The inputs are not modified but a new matrix is created.
|
|
*
|
|
* \param C Matrix, may be NULL for automatic creation
|
|
* \param A Matrix
|
|
* \param B Matrix
|
|
*/
|
|
|
|
mzd_t *mzd_stack(mzd_t *C, mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Copy a submatrix.
|
|
*
|
|
* Note that the upper bounds are not included.
|
|
*
|
|
* \param S Preallocated space for submatrix, may be NULL for automatic creation.
|
|
* \param M Matrix
|
|
* \param lowr start rows
|
|
* \param lowc start column
|
|
* \param highr stop row (this row is \em not included)
|
|
* \param highc stop column (this column is \em not included)
|
|
*/
|
|
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);
|
|
|
|
/**
|
|
* \brief Invert the matrix target using Gaussian elimination.
|
|
*
|
|
* To avoid recomputing the identity matrix over and over again, I may
|
|
* be passed in as identity parameter.
|
|
*
|
|
* \param INV Preallocated space for inversion matrix, may be NULL for automatic creation.
|
|
* \param A Matrix to be reduced.
|
|
* \param I Identity matrix.
|
|
*/
|
|
|
|
mzd_t *mzd_invert_naive(mzd_t *INV, mzd_t const *A, mzd_t const *I);
|
|
|
|
/**
|
|
* \brief Set C = A+B.
|
|
*
|
|
* C is also returned. If C is NULL then a new matrix is created which
|
|
* must be freed by mzd_free.
|
|
*
|
|
* \param C Preallocated sum matrix, may be NULL for automatic creation.
|
|
* \param A Matrix
|
|
* \param B Matrix
|
|
*/
|
|
|
|
mzd_t *mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Same as mzd_add but without any checks on the input.
|
|
*
|
|
* \param C Preallocated sum matrix, may be NULL for automatic creation.
|
|
* \param A Matrix
|
|
* \param B Matrix
|
|
*/
|
|
|
|
mzd_t *_mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);
|
|
|
|
/**
|
|
* \brief Same as mzd_add.
|
|
*
|
|
* \param C Preallocated difference matrix, may be NULL for automatic creation.
|
|
* \param A Matrix
|
|
* \param B Matrix
|
|
*/
|
|
|
|
#define mzd_sub mzd_add
|
|
|
|
/**
|
|
* \brief Same as mzd_sub but without any checks on the input.
|
|
*
|
|
* \param C Preallocated difference matrix, may be NULL for automatic creation.
|
|
* \param A Matrix
|
|
* \param B Matrix
|
|
*/
|
|
|
|
#define _mzd_sub _mzd_add
|
|
|
|
|
|
|
|
/**
|
|
* Get n bits starting a position (x,y) from the matrix M.
|
|
*
|
|
* \param M Source matrix.
|
|
* \param x Starting row.
|
|
* \param y Starting column.
|
|
* \param n Number of bits (<= m4ri_radix);
|
|
*/
|
|
|
|
static inline word mzd_read_bits(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
|
|
int const spot = y % m4ri_radix;
|
|
wi_t const block = y / m4ri_radix;
|
|
int const spill = spot + n - m4ri_radix;
|
|
word temp = (spill <= 0) ? M->rows[x][block] << -spill : (M->rows[x][block + 1] << (m4ri_radix - spill)) | (M->rows[x][block] >> spill);
|
|
return temp >> (m4ri_radix - n);
|
|
}
|
|
|
|
|
|
/**
|
|
* \brief a_row[a_startblock:] += b_row[b_startblock:] for offset 0
|
|
*
|
|
* Adds a_row of A, starting with a_startblock to the end, to
|
|
* b_row of B, starting with b_startblock to the end. This gets stored
|
|
* in A, in a_row, starting with a_startblock.
|
|
*
|
|
* \param A destination matrix
|
|
* \param a_row destination row for matrix C
|
|
* \param a_startblock starting block to work on in matrix C
|
|
* \param B source matrix
|
|
* \param b_row source row for matrix B
|
|
* \param b_startblock starting block to work on in matrix B
|
|
*
|
|
*/
|
|
|
|
static inline 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) {
|
|
|
|
wi_t wide = A->width - a_startblock - 1;
|
|
|
|
word *a = A->rows[a_row] + a_startblock;
|
|
word *b = B->rows[b_row] + b_startblock;
|
|
|
|
#if __M4RI_HAVE_SSE2
|
|
if(wide > 2) {
|
|
/** check alignments **/
|
|
if (__M4RI_ALIGNMENT(a,16)) {
|
|
*a++ ^= *b++;
|
|
wide--;
|
|
}
|
|
|
|
if (__M4RI_ALIGNMENT(a, 16) == 0 && __M4RI_ALIGNMENT(b, 16) == 0) {
|
|
__m128i *a128 = (__m128i*)a;
|
|
__m128i *b128 = (__m128i*)b;
|
|
const __m128i *eof = (__m128i*)((unsigned long)(a + wide) & ~0xFUL);
|
|
|
|
do {
|
|
*a128 = _mm_xor_si128(*a128, *b128);
|
|
++b128;
|
|
++a128;
|
|
} while(a128 < eof);
|
|
|
|
a = (word*)a128;
|
|
b = (word*)b128;
|
|
wide = ((sizeof(word) * wide) % 16) / sizeof(word);
|
|
}
|
|
}
|
|
#endif // __M4RI_HAVE_SSE2
|
|
|
|
if (wide > 0) {
|
|
wi_t n = (wide + 7) / 8;
|
|
switch (wide % 8) {
|
|
case 0: do { *(a++) ^= *(b++);
|
|
case 7: *(a++) ^= *(b++);
|
|
case 6: *(a++) ^= *(b++);
|
|
case 5: *(a++) ^= *(b++);
|
|
case 4: *(a++) ^= *(b++);
|
|
case 3: *(a++) ^= *(b++);
|
|
case 2: *(a++) ^= *(b++);
|
|
case 1: *(a++) ^= *(b++);
|
|
} while (--n > 0);
|
|
}
|
|
}
|
|
|
|
*a ^= *b & A->high_bitmask;
|
|
|
|
__M4RI_DD_MZD(A);
|
|
}
|
|
|
|
|
|
/**
|
|
* \brief c_row[c_startblock:] = a_row[a_startblock:] + b_row[b_startblock:] for offset 0
|
|
*
|
|
* Adds a_row of A, starting with a_startblock to the end, to
|
|
* b_row of B, starting with b_startblock to the end. This gets stored
|
|
* in C, in c_row, starting with c_startblock.
|
|
*
|
|
* \param C destination matrix
|
|
* \param c_row destination row for matrix C
|
|
* \param c_startblock starting block to work on in matrix C
|
|
* \param A source matrix
|
|
* \param a_row source row for matrix A
|
|
* \param a_startblock starting block to work on in matrix A
|
|
* \param B source matrix
|
|
* \param b_row source row for matrix B
|
|
* \param b_startblock starting block to work on in matrix B
|
|
*
|
|
*/
|
|
|
|
static inline 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) {
|
|
|
|
wi_t wide = A->width - a_startblock - 1;
|
|
word *a = A->rows[a_row] + a_startblock;
|
|
word *b = B->rows[b_row] + b_startblock;
|
|
word *c = C->rows[c_row] + c_startblock;
|
|
|
|
#if __M4RI_HAVE_SSE2
|
|
if(wide > 2) {
|
|
/** check alignments **/
|
|
if (__M4RI_ALIGNMENT(a,16)) {
|
|
*c++ = *b++ ^ *a++;
|
|
wide--;
|
|
}
|
|
|
|
if ( (__M4RI_ALIGNMENT(b, 16) | __M4RI_ALIGNMENT(c, 16)) == 0) {
|
|
__m128i *a128 = (__m128i*)a;
|
|
__m128i *b128 = (__m128i*)b;
|
|
__m128i *c128 = (__m128i*)c;
|
|
const __m128i *eof = (__m128i*)((unsigned long)(a + wide) & ~0xFUL);
|
|
|
|
do {
|
|
*c128 = _mm_xor_si128(*a128, *b128);
|
|
++c128;
|
|
++b128;
|
|
++a128;
|
|
} while(a128 < eof);
|
|
|
|
a = (word*)a128;
|
|
b = (word*)b128;
|
|
c = (word*)c128;
|
|
wide = ((sizeof(word) * wide) % 16) / sizeof(word);
|
|
}
|
|
}
|
|
#endif // __M4RI_HAVE_SSE2
|
|
|
|
if (wide > 0) {
|
|
wi_t n = (wide + 7) / 8;
|
|
switch (wide % 8) {
|
|
case 0: do { *(c++) = *(a++) ^ *(b++);
|
|
case 7: *(c++) = *(a++) ^ *(b++);
|
|
case 6: *(c++) = *(a++) ^ *(b++);
|
|
case 5: *(c++) = *(a++) ^ *(b++);
|
|
case 4: *(c++) = *(a++) ^ *(b++);
|
|
case 3: *(c++) = *(a++) ^ *(b++);
|
|
case 2: *(c++) = *(a++) ^ *(b++);
|
|
case 1: *(c++) = *(a++) ^ *(b++);
|
|
} while (--n > 0);
|
|
}
|
|
}
|
|
*c ^= ((*a ^ *b ^ *c) & C->high_bitmask);
|
|
|
|
__M4RI_DD_MZD(C);
|
|
}
|
|
|
|
|
|
/**
|
|
* \brief row3[col3:] = row1[col1:] + row2[col2:]
|
|
*
|
|
* Adds row1 of SC1, starting with startblock1 to the end, to
|
|
* row2 of SC2, starting with startblock2 to the end. This gets stored
|
|
* in DST, in row3, starting with startblock3.
|
|
*
|
|
* \param C destination matrix
|
|
* \param c_row destination row for matrix dst
|
|
* \param c_startblock starting block to work on in matrix dst
|
|
* \param A source matrix
|
|
* \param a_row source row for matrix sc1
|
|
* \param a_startblock starting block to work on in matrix sc1
|
|
* \param B source matrix
|
|
* \param b_row source row for matrix sc2
|
|
* \param b_startblock starting block to work on in matrix sc2
|
|
*
|
|
*/
|
|
static inline 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) {
|
|
|
|
if( (C == A) & (a_row == c_row) & (a_startblock == c_startblock) )
|
|
mzd_combine_even_in_place(C, c_row, c_startblock, B, b_row, b_startblock);
|
|
else
|
|
mzd_combine_even(C, c_row, c_startblock, A, a_row, a_startblock, B, b_row, b_startblock);
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* \brief Get n bits starting a position (x,y) from the matrix M.
|
|
*
|
|
* This function is in principle the same as mzd_read_bits,
|
|
* but it explicitely returns an 'int' and is used as
|
|
* index into an array (Gray code).
|
|
*/
|
|
|
|
static inline int mzd_read_bits_int(mzd_t const *M, rci_t const x, rci_t const y, int const n) {
|
|
return __M4RI_CONVERT_TO_INT(mzd_read_bits(M, x, y, n));
|
|
}
|
|
|
|
|
|
/**
|
|
* \brief Zero test for matrix.
|
|
*
|
|
* \param A Input matrix.
|
|
*
|
|
*/
|
|
int mzd_is_zero(mzd_t const *A);
|
|
|
|
/**
|
|
* \brief Clear the given row, but only begins at the column coloffset.
|
|
*
|
|
* \param M Matrix
|
|
* \param row Index of row
|
|
* \param coloffset Column offset
|
|
*/
|
|
|
|
void mzd_row_clear_offset(mzd_t *M, rci_t const row, rci_t const coloffset);
|
|
|
|
/**
|
|
* \brief Find the next nonzero entry in M starting at start_row and start_col.
|
|
*
|
|
* This function walks down rows in the inner loop and columns in the
|
|
* outer loop. If a nonzero entry is found this function returns 1 and
|
|
* zero otherwise.
|
|
*
|
|
* If and only if a nonzero entry is found r and c are updated.
|
|
*
|
|
* \param M Matrix
|
|
* \param start_row Index of row where to start search
|
|
* \param start_col Index of column where to start search
|
|
* \param r Row index updated if pivot is found
|
|
* \param c Column index updated if pivot is found
|
|
*/
|
|
|
|
int mzd_find_pivot(mzd_t const *M, rci_t start_row, rci_t start_col, rci_t *r, rci_t *c);
|
|
|
|
|
|
/**
|
|
* \brief Return the number of nonzero entries divided by nrows *
|
|
* ncols
|
|
*
|
|
* If res = 0 then 100 samples per row are made, if res > 0 the
|
|
* function takes res sized steps within each row (res = 1 uses every
|
|
* word).
|
|
*
|
|
* \param A Matrix
|
|
* \param res Resolution of sampling (in words)
|
|
*/
|
|
|
|
double mzd_density(mzd_t const *A, wi_t res);
|
|
|
|
/**
|
|
* \brief Return the number of nonzero entries divided by nrows *
|
|
* ncols considering only the submatrix starting at (r,c).
|
|
*
|
|
* If res = 0 then 100 samples per row are made, if res > 0 the
|
|
* function takes res sized steps within each row (res = 1 uses every
|
|
* word).
|
|
*
|
|
* \param A Matrix
|
|
* \param res Resolution of sampling (in words)
|
|
* \param r Row to start counting
|
|
* \param c Column to start counting
|
|
*/
|
|
|
|
double _mzd_density(mzd_t const *A, wi_t res, rci_t r, rci_t c);
|
|
|
|
|
|
/**
|
|
* \brief Return the first row with all zero entries.
|
|
*
|
|
* If no such row can be found returns nrows.
|
|
*
|
|
* \param A Matrix
|
|
*/
|
|
|
|
rci_t mzd_first_zero_row(mzd_t const *A);
|
|
|
|
/**
|
|
* \brief Return hash value for matrix.
|
|
*
|
|
* \param A Matrix
|
|
*/
|
|
|
|
static inline word mzd_hash(mzd_t const *A) {
|
|
word hash = 0;
|
|
for (rci_t r = 0; r < A->nrows; ++r)
|
|
hash ^= rotate_word(calculate_hash(A->rows[r], A->width), r % m4ri_radix);
|
|
return hash;
|
|
}
|
|
|
|
/**
|
|
* Return upper triangular submatrix of A
|
|
*
|
|
* \param U Output matrix, if NULL a new matrix will be returned
|
|
* \param A Source matrix
|
|
*
|
|
* \return U
|
|
*/
|
|
|
|
mzd_t *mzd_extract_u(mzd_t *U, mzd_t const *A);
|
|
|
|
/**
|
|
* Return lower triangular submatrix of A
|
|
*
|
|
* \param L Output matrix, if NULL a new matrix will be returned
|
|
* \param A Source matrix
|
|
*
|
|
* \return L
|
|
*/
|
|
|
|
mzd_t *mzd_extract_l(mzd_t *L, mzd_t const *A);
|
|
|
|
#endif // M4RI_MZD
|