makay的LDPC码程序

/* Invert sparse binary H for LDPC*/ /* Author : Igor Kozintsev igor@ifp.uiuc.edu Please let me know if you find bugs in this code (I did test it but I still have some doubts). All other comments are welcome too :) ! I use a simple algorithm to invert H. We convert H to [I | A] [junk ] using column reodering and row operations (junk - a few rows of H which are linearly dependent on the previous ones) G is then found as G = [A'|I] G is stored as array of doubles in Matlab which is very inefficient. Internal representation in this programm is unsigned char. Please modify the part which writes G if you wish. */ #include <math.h> #include "mex.h" /* Input Arguments: tentative H matrix*/ #define H_IN prhs[0] #define Q_IN prhs[1] /* field base */ /* Output Arguments: final matrices*/ #define H_OUT plhs[0] #define G_OUT plhs[1] /************************************ GFq math *******************************/ /* This file contains lookup tables and routines required * to perform the field operations over GF(q), i.e. addition * and multiplication. * * Addition is easy, we use exclusive-or operation. * For multiplication we need two tables for each q. * The first is the logarithm table, and the second * is the exponential table. We catch multiplication * by zero and one separately. * * Source for tables: MacWilliams and Sloane, fig 4.5 * * WARNING: for speed, no check is made that legal values * are supplied. */ const int log4[4] = {0,0,1,2}; /* stores i at address \alpha^i */ const int log8[8] = {0,0,1,3,2,6,4,5}; const int log16[16] = {0,0,1,4,2,8,5,10,3,14,9,7,6,13,11,12}; const int log32[32] = {0,0,1,18,2,5,19,11,3,29,6,27,20,8,12,23,4,\ 10,30,17,7,22,28,26,21,25,9,16,13,14,24,15}; const int log64[64] = {0,0,1,6,2,12,7,26,3,32,13,35,8,48,27,18,4,24,\ 33,16,14,52,36,54,9,45,49,38,28,41,19,56,5,62,\ 25,11,34,31,17,47,15,23,53,51,37,44,55,40,10,\ 61,46,30,50,22,39,43,29,60,42,21,20,59,57,58}; const int log128[128] = {0,0,1,31,2,62,32,103,3,7,63,15,33,84,104,\ 93, 4,124,8,121,64,79,16,115,34,11,85,38,105,46,94,51,\ 5,82,125,60,9,44,122,77,65,67,80,42,17,69,116,23,35,118,\ 12,28,86,25,39,57,106,19,47,89,95,71,52,110,6,14,83,92,126,\ 30,61,102,10,37,45,50,123,120,78,114,66,41,68,22,81,59,43,76,\ 18,88,70,109,117,27,24,56,36,49,119,113,13,91,29,101,87,108,\ 26,55,40,21,58,75,107,54,20,74,48,112,90,100,96,97,72,98,53,73,111,99}; const int log256[256] = {0,0,1,25,2,50,26,198,3,223,51,238,27,104,199,75,4,100,\ 224,14,52,141,239,129,28,193,105,248,200,8,76,113,5,138,101,47,225,\ 36,15,33,53,147,142,218,240,18,130,69,29,181,194,125,106,39,249,185,\ 201,154,9,120,77,228,114,166,6,191,139,98,102,221,48,253,226,152,37,\ 179,16,145,34,136,54,208,148,206,143,150,219,189,241,210,19,92,131,\ 56,70,64,30,66,182,163,195,72,126,110,107,58,40,84,250,133,186,61,202,\ 94,155,159,10,21,121,43,78,212,229,172,115,243,167,87,7,112,192,247,\ 140,128,99,13,103,74,222,237,49,197,254,24,227,165,153,119,38,184,180,\ 124,17,68,146,217,35,32,137,46,55,63,209,91,149,188,207,205,144,135,151,\ 178,220,252,190,97,242,86,211,171,20,42,93,158,132,60,57,83,71,109,65,\ 162,31,45,67,216,183,123,164,118,196,23,73,236,127,12,111,246,108,161,59,\ 82,41,157,85,170,251,96,134,177,187,204,62,90,203,89,95,176,156,169,160,\ 81,11,245,22,235,122,117,44,215,79,174,213,233,230,231,173,232,116,214,\ 244,234,168,80,88,175}; const int exp4[3] = {1,2,3}; /* stores \alpha^i at address i */ const int exp8[7] = {1,2,4,3,6,7,5}; const int exp16[15] = {1,2,4,8,3,6,12,11,5,10,7,14,15,13,9}; const int exp32[31] = {1,2,4,8,16,5,10,20,13,26,17,7,14,28,29,31,\ 27,19,3,6,12,24,21,15,30,25,23,11,22,9,18}; const int exp64[63] = {1,2,4,8,16,32,3,6,12,24,48,35,5,10,20,40,19,\ 38,15,30,60,59,53,41,17,34,7,14,28,56,51,37,\ 9,18,36,11,22,44,27,54,47,29,58,55,45,25,50,\ 39,13,26,52,43,21,42,23,46,31,62,63,61,57,49,33}; const int exp128[127] = {1,2,4,8,16,32,64,9,18,36,72,25,50,100,65,11,\ 22,44,88,57,114,109,83,47,94,53,106,93,51,102,69,3,6,12,24,\ 48,96,73,27,54,108,81,43,86,37,74,29,58,116,97,75,31,62,124,\ 113,107,95,55,110,85,35,70,5,10,20,40,80,41,82,45,90,61,122,\ 125,115,111,87,39,78,21,42,84,33,66,13,26,52,104,89,59,118,101,\ 67,15,30,60,120,121,123,127,119,103,71,7,14,28,56,112,105,91,63,\ 126,117,99,79,23,46,92,49,98,77,19,38,76,17,34,68}; const int exp256[255] = {1,2,4,8,16,32,64,128,29,58,116,232,205,135,19,38,76,\ 152,45,90,180,117,234,201,143,3,6,12,24,48,96,192,157,39,78,156,\ 37,74,148,53,106,212,181,119,238,193,159,35,70,140,5,10,20,40,80,\ 160,93,186,105,210,185,111,222,161,95,190,97,194,153,47,94,188,101,\ 202,137,15,30,60,120,240,253,231,211,187,107,214,177,127,254,\ 225,223,163,91,182,113,226,217,175,67,134,17,34,68,136,13,26,52,104,\ 208,189,103,206,129,31,62,124,248,237,199,147,59,118,236,197,151,51,\ 102,204,133,23,46,92,184,109,218,169,79,158,33,66,132,21,42,84,168,\ 77,154,41,82,164,85,170,73,146,57,114,228,213,183,115,230,209,191,99,\ 198,145,63,126,252,229,215,179,123,246,241,255,227,219,171,75,150,49,\ 98,196,149,55,110,220,165,87,174,65,130,25,50,100,200,141,7,14,28,56,\ 112,224,221,167,83,166,81,162,89,178,121,242,249,239,195,155,43,86,172,\ 69,138,9,18,36,72,144,61,122,244,245,247,243,251,235,203,139,11,22,44,\ 88,176,125,250,233,207,131,27,54,108,216,173,71,142}; /* For testing: main(){ u_int i,j,q; while(1){ printf("please enter a b q:"); scanf("%u %u %u",&i,&j,&q); printf("\n a * b in GF(q) = %u\n",GFq_m(i,j,q)); } } */ int GFq_m(int a, int b, int q) { if ( a == 0 || b == 0 ) return 0 ; if ( a == 1 ) return b ; if ( b == 1 ) return a ; switch (q){ case 256: return exp256[(log256[a]+log256[b])%255]; case 128: return exp128[(log128[a]+log128[b])%127]; case 64: return exp64[(log64[a]+log64[b])%63]; case 32: return exp32[(log32[a]+log32[b])%31]; case 16: return exp16[(log16[a]+log16[b])%15]; case 8: return exp8[(log8[a]+log8[b])%7]; case 4: return exp4[(log4[a]+log4[b])%3]; } mexErrMsgTxt(1,"GFq_m: I'm afraid I don't know how to multiply in GFq\n"); return 0 ; } int GFq_inv(int a, int q) { if ( a == 0) mexErrMsgTxt(1,"GFq_inv: no inverse for 0!\n");; if ( a == 1 ) return 1 ; switch (q){ case 256: return exp256[(255-log256[a])]; case 128: return exp128[(127-log128[a])]; case 64: return exp64[(63-log64[a])]; case 32: return exp32[(31-log32[a])]; case 16: return exp16[(15-log16[a])]; case 8: return exp8[(7-log8[a])]; case 4: return exp4[(3-log4[a])]; } mexErrMsgTxt(1,"GFq_inv: not defined inverse for GFq\n"); return 0 ; } int GFq_a(int a, int b) { return a^b; } /************************************ end GFq math *******************************/ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) { unsigned char **HH, **GG; int ii, jj, *ir, *jc, rdep, tmp, d, q, scale; double *sr1, *sr2, *g; int N,M,K,i,j,k,kk,nz,*irs1,*jcs1, *irs2, *jcs2; /* Check for proper number of arguments */ if (nrhs != 2) { mexErrMsgTxt("h2g requires two input arguments."); } else if (nlhs != 2) { mexErrMsgTxt("h2g requires two output arguments."); } else if (!mxIsSparse(H_IN)) { mexErrMsgTxt("h2g requires sparse H matrix."); } /* get th