packmat.cc.txt

压缩文件中是Error Correction Coding - Mathematical Methods and Algorithms(Wiley 2005)作者:(Todd K. Moon )的配

   int i, j;
   if(packmode != newpackmode) { // switching mode
	  for(i = 0; i < numpackrows; i++) {
		 for(j = 0; j < numpackcols; j++) {
			NewMat[j][i] = Matrix[i][j];
		 }
	  }
   }
   else {						// same mode
	  for(i = 0; i < numrows; i++) {
		 for(j = 0;j < numcols; j++) {
			setval(NewMat,j,i,GetVal(i,j));
		 }
	  }
   }
   setsize(T,numcols, numrows, newpackmode);
   if(T.Matrix) FREEMATRIX(T.Matrix);
   T.Matrix = NewMat;
}




void PackMat::clearfringe()
{
   int i,j;
   PKSIZE mask;
   if(packmode == HORIZPACK) {
	  mask = (1<<fracbits)-1;
	  for(i = 0; i < numrows; i++) {
		 Matrix[i][numpackcols1] &= mask;
	  }
   }
   else if(packmode == VERTPACK && fracbits) {
	  mask = (1<<fracbits)-1;
	  for(j = 0; j < numcols; j++) {
		 Matrix[numpackrows1][j] &= mask;
	  }
   }
}

int PackMat::getWeight() 		// get the weight of the entire matrix
{
   int i,j, wt;

   wt = 0;
   if(fracbits) {				// make sure there is no extra stuff on edges
	  clearfringe();
   }
   for(i = 0; i < numpackrows; i++) {
	  for(j = 0; j < numpackcols; j++) {
		 wt += getWeight(Matrix[i][j]);
	  }
   }
   return wt;
}

inline
int PackMat::getWeight(PKSIZE blob)
{
   int wt;
   int i;
   CTSIZE *ptr;					// point to the pieces of blob

   ptr = (CTSIZE *)&blob;
   wt = 0;
   for(i = 0; i < PkperCt; i++) {
	  wt += wtct[*ptr];
	  ptr++;					// point to the next bit
   }
   return wt;
}


// print function
ostream&
PackMat:: printpackmat(ostream &os) const 
{
   int i,j,k,i1;


   if(packmode == HORIZPACK) {
	  for(i = 0; i< numrows; i++) {
		 for(j = 0; j < numpackcols1; j++) {
			for(k = 0; k < bitsper; k++) {
			   if(Matrix[i][j] & (1<<k))
				  os << "1";
			   else
				  os << "0";
			   os << printspace;
			}
		 }
		 for(j = 0; j < numpackcols2; j++) {
			for(k = 0; k < fracbits; k++) {
			   if(Matrix[i][j+numpackcols1] & (1<<k))
				  os << "1";
			   else
				  os << "0";
			   os << printspace;
			}
		 }
		 os << endl;
	  }
   }
   else if(packmode == VERTPACK) {
	  for(i = 0; i < numpackrows1; i++) {
		 for(i1 = 0; i1 < bitsper; i1++) {
			for(j = 0; j < numcols; j++) {
			   if(Matrix[i][j] & (1<<i1)) 
				  os << "1";
			   else
				  os << "0";
			   os << printspace;
			}
			os << endl;
		 }
	  }
	  for(i = 0; i < numpackrows2; i++) {
		 for(i1 = 0; i1 < fracbits; i1++) {
			for(j = 0; j < numcols; j++) {
			   if(Matrix[i+numpackrows1][j] & (1<<i1)) 
				  os << "1";
			   else
				  os << "0";
			   os << printspace;
			}
			os << endl;
		 }
	  }
   }
   return os;
}


void PackMat::PrintInfo(void)
{
   
   if(packmode == HORIZPACK)
	  cout<<"Pack Mode = Column Pack"<<endl;
   else if(packmode == VERTPACK)
	  cout<<"PackMode = Row Pack"<<endl;
   cout<<"Rows = "<<numrows<<",  Cols = "<<numcols<<","<<endl;
   cout<<"Packrows: "<< numpackrows<< "  " << numpackrows1 << "  " <<
	  numpackrows2 << endl;
   cout<<"Packcols: "<< numpackcols<< "  " << numpackcols1 << "  " <<
	  numpackcols2 << endl;
   cout << "fracbits: " << fracbits << endl;

   cout << "bitsper=" << bitsper << "  PkperCt=" << PkperCt << 
	  "  shifter=" << shifter << endl;
}

void
PackMat::MakeRandom(double p)
{
   int i,j,k;
   unsigned int rno;
   PKSIZE blob;
   PKSIZE mask;

   if(p == 0 || p == 0.5) {	// if equal prob of 0 or 1
	  for(i = 0; i < numpackrows; i++) {
		 for(j = 0; j < numpackcols; j++) {
			blob = 0;
			for(k = 0; k < sizeof(PKSIZE); k++) {
			   rno = 0xFF&(rand()>>8); // pick off a byte
			   blob = (blob << 8) + rno;
			}
			Matrix[i][j] = blob;
		 }
	  }
	  // now make sure there is nothing left over 
	  if(fracbits) {
		 clearfringe();
	  }
   }
   else {						// else go by the probabilities
	  double rnod;
	  for(i = 0; i < numrows; i++) {
		 for(j = 0; j < numcols;  j++) {
			rnod = double(rand())/double(RAND_MAX);
			if(rnod < p)
			   SetVal(i,j,1);
			else
			   SetVal(i,j,0);
		 }
	  }
   }
}


void
PackMat::SetToZero()
{
   int i,j,k;
   unsigned int rno;
   PKSIZE blob;

   for(i = 0; i < numpackrows; i++) {
	  for(j = 0; j < numpackcols; j++) {
		 Matrix[i][j] = 0;
	  }
   }
}

void
PackMat::SetToOne()
{
   int i,j,k;
   unsigned int rno;
   PKSIZE blob;

   for(i = 0; i < numpackrows; i++) {
	  for(j = 0; j < numpackcols; j++) {
		 Matrix[i][j] = masker;
	  }
   }
   // now make sure there is nothing left over 
   clearfringe();
}



void PackMat::Add(PackMat &DM, PackMat &Sum)
{
// This function tests to make sure things are the right size
// and then calls the addition function

   if(packmode == DM.packmode) {
	  if(!(numcols == DM.numcols && numrows == DM.numrows)) {  
          // correct size to add
		 cerr << "Error: Matrices not same size in sum\n";
		 exit(-1);
	  }
	  checkandset(Sum,packmode,numrows,numcols);
	  Addnocheck(DM,Sum);
   }
   else {
	  cerr << "Addition not implemented for these matrix modes\n";
	  exit(-1);
   }
}


void PackMat::Addnocheck(PackMat &DM, PackMat &Sum)
{
// This function does no testing --- it assumes that
// everything is the right size and shape for addition

   int i,j;
   PKSIZE temp;
   int prodwt;

   for(i = 0; i < numpackrows; i++) {
	  for(j = 0; j < numpackcols; j++) {
		 Sum.Matrix[i][j] = Matrix[i][j] ^ DM.Matrix[i][j];
	  }
   }
}

void
PackMat::Identity(int n, PackMode pkmode)
{
   int i,j;

   checkandset(*this,pkmode,n,n);
   for(i = 0; i < numpackrows; i++) {
	  for(j = 0; j< numpackcols; j++) {
		 Matrix[i][j] = 0;
	  }
   }
   for(i = 0; i < n; i++) {
	  SetVal(i,i,1);
   }
}

int PackMat::IsIdentity(void)
{
   int i;
   if(numrows != numcols) return 0;
   int wt = getWeight();
   if(wt != numrows) return 0;
   for(i = 0; i < numrows; i++) {
	  if(GetVal(i,i) != 1) return 0;
   }
   return 1;
}


int
PackMat::ludcmp(int &numindout)
// returns 0 if matrix is singular, and 1 if not singular
// numindout is the number of linearly independent columns
{
   int i, j, k,k1,j1, jbig, jlit;
   PKSIZE tmp2,jlitmask;
   int tmp1;
   int pivotfound;

   if(packmode != HORIZPACK) {
	  cerr << "LU decomposition for horizontally packed data only\n";
	  exit(-1);
   }
   if(luindex) delete[] luindex;
   luindex = new int[numrows];
   for(j = 0; j < numrows; j++) {
	  luindex[j] = j;
   }

   U = *this;
   L.Size(numrows,numrows,HORIZPACK);
   L.SetToZero();

//   for(j = 0; j < numcols-1; j++) {		// loop over columns
   for(j = 0; j < numcols; j++) {		// loop over columns
	  jbig = j/bitsper;
	  jlit = j % bitsper;
	  jlitmask = (1<<jlit);
	  // determine a pivot row
	  pivotfound = 0;
	  for(k1 = j; k1 < numrows; k1++) {
		 if(U.Matrix[k1][jbig]&jlitmask) {// nonzero bit --- pivot on this
			pivotfound = 1;
			break;
		 }
	  }
	  if(!pivotfound) {
//		 numind = MIN(numrows,numcols);
		 numind = j;
		 numindout = numind;
// cout <<"Bailing: numind(1)=" << numind << "\n";
// cout << "L\n" << L << endl;
// cout << "U\n" << U << endl;
// for(j = 0; j < numrows; j++) {
//  L.SetVal(j,j,1);
//  }
//  PackMat Uc;
//  U.SetMode(VERTPACK,Uc);
//  PackMat LU;
//  L.Mult(Uc,LU);
// cout << "LU\n" << LU << endl;
		 return 0;
	  }
//cout << "j=" << j << "  pivotrow=" << k1 << endl;
	  if(j != k1) {				// necessary to swap rows j and k1
		 for(j1 = jbig; j1 < numpackcols; j1++) {
			tmp2 = U.Matrix[k1][j1];
			U.Matrix[k1][j1] = U.Matrix[j][j1];
			U.Matrix[j][j1] = tmp2;
		 }
		 for(j1 = 0; j1 <= jbig; j1++) {
			tmp2 = L.Matrix[k1][j1];
			L.Matrix[k1][j1] = L.Matrix[j][j1];
			L.Matrix[j][j1] = tmp2;
		 }

//		 cout << "swapped rows\n" << U << endl;
		 tmp1 = luindex[j];
		 luindex[j] = luindex[k1];
		 luindex[k1] = tmp1;
	  }
	  for(i = j+1; i < numrows; i++) {
		 if(U.Matrix[i][jbig]&jlitmask) {
//			cout << "l("<<i<<","<<j<<")=" << 1 <<endl;
			L.Matrix[i][jbig] |= (jlitmask | L.Matrix[i][jbig]);
			// row i <- row i - mult * row j
			for(k = jbig/bitsper; k < numpackcols; k++) {
			   U.Matrix[i][k] = U.Matrix[i][k] ^ U.Matrix[j][k];
			}
		 }
	  }
// cout << "j=" << j << endl << "U=\n" << U << endl;
// cout << "L=\n" << L << endl;

   }
//    for(j = 0; j < numrows; j++) {
// 	  L.SetVal(j,j,1);
//    }
//  PackMat Uc;
//  U.SetMode(VERTPACK,Uc);
//  PackMat LU;
//  L.Mult(Uc,LU);
// cout << "LU\n" << LU << endl;

   numindout = numind = MIN(numrows,numcols);
   return 1;
}	 


void
PackMat::lubacksub(PackVec &b,PackVec &y)
{  
   // make sure there is enough room for y
   y.checkandset(y,numrows);
   lubacksub(b.Vec,y.Vec);
}


void PackMat::lubacksub(PKSIZE *b,PKSIZE *y)
// b is a packed one-dimensional array representing a RHS
{
   int i,ibig,ilit,sum;
   int indx,indxbig, indxlit;
   int bi;
   for(i = 0; i < numpackcols; i++) y[i] = 0;
   // do the backsubstitution to solve Ly = Pb
   for(i = 0; i < numrows; i++) { 
	  ibig = i/bitsper;
	  ilit = i % bitsper;
	  indx = luindex[i];
	  indxbig = indx/bitsper;
	  indxlit = indx % bitsper;
	  bi = ((b[indxbig]&(1<<indxlit)))>>indxlit;
	  sum = (bi + innerprod(L.Matrix[i],y,i)) % 2;  
//cout << "i=" << i << "  bi=" << bi << "  y=" << sum << endl;
	  if(sum) {
		 y[ibig] |= (1<<ilit);
	  }
   }
//cout << endl;
//cout << endl;
   // do forward subsitution to solve Uy = x
   i = numrows-1; 				// do the last row by itself
   ibig = i/bitsper;
   ibig = i % bitsper;
   bi = (y[ibig] & (1<<ilit)) >> ilit;
   if(bi)
	  y[ibig] |= (1<<ilit);
   else
	  y[ibig] &= ~(1<<ilit);
//cout << "i=" << numrows-1 << "  bi=" << bi << "  y=" << int(y[ibig]) << endl;

   for(i = numrows-2; i>= 0; i--) {	// now do all the rest of the rows
//cout << "i=" << i << " ";
	  ibig = i/bitsper;
	  ilit = i % bitsper;
	  bi = (y[ibig] & (1<<ilit)) >> ilit;
	  sum = (bi + innerprod2(U.Matrix[i],y,i+1,numrows-1)) % 2;
//cout << "i=" << i << "  yi=" << bi << "  newy=" << sum << endl;
	  if(sum)
		 y[ibig] |= (1<<ilit);
	  else
		 y[ibig] &= ~(1<<ilit);
   }

}

int PackMat::innerprod2(PKSIZE *row, PKSIZE *col, int start, int end)
// do the inner product between two row-packed quantities, 
// ending at the last bit
{
   int stbig, stlit, ebig, elit, i, nl;
   int sum=0;
   PKSIZE mask,mask1;
   int bitsdone=0;

   stbig = start/bitsper;
   stlit = start % bitsper;
   ebig = end/bitsper;
   elit = end % bitsper;
   mask1 = ((1<<(bitsper-stlit))-1) << stlit;
   mask = mask1;

//   cout << "start=" << start << "  end=" << end << endl;

   for(i = stbig; i < ebig; i++) {
//cout<<"a: i=" << i << "   row=" << int(row[i]) << "  col=" << int(col[i]) << 
//  "  mask=" << int(mask) << "  ip=" << getWeight(row[i]&col[i]&mask) << endl;
	  sum += getWeight(row[i] & col[i] & mask);
	  mask = masker;
	  stlit = 0;
   }
   nl = elit-stlit+1;
   mask = ((1<<nl)-1)<<stlit;
//cout<<"b: i=" << i << "  mask=" << int(mask) << "   row=" << int(row[i]) 
//<<"  col="<<int(col[i]) <<  "  ip=" << getWeight(row[i]&col[i]&mask) << endl;
 sum += getWeight(row[i] & col[i] & mask);
   return sum % 2;
}

int PackMat::innerprod(PKSIZE *row, PKSIZE *col, int end)
// do the inner product between two row-packed quantities, starting
// at the first bit and covering end bits.
{
   int i, ebig, elit;
   int sum= 0;

   ebig = end/bitsper; elit = end % bitsper;
   for(i = 0; i < ebig; i++) {
	  sum += getWeight(row[i] & col[i]);
   }
   if(elit) {
	  sum += getWeight((row[i]&col[i]) & ((1<<elit)-1));
   }
//cout << "ip:  row=" << int(row[0]) << " col=" << int(col[0]) << "  sum=" << 
//   int(sum) << "  elit=" << elit << "  mask=" << ((1<<elit)-1) << endl;
   return sum % 2;
}



void 
PackMat::matinv(PackMat &Inv)
{
   int i,ibig,ilit,j;
   int numind;

   setsize(Inv,numrows,numcols,VERTPACK);
   callocmat(Inv);

   if(packmode != HORIZPACK) {
	  cerr << "Matrix must be horizontal packed";
	  exit(-1);
   }
   PKSIZE *b = new PKSIZE[numpackcols+1];
   PKSIZE *y = new PKSIZE[numpackcols+1];
   if(!ludcmp(numind)) {
	  cerr << "Cannot invert matrix\n";
	  exit(-1);
   }
   for(i = 0; i < numpackcols+1; i++) {
	  b[i] = y[i] = 0;
   }
   for(i = 0; i < numrows; i++) {
	  ibig = i/bitsper;
	  ilit = i % bitsper;
	  b[ibig] = (1<<ilit);
	  lubacksub(b,y);
	  // copy into the columns of the inverse
	  for(j = 0; j < numpackcols; j++) {
		 Inv.Matrix[j][i] = y[j];
	  }
	  b[ibig] = 0;
   }
}

void 
PackMat::reducetosystematic(PackMat &G, int &numind)
// Given a matrix, reduce to a systematic form.  The matrix
// returned may have row permutations from the original
// The number of linearly independent rows (the dimension of the code)
// is returned in numind
{
   int i,j,k;
   PKSIZE bi;
   int jbig, jlit;

   G = *this;
   G.ludcmp(numind);
//   cout << "numind=" << numind << endl;

   // now work with the upper matrix
   G = U;
//cout << "G=\n" << G << endl;
   for(j = numind-1; j > 0; j--) {	// column
	  jbig = j/bitsper;
	  jlit = j % bitsper;
	  for(i = j-1; i >= 0; i--) { // row
		 bi = G.Matrix[i][jbig] & (1<<jlit);
//cout << "i=" << i << "  j=" << j << "  " << int(bi) << 
//			"  G[i][jbig]=" << int(G.Matrix[i][jbig]) << endl;
		 if(bi) {				// row[i] <- row[i] + row[j]
			for(k = jbig; k < numpackcols; k++) {
//cout << "k=" << k << " ";
			   G.Matrix[i][k] ^= G.Matrix[j][k];
			}
//cout << endl;
		 }
	  }
//	  cout << "j=" << j << "G=\n" << G << endl;
   }
}

//**********************************************************************


void
PackVec::initPackVec(void )
{
   PackMat::initPackMat();
}

PackVec::PackVec(void)
{
   numrows = 0;
   numpackrows = 0;
   numpackrows1 = numpackrows2 = 0;
   fracbits = 0;
   Vec = 0;
}

inline void
PackVec::setsize(PackVec &A,int rows)
{
   int i,j;
   A.numrows = rows;

   A.numpackrows1 = A.numpackrows = A.numrows/PackMat::bitsper;
   A.numpackrows2 = 0;
   A.fracbits = A.numrows%PackMat::bitsper;
   if(A.fracbits) {
	  A.numpackrows++;			// increment to hold the fractional int
	  A.numpackrows2 = 1;
   }
}

inline void
PackVec::callocvec(PackVec &A)
{
   int i;
   A.Vec = new PKSIZE[A.numpackrows];
   // clear out the matrix
   for(i = 0; i < A.numpackrows; i++) {
	  A.Vec[i] = 0;
   }
}



PackVec::PackVec(int rows)
{
   int i,j;

   double temp1, temp2;

   setsize(*this,rows);
   callocvec(*this);
}


PackVec::PackVec(PackMat &mat)
{
   int i;

   if(mat.numrows==1) {			// copy the first row into the vector
	  setsize(*this,mat.numcols);
	  Vec = new PKSIZE[numpackrows];
	  for(i = 0; i < numpackrows; i++) {
		 Vec[i] = mat.Matrix[0][i];
	  }
   }
   else if(mat.numcols == 1) {	// copy the first column into the vector
	  setsize(*this, mat.numrows);
	  Vec = new PKSIZE[numpackrows];
	  for(i = 0; i < numpackrows; i++) {
		 Vec[i] = mat.Matrix[i][0];
	  }
   }
   else {						// ambiguous
	  cerr << "Cannot uniquely build vector from matrix\n";
	  exit(-1);
   }

}


PackVec::PackVec(char *fname)
{

   ifstream Densefile(fname);
   if(!Densefile) {
	  cerr<<"Error: Unable to open input file"<<fname<<endl;
	  exit(-1);
   }
   PackVec_is(Densefile);
   Densefile.close();
}

PackVec::PackVec(string desc)
{
   istringstream Densefile(desc);
   PackVec_is(Densefile);
}