mzd.h

Go to the documentation of this file.

#ifndef M4RI_MZD
#define M4RI_MZD

/*******************************************************************
*
*                M4RI: Linear Algebra over GF(2)
*
*    Copyright (C) 2007, 2008 Gregory Bard <bard@fordham.edu>
*    Copyright (C) 2008-2010 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/
*
********************************************************************/

#include "m4ri_config.h"

#include <math.h>
#include <assert.h>
#include <stdio.h>

#if __M4RI_HAVE_SSE2
#include <emmintrin.h>
#endif

#include "misc.h"
#include "debug_dump.h"

#if __M4RI_HAVE_SSE2

#define __M4RI_SSE2_CUTOFF 10
#endif

#define __M4RI_MAX_MZD_BLOCKSIZE (((size_t)1) << 27)

#define __M4RI_MUL_BLOCKSIZE MIN(((int)sqrt((double)(4 * __M4RI_CPU_L2_CACHE))) / 2, 2048)

typedef struct {
  size_t size;
  word* begin;
  word* end;
} mzd_block_t;

typedef struct mzd_t {
  rci_t nrows;

  rci_t ncols;

  wi_t width; 

  wi_t rowstride;

  wi_t offset_vector;

  wi_t row_offset;

  uint16_t offset;

  uint8_t flags;

  uint8_t blockrows_log;

#if 0   // Commented out in order to keep the size of mzd_t 64 bytes (one cache line). This could be added back if rows was ever removed.

  int blockrows_mask;
#endif

  word high_bitmask;

  word low_bitmask;

  mzd_block_t *blocks;

  word **rows;

} mzd_t;

static wi_t const mzd_paddingwidth = 3;

static uint8_t const mzd_flag_nonzero_offset = 0x1;
static uint8_t const mzd_flag_nonzero_excess = 0x2;
static uint8_t const mzd_flag_windowed_zerooffset = 0x4;
static uint8_t const mzd_flag_windowed_zeroexcess = 0x8;
static uint8_t const mzd_flag_windowed_ownsblocks = 0x10;
static uint8_t const mzd_flag_multiple_blocks = 0x20;

static inline int mzd_is_windowed(mzd_t const *M) {
  return M->flags & (mzd_flag_nonzero_offset | mzd_flag_windowed_zerooffset);
}

static inline int mzd_owns_blocks(mzd_t const *M) {
  return M->blocks && (!mzd_is_windowed(M) || ((M->flags & mzd_flag_windowed_ownsblocks)));
}

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;
}

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;
}

static inline int mzd_row_to_block(mzd_t const* M, rci_t row) {
  return (M->row_offset + row) >> M->blockrows_log;
}

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;
}

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;
}

mzd_t *mzd_init(rci_t const r, rci_t const c);

void mzd_free(mzd_t *A);


mzd_t *mzd_init_window(mzd_t *M, rci_t const lowr, rci_t const lowc, rci_t const highr, rci_t const highc);

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);
}

#define mzd_free_window mzd_free

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;

  /* This is the case since we're only called from _mzd_ple_mmpf,
   * which makes the same assumption. Therefore we don't need
   * to take a mask_begin into account. */
  assert(M->offset == 0);

  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 = __M4RI_LEFT_BITMASK((M->ncols + M->offset) % m4ri_radix);

  if (width != 0) {
    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);
}

static inline void mzd_row_swap(mzd_t *M, rci_t const rowa, rci_t const rowb) {
  if(rowa == rowb)
    return;

  wi_t width = M->width - 1;
  word *a = M->rows[rowa];
  word *b = M->rows[rowb];
  word const mask_begin = __M4RI_RIGHT_BITMASK(m4ri_radix - M->offset);
  word const mask_end = __M4RI_LEFT_BITMASK((M->ncols + M->offset) % m4ri_radix);

  word tmp = (a[0] ^ b[0]) & mask_begin;
  if (width != 0) {
    a[0] ^= tmp;
    b[0] ^= tmp;
    
    for(wi_t i = 1; 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;
    
  } else {
    tmp &= mask_end;
    a[0] ^= tmp;
    b[0] ^= tmp;
  }

  __M4RI_DD_ROW(M, rowa);
  __M4RI_DD_ROW(M, rowb);
}

void mzd_copy_row(mzd_t *B, rci_t i, mzd_t const *A, rci_t j);

void mzd_col_swap(mzd_t *M, rci_t const cola, rci_t const colb);

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 + M->offset;
  rci_t const _colb = colb + M->offset;

  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_rows_in_block(M, 0), 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;
      }
      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;
      }
      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);
}

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+M->offset)/m4ri_radix], (col+M->offset) % m4ri_radix);
}

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 + M->offset) / m4ri_radix], (col + M->offset) % m4ri_radix, value);
}


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 + M->offset) % m4ri_radix;
  wi_t const block = (y + M->offset) / m4ri_radix;
  M->rows[x][block] ^= values << spot;
  int const space = m4ri_radix - spot;
  if (n > space)
    M->rows[x][block + 1] ^= values >> space;
}

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 + M->offset) % m4ri_radix;
  wi_t const block = (y + M->offset) / m4ri_radix;
  M->rows[x][block] &= values << spot;
  int const space = m4ri_radix - spot;
  if (n > space)
    M->rows[x][block + 1] &= values >> space;
}

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 + M->offset) % m4ri_radix;
  wi_t const block = (y + M->offset) / m4ri_radix;
  M->rows[x][block] &= ~(values << spot);
  int const space = m4ri_radix - spot;
  if (n > space)
    M->rows[x][block + 1] &= ~(values >> space);
}

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);
  coloffset += M->offset;
  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 = __M4RI_LEFT_BITMASK((M->ncols + M->offset) % m4ri_radix);

  *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);
}

void mzd_row_add(mzd_t *M, rci_t const sourcerow, rci_t const destrow);

mzd_t *mzd_transpose(mzd_t *DST, mzd_t const *A);

mzd_t *mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);

mzd_t *mzd_addmul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B);

mzd_t *_mzd_mul_naive(mzd_t *C, mzd_t const *A, mzd_t const *B, int const clear);

mzd_t *_mzd_mul_va(mzd_t *C, mzd_t const *v, mzd_t const *A, int const clear);

void mzd_randomize(mzd_t *M);

void mzd_set_ui(mzd_t *M, unsigned int const value);

rci_t mzd_gauss_delayed(mzd_t *M, rci_t const startcol, int const full);

rci_t mzd_echelonize_naive(mzd_t *M, int const full);

int mzd_equal(mzd_t const *A, mzd_t const *B);

int mzd_cmp(mzd_t const *A, mzd_t const *B);

mzd_t *mzd_copy(mzd_t *DST, mzd_t const *A);

mzd_t *mzd_concat(mzd_t *C, mzd_t const *A, mzd_t const *B);

mzd_t *mzd_stack(mzd_t *C, mzd_t const *A, mzd_t const *B);

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);

mzd_t *mzd_invert_naive(mzd_t *INV, mzd_t const *A, mzd_t const *I);

mzd_t *mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);

mzd_t *_mzd_add(mzd_t *C, mzd_t const *A, mzd_t const *B);

#define mzd_sub mzd_add

#define _mzd_sub _mzd_add



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 + M->offset) % m4ri_radix;
  wi_t const block = (y + M->offset) / 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);
}


void mzd_combine(mzd_t *DST, rci_t const row3, wi_t const startblock3,
                 mzd_t const *SC1, rci_t const row1, wi_t const startblock1, 
                 mzd_t const *SC2, rci_t const row2, wi_t const startblock2);


static inline void mzd_combine_weird(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) {
  word tmp;
  rci_t i = 0;


  for(; i + m4ri_radix <= A->ncols; i += m4ri_radix) {
    tmp = mzd_read_bits(A, a_row, i, m4ri_radix) ^ mzd_read_bits(B, b_row, i, m4ri_radix);
    mzd_clear_bits(C, c_row, i, m4ri_radix);
    mzd_xor_bits(C, c_row, i, m4ri_radix, tmp);
  }
  if(A->ncols - i) {
    tmp = mzd_read_bits(A, a_row, i, (A->ncols - i)) ^ mzd_read_bits(B, b_row, i, (B->ncols - i));
    mzd_clear_bits(C, c_row, i, (C->ncols - i));
    mzd_xor_bits(C, c_row, i, (C->ncols - i), tmp);
  }

  __M4RI_DD_MZD(C);
}

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 > __M4RI_SSE2_CUTOFF) {
    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 & __M4RI_LEFT_BITMASK(A->ncols%m4ri_radix);

  __M4RI_DD_MZD(A);
}


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;
  
  /* /\* this is a corner case triggered by Strassen multiplication */
  /*  * which assumes certain (virtual) matrix sizes  */
  /*  * 2011/03/07: I don't think this was ever correct *\/ */
  /* if (a_row >= A->nrows) { */
  /*   assert(a_row < A->nrows); */
  /*   for(wi_t i = 0; i < wide; ++i) { */
  /*     c[i] = b[i]; */
  /*   } */
  /* } else { */
#if __M4RI_HAVE_SSE2
  if(wide > __M4RI_SSE2_CUTOFF) {
    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) & __M4RI_LEFT_BITMASK(C->ncols%m4ri_radix));

  __M4RI_DD_MZD(C);
}


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));
}


int mzd_is_zero(mzd_t const *A);

void mzd_row_clear_offset(mzd_t *M, rci_t const row, rci_t const coloffset);

int mzd_find_pivot(mzd_t const *M, rci_t start_row, rci_t start_col, rci_t *r, rci_t *c);


double mzd_density(mzd_t const *A, wi_t res);

double _mzd_density(mzd_t const *A, wi_t res, rci_t r, rci_t c);


rci_t mzd_first_zero_row(mzd_t const *A);

static inline unsigned long long mzd_hash(mzd_t const *A) {
  unsigned long long 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;
}

mzd_t *mzd_extract_u(mzd_t *U, mzd_t const *A);

mzd_t *mzd_extract_l(mzd_t *L, mzd_t const *A);

#endif // M4RI_MZD