src/mat.c

   1 /*
   2  * transsip - the telephony toolkit
   3  * By Daniel Borkmann <daniel@transsip.org>
   4  * Copyright 2012 Daniel Borkmann <dborkma@tik.ee.ethz.ch>
   5  * Subject to the GPL, version 2.
   6  * Based on Bhaskar Biswas and Nicolas Sendrier McEliece
   7  * implementation, LGPL 2.1.
   8  */
   9
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12 #include <string.h>
  13
  14 #include "mat.h"
  15 #include "xmalloc.h"
  16
  17 struct matrix *matrix_init(size_t rown, size_t coln)
  18 {
  19         struct matrix *A = NULL;
  20
  21         A = xzmalloc(sizeof(*A));
  22         A->coln = coln;
  23         A->rown = rown;
  24         A->rwdcnt = (1 + (coln - 1) / BITS_PER_LONG);
  25         A->alloc_size = rown * A->rwdcnt * sizeof(uint32_t);
  26         A->elem = xzmalloc(A->alloc_size);
  27
  28         return A;
  29 }
  30
  31 struct matrix *matrix_init_from_string(size_t rown, size_t coln,
  32                                        const unsigned char *str, size_t len)
  33 {
  34         struct matrix *A = NULL;
  35
  36         A = xzmalloc(sizeof(*A));
  37         A->coln = coln;
  38         A->rown = rown;
  39         A->rwdcnt = (1 + (coln - 1) / BITS_PER_LONG);
  40         A->alloc_size = rown * A->rwdcnt * sizeof(uint32_t);
  41         A->elem = xmemdup(str, len);
  42
  43         bug_on(A->alloc_size != len);
  44
  45         return A;
  46 }
  47
  48 void matrix_free(struct matrix *A)
  49 {
  50         xfree(A->elem);
  51         xfree(A);
  52 }
  53
  54 struct matrix *matrix_copy(struct matrix *A)
  55 {
  56         int i;
  57         struct matrix *X;
  58
  59         X = matrix_init(A->rown, A->coln);
  60
  61         for (i = 0; i < A->rwdcnt * A->rown; ++i)
  62                 X->elem[i] = A->elem[i];
  63
  64         return X;
  65 }
  66
  67 struct matrix *matrix_rowxor(struct matrix *A, int a, int b)
  68 {
  69         int i;
  70
  71         for (i = 0; i < A->rwdcnt; ++i)
  72                 A->elem[a * A->rwdcnt + i] ^= A->elem[b * A->rwdcnt + i];
  73
  74         return A;
  75 }
  76
  77 void matrix_vec_mul(uint32_t *cR, uint8_t *x, struct matrix *A)
  78 {
  79         int i,j;
  80         uint32_t *pt;
  81
  82         memset(cR, 0, A->rwdcnt * sizeof(uint32_t));
  83         pt = A->elem;
  84
  85         /* extract the first column in the form of char array */
  86         for (i = 0; i < A->rown; ++i) {
  87                 if ((x[i / 8] >> (i % 8)) & 1) {
  88                         for (j = 0; j < A->rwdcnt; ++j)
  89                                 cR[j] ^= *pt++;
  90                 } else {
  91                         pt += A->rwdcnt;
  92                 }
  93         }
  94 }
  95
  96 struct matrix *matrix_mul(struct matrix *A, struct matrix *B)
  97 {
  98         int i, j, k;
  99         struct matrix *C;
 100
 101         bug_on(A->coln != B->rown);
 102
 103         C = matrix_init(A->rown, B->coln);
 104         memset(C->elem, 0, C->alloc_size);
 105
 106         for (i = 0; i <A->rown; ++i)
 107                 for (j = 0; j < B->coln; ++j)
 108                         for (k = 0; k < A->coln; ++k)
 109                                 if (matrix_coeff(A, i, k) &&
 110                                     matrix_coeff(B, k, j))
 111                                         matrix_change_coeff(C, i, j);
 112         return C;
 113 }
 114
 115 uint32_t *matrix_rref(struct matrix *A)
 116 {
 117         /* the matrix is reduced from LSB...(from right) */
 118         uint32_t *perm;
 119         int i, j, failcnt, findrow, max;
 120
 121         max = A->coln - 1;
 122         perm = xmalloc(A->coln * sizeof(uint32_t));
 123
 124         for (i = 0; i < A->coln; ++i)
 125                 perm[i] = i; /* initialize permutation */
 126         failcnt = 0;
 127
 128         for (i = 0; i < A->rown; ++i,--max) {
 129                 findrow = 0;
 130
 131                 for (j = i; j < A->rown; ++j) {
 132                         if(matrix_coeff(A, j, max)) {
 133                                 /* not needed as ith row is 0 and jth row is 1 */
 134                                 if (i != j) /* xor to the row.(swap)? */
 135                                         A = matrix_rowxor(A, i, j);
 136                                 findrow = 1;
 137                                 break;
 138                         }
 139                 }
 140
 141                 /* if no row with a 1 found then swap last column and
 142                  * the column with no 1 down */
 143                 if (!findrow) {
 144                         perm[A->coln - A->rown - 1 - failcnt] = max;
 145                         failcnt++;
 146
 147                         if (!max)
 148                                 return NULL;
 149                         i--;
 150                 } else {
 151                         perm[i + A->coln - A->rown] = max;
 152                         /* fill the column downwards with 0's */
 153                         for (j = i + 1; j < A->rown; ++j) {
 154                                 if (matrix_coeff(A, j, max))
 155                                         /* check the arg. order */
 156                                         A = matrix_rowxor(A, j, i);
 157                         }
 158
 159                         /*  fill the column with 0's upwards too */
 160                         for (j = i - 1; j >= 0; --j) {
 161                                 if (matrix_coeff(A, j, max))
 162                                         A = matrix_rowxor(A, j, i);
 163                         }
 164                 }
 165         }
 166
 167         return perm;
 168 }