mzd_slice.h

Go to the documentation of this file.

#ifndef M4RIE_MZD_SLICE
#define M4RIE_MZD_SLICE

/******************************************************************************
*
*            M4RIE: Linear Algebra over GF(2^e)
*
*    Copyright (C) 2010,2011 Martin Albrecht <martinralbrecht@googlemail.com>
*
*  Distributed under the terms of the GNU General Public License (GEL)
*  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/m4ri.h>
#include <m4rie/mzd_poly.h>
#include <m4rie/mzed.h>
#include <m4rie/blm.h>

typedef struct {
  mzd_t *x[16]; 
  rci_t nrows; 
  rci_t ncols; 
  unsigned int depth;   
  const gf2e *finite_field; 
} mzd_slice_t;


static inline mzd_slice_t *mzd_slice_init(const gf2e *ff, const rci_t m, const rci_t n) {
  mzd_slice_t *A = (mzd_slice_t*)m4ri_mm_malloc(sizeof(mzd_slice_t));

  A->finite_field = ff;
  A->nrows = m;
  A->ncols = n;
  A->depth = ff->degree;

  for(int i=0; i<A->depth; i++)
    A->x[i] = mzd_init(m,n);
  return A;
}

void mzd_slice_set_ui(mzd_slice_t *A, word value);

static inline mzd_slice_t *_mzd_slice_adapt_depth(mzd_slice_t *A, const unsigned int new_depth) {
  assert(A->finite_field->degree <= new_depth);

  if (new_depth < A->depth) {
    for(unsigned int i=new_depth; i<A->depth; i++) {
      mzd_free(A->x[i]);
      A->x[i] = NULL;
    }
  } else {
    for(unsigned int i=A->depth; i<new_depth; i++) {
      A->x[i] = mzd_init(A->nrows,A->ncols);
    }
  }
  A->depth = new_depth;
  return A;
}


static inline void mzd_slice_free(mzd_slice_t *A) {
  for(int i=0; i<A->depth; i++)
   mzd_free(A->x[i]);
#if __M4RI_USE_MM_MALLOC
  _mm_free(A);
#else
  free(A);
#endif
}

static inline mzd_slice_t *mzd_slice_copy(mzd_slice_t *B, const mzd_slice_t *A) {
  if(B == NULL)
    B = mzd_slice_init(A->finite_field, A->nrows, A->ncols);

  for(int i=0; i<A->depth; i++) {
    mzd_copy(B->x[i],A->x[i]);
  }
  return B;
}

static inline word mzd_slice_read_elem(const mzd_slice_t *A, const rci_t row, const rci_t col) {
  word ret = 0;
  for(int i=0; i<A->depth; i++) {
    ret |= mzd_read_bit(A->x[i], row, col)<<i;
  }
  return ret;
}

static inline void mzd_slice_add_elem(mzd_slice_t *A, const rci_t row, const rci_t col, word elem) {
  for(int i=0; i<A->depth; i++) {
    __mzd_xor_bits(A->x[i], row, col, 1, elem&1);
    elem=elem>>1;
  }
}

static inline void mzd_slice_write_elem(mzd_slice_t *A, const rci_t row, const rci_t col, word elem) {
  for(int i=0; i<A->depth; i++) {
    mzd_write_bit(A->x[i], row, col, elem&1);
    elem=elem>>1;
  }
}

static inline int mzd_slice_cmp(mzd_slice_t *A, mzd_slice_t *B) {
  int r = 0;
  if ((A->finite_field != B->finite_field) | (A->depth != B->depth) )
    return -1;
  for(int i=0; i<A->depth; i++)
    r |= mzd_cmp(A->x[i],B->x[i]);
  return r;
}

static inline int mzd_slice_is_zero(const mzd_slice_t *A) {
  for(int i=0; i<A->depth; i++) {
    if (!mzd_is_zero(A->x[i]))
      return 0;
  }
  return 1;
}

static inline void mzd_slice_row_swap(mzd_slice_t *A, const rci_t rowa, const rci_t rowb) {
  for(int i=0; i<A->depth; i++) {
    mzd_row_swap(A->x[i], rowa, rowb);
  }
}

static inline void mzd_slice_copy_row(mzd_slice_t* B, size_t i, const mzd_slice_t* A, size_t j) {
  for(int ii=0; ii<A->depth; ii++)
    mzd_copy_row(B->x[ii], i, A->x[ii], j);
}

static inline void mzd_slice_col_swap(mzd_slice_t *A, const rci_t cola, const rci_t colb) {
  for(int i=0; i<A->depth; i++)
    mzd_col_swap(A->x[i], cola, colb);
}

static inline void mzd_slice_col_swap_in_rows(mzd_slice_t *A, const rci_t cola, const rci_t colb, const rci_t start_row, rci_t stop_row) {
  for(unsigned int e=0; e < A->finite_field->degree; e++) {
    mzd_col_swap_in_rows(A->x[e], cola, colb, start_row, stop_row);
  };
}

static inline void mzd_slice_row_add(mzd_slice_t *A, const rci_t sourcerow, const rci_t destrow) {
  for(int i=0; i<A->depth; i++)
    mzd_row_add(A->x[i], sourcerow, destrow);
}

void mzd_slice_print(const mzd_slice_t *A);

static inline void _mzd_slice_compress_l(mzd_slice_t *A, const rci_t r1, const rci_t n1, const rci_t r2) {
  for(int i=0; i<A->depth; i++)
    _mzd_compress_l(A->x[i], r1, n1, r2);
}

static inline mzd_slice_t *mzd_slice_concat(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  if(C == NULL)
    C = mzd_slice_init(A->finite_field, A->nrows, A->ncols + B->ncols);

  for(int i=0; i<A->depth; i++) {
    mzd_concat(C->x[i], A->x[i], B->x[i]);
  }
  return C;
}

static inline mzd_slice_t *mzd_slice_stack(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  if(C == NULL)
    C = mzd_slice_init(A->finite_field, A->nrows + B->nrows, A->ncols);

  for(int i=0; i<A->depth; i++) {
    mzd_stack(C->x[i], A->x[i], B->x[i]);
  }
  return C;
}

static inline mzd_slice_t *mzd_slice_submatrix(mzd_slice_t *S, const mzd_slice_t *A,
                                               const size_t lowr, const size_t lowc, const size_t highr, const size_t highc) {
  if(S==NULL)
    S = mzd_slice_init(A->finite_field, highr - lowr, highc - lowc);

  for(int i=0; i<A->depth; i++) {
    mzd_submatrix(S->x[i], A->x[i], lowr, lowc, highr, highc);
  }
  return S;
}

static inline mzd_slice_t *mzd_slice_init_window(const mzd_slice_t *A,
                                                 const size_t lowr, const size_t lowc,
                                                 const size_t highr, const size_t highc) {
  mzd_slice_t *B = (mzd_slice_t *)m4ri_mm_malloc(sizeof(mzd_slice_t));
  B->finite_field = A->finite_field;
  B->depth = A->depth;
  B->nrows = highr - lowr;
  B->ncols = highc - lowc;
  for(int i=0; i<A->depth; i++) {
    B->x[i] = mzd_init_window(A->x[i], lowr, lowc, highr, highc);
  }
  return B;
}

static inline void mzd_slice_free_window(mzd_slice_t *A) {
  for(int i=0; i<A->depth; i++) {
    mzd_free_window(A->x[i]);
  }
  m4ri_mm_free(A);
}

static inline mzd_slice_t *_mzd_slice_add(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  _mzd_ptr_add(C->x, (const mzd_t**)A->x, (const mzd_t**)B->x, A->depth);
  return C;
}

static inline mzd_slice_t *mzd_slice_add(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  if ( (A->finite_field != B->finite_field) | (A->nrows != B->nrows) | (A->ncols != B->ncols) ) 
    m4ri_die("mzd_slice_add: input matrices A (%d x %d) and B (%d x %d) do not match.\n",A->nrows,A->ncols, B->nrows,B->ncols);

  if(C == NULL)
    C = mzd_slice_init(A->finite_field, A->nrows, A->ncols);
  else if ( (A->finite_field != C->finite_field) | (A->nrows != C->nrows) | (A->ncols != C->ncols) )
    m4ri_die("mzd_slice_add: input matrix A (%d x %d) and output matrix (%d x %d) do not match.\n",A->nrows,A->ncols, C->nrows, C->ncols);

  return _mzd_slice_add(C,A,B);
}

#define mzd_slice_sub mzd_slice_add

#define _mzd_slice_sub _mzd_slice_add

mzd_slice_t *_mzd_slice_addmul_naive(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B);

static inline mzd_slice_t *_mzd_slice_addmul_karatsuba(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  if (C == NULL)
    C = mzd_slice_init(A->finite_field, A->nrows, B->ncols);
  switch(A->finite_field->degree) {
  case  2:  _mzd_ptr_addmul_karatsuba2(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case  3:  _mzd_ptr_addmul_karatsuba3(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case  4:  _mzd_ptr_addmul_karatsuba4(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case  5:  _mzd_ptr_addmul_karatsuba5(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case  6:  _mzd_ptr_addmul_karatsuba6(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case  7:  _mzd_ptr_addmul_karatsuba7(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case  8:  _mzd_ptr_addmul_karatsuba8(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case  9:  _mzd_ptr_addmul_karatsuba9(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case 10: _mzd_ptr_addmul_karatsuba10(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case 11: _mzd_ptr_addmul_karatsuba11(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case 12: _mzd_ptr_addmul_karatsuba12(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case 13: _mzd_ptr_addmul_karatsuba13(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case 14: _mzd_ptr_addmul_karatsuba14(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case 15: _mzd_ptr_addmul_karatsuba15(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  case 16: _mzd_ptr_addmul_karatsuba16(A->finite_field, C->x, (const mzd_t**)A->x, (const mzd_t**)B->x); break;
  default:
    C = _mzd_slice_addmul_naive(C, A, B); break;
  }
  return C;
}

static inline mzd_slice_t *mzd_slice_mul_karatsuba(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  if (A->ncols != B->nrows || A->finite_field != B->finite_field)
    m4ri_die("mzd_slice_mul_karatsuba: rows, columns and fields must match.\n");
  if (C != NULL) {
    if (C->finite_field != A->finite_field || C->nrows != A->nrows || C->ncols != B->ncols)
      m4ri_die("mzd_slice_mul_karatsuba: rows and columns of returned matrix must match.\n");
    mzd_slice_set_ui(C,0);
  }
  return _mzd_slice_addmul_karatsuba(C, A, B);
}

static inline mzd_slice_t *mzd_slice_addmul_karatsuba(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  assert(C != NULL);
  if (A->ncols != B->nrows || A->finite_field != B->finite_field)
    m4ri_die("mzd_slice_addmul_karatsuba: rows, columns and fields must match.\n");
  if (C->finite_field != A->finite_field || C->nrows != A->nrows || C->ncols != B->ncols)
    m4ri_die("mzd_slice_addmul_karatsuba: rows and columns of returned matrix must match.\n");
  return _mzd_slice_addmul_karatsuba(C, A, B);
}

static inline mzd_slice_t *_mzd_slice_mul_blm(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B, blm_t *f) {  
  if (C == NULL)
    C = mzd_slice_init(A->finite_field, A->nrows, B->ncols);

  int free_f = 0;

  if (f == NULL) {
    const deg_t d = C->finite_field->degree;
    if (d > 16) 
      m4ri_die("degrees > 16 unsupported.\n");
    int *p = (int *)m4ri_mm_calloc(M4RIE_MAX_DEGREE+1, sizeof(int));
    p[d] = 1;
    free_f = 1;
    f = blm_init_crt(C->finite_field, d, d, p, 1);
    m4ri_mm_free(p);
  }

  _mzd_ptr_apply_blm(C->x, (const mzd_t**)A->x, (const mzd_t**)B->x, f);

  if (free_f)
    blm_free(f);
  
  return C;
}

static inline mzd_slice_t *mzd_slice_mul_blm(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B, blm_t *f) {
  if (A->ncols != B->nrows || A->finite_field != B->finite_field)
    m4ri_die("mzd_slice_mul_karatsuba: rows, columns and fields must match.\n");
  if (C != NULL) {
    if (C->finite_field != A->finite_field || C->nrows != A->nrows || C->ncols != B->ncols)
      m4ri_die("mzd_slice_mul_karatsuba: rows and columns of returned matrix must match.\n");
    mzd_slice_set_ui(C,0);
  }
  return _mzd_slice_mul_blm(C, A, B, f);
}

static inline mzd_slice_t *mzd_slice_addmul_blm(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B, blm_t *f) {
  assert(C != NULL);
  if (A->ncols != B->nrows || A->finite_field != B->finite_field)
    m4ri_die("mzd_slice_addmul_karatsuba: rows, columns and fields must match.\n");
  if (C->finite_field != A->finite_field || C->nrows != A->nrows || C->ncols != B->ncols)
    m4ri_die("mzd_slice_addmul_karatsuba: rows and columns of returned matrix must match.\n");
  mzd_slice_t *T = _mzd_slice_mul_blm(NULL, A, B, f);
  mzd_slice_add(C, C, T);
  mzd_slice_free(T);
  return C;
}


mzd_slice_t *mzd_slice_mul_scalar(mzd_slice_t *C, const word a, const mzd_slice_t *B);

mzd_slice_t *mzd_slice_addmul_scalar(mzd_slice_t *C, const word a, const mzd_slice_t *B);


static inline mzd_slice_t *mzd_slice_mul(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  return mzd_slice_mul_karatsuba(C,A,B);
}

static inline mzd_slice_t *mzd_slice_addmul(mzd_slice_t *C, const mzd_slice_t *A, const mzd_slice_t *B) {
  return mzd_slice_addmul_karatsuba(C,A,B);
}

static inline void mzd_slice_randomize(mzd_slice_t *A) {
  for(int i=0; i<A->depth; i++)
    mzd_randomize(A->x[i]);
}

#endif //M4RIE_MZD_SLICE