/* 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 the field base */

  q = (int)mxGetScalar(Q_IN);



/* read sparse matrix H */

    sr1  = mxGetPr(H_IN);

    irs1 = mxGetIr(H_IN);  /* row */

    jcs1 = mxGetJc(H_IN);  /* column */

    nz = mxGetNzmax(H_IN); /* number of nonzero elements (they are ones)*/

    M = mxGetM(H_IN);   

    N = mxGetN(H_IN);



	

/* create working array HH[row][column]*/

    HH = (unsigned char **)mxMalloc(M*sizeof(unsigned char *));

    for(i=0 ; i<M ; i++){

      HH[i] = (unsigned char *)mxMalloc(N*sizeof(unsigned char));

    }

    for(i=0 ; i<M ; i++)

      for(j=0 ; j<N ; j++)

	    HH[i][j] = 0; /* initialize all to zero */



    k=0;

    for(j=0 ; j<N ; j++) {

      for(i=0 ; i<(jcs1[j+1]-jcs1[j]) ; i++) {

        ii = irs1[k]; /* index in column j*/ 

        HH[ii][j] = (unsigned char)sr1[k]; /* put  nonzeros */

        k++;

      }

    }



/* invert HH matrix here */

	/* row and column indices */

	ir = (int *)mxMalloc(M*sizeof(int));

        jc = (int *)mxMalloc(N*sizeof(int));

	for( i=0 ; i<M ; i++)

		ir[i] = i;

	for( j=0 ; j<N ; j++)

		jc[j] = j;







    /* perform Gaussian elimination on H, store reodering operations */

    rdep = 0; /* number of dependent rows in H*/

    d = 0;    /* current diagonal element */



    while( (d+rdep) < M) { /* cycle through independent rows of H */

    

		j = d; /* current column index along row ir[d] */

		while( (HH[ir[d]][jc[j]] == 0) && (j<(N-1)) )

			j++;            /* find first nonzero element in row i */

		if( HH[ir[d]][jc[j]] ) { /* found nonzero element */





            /* swap columns */

			tmp = jc[d]; jc[d] = jc[j]; jc[j] = tmp;



			if(q==2) { /* GF2 */

			  /* eliminate current column using row operations */

		    	for(ii=0 ; ii<M ; ii++) 

			     if(HH[ir[ii]][jc[d]] && (ii != d)) /* nonzero and non-diagonal */

			    	for(jj=d ; jj<N ; jj++) 

				     HH[ir[ii]][jc[jj]] = (HH[ir[ii]][jc[jj]]+HH[ir[d]][jc[jj]])%2;

			}

			else { /* GFq */

			  	

				scale = GFq_inv(HH[ir[d]][jc[d]],q); /* inverse of the diag. element */

                /* scale the current row to make the first element 1 */

				for(jj=0 ; jj<N ; jj++) 

					HH[ir[d]][jc[jj]] = GFq_m(HH[ir[d]][jc[jj]],scale,q);



				/* eliminate current column using row operations */

				for(ii=0 ; ii<M ; ii++) {

					if(HH[ir[ii]][jc[d]] && (ii != d)) {

						scale = HH[ir[ii]][jc[d]];

						for(jj=d ; jj<N ; jj++) {

							tmp = GFq_m(HH[ir[d]][jc[jj]],scale,q);

							HH[ir[ii]][jc[jj]] = GFq_a(HH[ir[ii]][jc[jj]],tmp);

						}

					}

				}

			}

		}

		else { /* all zeros -  need to delete this row and update indices */

            rdep++; /* increase number of dependent rows */

			tmp = ir[d];

			ir[d] = ir[M-rdep];

			ir[M-rdep] = tmp;

			d--; /* no diagonal element is found */

		}

		d++; /* increase the number of diagonal elements */

  

	}/*while i+rdep*/

/* done inverting HH */





    K = N-M+rdep; /* true K */



/* create G matrix  G = [A'| I] if H = [I|A]*/

	GG = (unsigned char **)mxMalloc(K*sizeof(unsigned char *));

    for(i=0 ; i<K ; i++){

        GG[i] = (unsigned char *)mxMalloc(N*sizeof(unsigned char));

	}

    for(i=0 ; i<K ; i++)

		for(j=0 ; j<(N-K) ; j++) {

			tmp = (N-K+i);

			GG[i][j] = HH[ir[j]][jc[tmp]];

		}



    for(i=0 ; i<K ; i++)

		for(j=(N-K); j<N ; j++)

			if(i == (j-N+K) ) /* diagonal */

			    GG[i][j] = 1;

			else

				GG[i][j] = 0;



/* NOTE, it is a very inefficient way to store G. Change to taste!*/

    G_OUT = mxCreateDoubleMatrix(K, N, mxREAL);

    /* Assign pointers to the output matrix */

    g = mxGetPr(G_OUT);

    for(i=0 ; i<K ; i++)

		for(j=0 ; j<N; j++)

            g[i+j*K] = GG[i][j];





    H_OUT = mxCreateSparse(M,N,nz,mxREAL);

    sr2  = mxGetPr(H_OUT);

    irs2 = mxGetIr(H_OUT);  /* row */

    jcs2 = mxGetJc(H_OUT);  /* column */

    /* Write H_OUT swapping columns according to jc */

    k = 0; 

    for (j=0; (j<N ); j++) {

	  jcs2[j] = k;

      tmp = jcs1[jc[j]+1]-jcs1[jc[j]];

	  for (i=0; i<tmp ; i++) {

        kk = jcs1[jc[j]]+i;

		sr2[k] = sr1[kk];

		irs2[k] = irs1[kk];

		k++;	    

	  }

    }

    jcs2[N] = k;

    



/* free the memory */

  for( j=0 ; j<M ; j++) {

    mxFree(HH[j]);

  }

  mxFree(HH);

  mxFree(ir);

  mxFree(jc);  

  for(i=0;i<K;i++){

    mxFree(GG[i]);

  }

  mxFree(GG);

  return;

}