#include #include #include #include typedef unsigned int word; #define WSIZE (8*sizeof(word)) typedef struct pMatrix { size_t m, n, nsz; /* nsz = (n + WSIZE-1)/WSIZE */ word **data; } pMatrix; pMatrix *create_pMatrix(size_t m, size_t n); void free_pMatrix(pMatrix *pm); void pack_matrix(mxLogical *data, pMatrix *pm); void unpack_matrix(pMatrix *pm, mxLogical *col, mxLogical *data); void upper_triangularize(pMatrix *primal, mxLogical *col); void back_substitute(pMatrix *primal, pMatrix *dual, mxLogical *col); void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* function [dual,cols] = FormDualMod2(primal) * This function computes the dual of a matrix over GF(2). * If primal is a k by n matrix of full rank, dual will be * an n by (n-k) matrix of full rank with the property * primal * dual = 0 * Thus the row span of the primal and the column span of the dual are * "orthogonal" spaces. cols is set to a logical vector that indexes * the (n-k) columns of the dual that form an identity matrix. */ { pMatrix *primal, *dual; mxArray *cols; mxLogical *col; size_t k, n; if (nrhs != 1 || nlhs> 2) { mexErrMsgTxt("Usage: [dual, cols] = FormDualMod2(primal)"); } if (!mxIsLogical(prhs[0]) || mxGetNumberOfDimensions(prhs[0]) != 2) { mexErrMsgTxt("Usage: primal must be a logical matrix"); } k = mxGetM(prhs[0]); n = mxGetN(prhs[0]); if (k> n) { mxErrMsgTxt("Usage: primal must have fewer rows than columns"); } plhs[0] = mxCreateLogicalMatrix(n, n-k); /* convert boolean matrix to packed boolean */ if ((primal = create_pMatrix(k, n)) == NULL) { mxErrMsgTxt("Out of memory creating pMatrix"); } pack_matrix(mxGetLogicals(prhs[0]), primal); cols = mxCreateLogicalMatrix(1, n); col = mxGetLogicals(cols); upper_triangularize(primal, col); if ((dual = create_pMatrix(n, n)) == NULL) { free_pMatrix(primal); mxErrMsgTxt("Out of memory creating pMatrix"); } back_substitute(primal, dual, col); unpack_matrix(dual, col, mxGetLogicals(plhs[0])); free_pMatrix(primal); free_pMatrix(dual); if (nlhs == 2) plhs[1] = cols; else mxDestroyArray(cols); } pMatrix *create_pMatrix(size_t m, size_t n) /* create a packed matrix object on the heap */ { pMatrix *pm; size_t i; if ((pm = (pMatrix *)malloc(sizeof(pMatrix))) == NULL) return NULL; pm->m = m, pm->n = n, pm->nsz = (pm->n + WSIZE-1)/WSIZE; if ((pm->data = (word **)calloc(m, sizeof(word *))) == NULL) { free(pm); return NULL; } for (i = 0; i < m; ++i) { if ((pm->data[i] = (word *)calloc(n, sizeof(word *))) == NULL) { size_t j; for (j = 0; j < i; ++j) free(pm->data[j]); free(pm->data); free(pm); return NULL; } } return pm; } void free_pMatrix(pMatrix *pm) /* delete a packed matrix object */ { size_t i; for (i = 0; i < pm->m; ++i) free(pm->data[i]); free(pm->data); free(pm); } /* for pack & unpack, note that Matlab stores matrices by columns */ void pack_matrix(mxLogical *data, pMatrix *pm) /* converts a logical array into a packed boolean format */ { size_t i, j; word jm; for (j = 0, jm = 1; j < pm->n; ++j%WSIZE ? (jm <<= 1): (jm = 1)) for (i = 0; i < pm->m; ++i) if (*data++) pm->data[i][j/WSIZE] |= jm; } void unpack_matrix(pMatrix *pm, mxLogical *col, mxLogical *data) /* unpacks the columns specified by col into data */ { size_t i, j; word jm; for (j = 0, jm = 1; j < pm->n; ++j%WSIZE ? (jm <<= 1): (jm = 1)) if (!col || col[j]) for (i = 0; i < pm->m; ++i) *data++ = (pm->data[i][j/WSIZE] & jm) ? true: false; } void xor(word *src, word *dst, size_t n); void upper_triangularize(pMatrix *primal, mxLogical *col) /* Implements the first phase of Gaussian elimination. * col[j] = false if j is a column on the "diagonal" after making the * primal matrix "upper triangular" */ { size_t i, j, k; word jm; word **p, *q; p = primal->data; for (i = j = 0, jm = 1; i < primal->m && j < primal->n; ++j % WSIZE ? (jm <<= 1): (jm = 1)) { /* look for a pivot element */ for (k = i; k < primal->m; ++k) if (p[k][j/WSIZE] & jm) break; if (k == primal->m) { /* entire column is zero: this means that this column will * correspond to a free variable (or, in coding parlance, a * message bit) in the dual */ col[j] = true; continue; } col[j] = false; /* swap rows if we have to */ if (i != k) q = p[i], p[i] = p[k], p[k] = q; else k = i+1; /* eliminate the rest of the column */ while (k < primal->m) { if (p[k][j/WSIZE] & jm) xor(&p[i][j/WSIZE], &p[k][j/WSIZE], primal->nsz-j/WSIZE); ++k; } ++i; } /* we couldn't finish the process? */ if (i < primal->m) { mxErrMsgTxt("primal is not of full rank"); } while (j < primal->n) col[j++] = true; } void back_substitute(pMatrix *primal, pMatrix *dual, mxLogical *col) /* Implements the second phase. * dual is set to the identity in the rows specified by col, and the * rest is determined by back-substituting using primal. */ { size_t i, j, k; word **p, **d, jm, km; p = primal->data; d = dual->data; i = primal->m, j = primal->n; jm = 1 << (j%WSIZE); do { k = j, km = jm; j-- % WSIZE ? (jm>>= 1): (jm = 1u<<(wsize-1)); if (!col[j]) { i--; /* for the non-free variables, compute the row in the dual * required to satisfy the ith constraint in the primal, that is, * the equation (ith row of primal) * dual = 0. * The calculation below takes advantage of the fact that the bits * in the columns corresponding to message bits will multiply out * to themselves, since they are multiplied by the identity in the * dual. */ memcpy(&d[j][j/WSIZE], &p[i][j/WSIZE], /* this does d[i] = p[i] */ sizeof(word)*(primal->nsz-j/WSIZE)); /* this loop adds the rows from the dual corresponding to * non-message bits */ while (k < primal->n) { if (!col[k] && (p[i][k/WSIZE] & km)) xor(&d[k][k/WSIZE], &d[j][k/WSIZE], dual->nsz-k/WSIZE); ++k % WSIZE ? (km <<= 1): (km = 1); } } else { /* for the free variables, store the identity matrix */ dual->data[j][j/WSIZE] |= jm; } } while (j> 0); } void xor(word *src, word *dst, size_t n) /* add the src to the dst mod 2 */ { while (n--> 0) *dst++ ^= *src++; }

AltStyle によって変換されたページ (->オリジナル) /