00001 00030 #include <itpp/base/gf2mat.h> 00031 #include <itpp/base/specmat.h> 00032 #include <itpp/base/matfunc.h> 00033 #include <itpp/base/converters.h> 00034 #include <iostream> 00035 00036 namespace itpp { 00037 00038 // ====== IMPLEMENTATION OF THE ALIST CLASS ========== 00039 00040 GF2mat_sparse_alist::GF2mat_sparse_alist(const std::string &fname) 00041 : data_ok(false) 00042 { 00043 read(fname); 00044 } 00045 00046 void GF2mat_sparse_alist::read(const std::string &fname) 00047 { 00048 std::ifstream file; 00049 std::string line; 00050 std::stringstream ss; 00051 00052 file.open(fname.c_str()); 00053 it_assert(file.is_open(), 00054 "GF2mat_sparse_alist::read(): Could not open file \"" 00055 << fname << "\" for reading"); 00056 00057 // parse N and M: 00058 getline(file, line); 00059 ss << line; 00060 ss >> N >> M; 00061 it_assert(!ss.fail(), 00062 "GF2mat_sparse_alist::read(): Wrong alist data (N or M)"); 00063 it_assert((N > 0) && (M > 0), 00064 "GF2mat_sparse_alist::read(): Wrong alist data"); 00065 ss.seekg(0, std::ios::end); 00066 ss.clear(); 00067 00068 // parse max_num_n and max_num_m 00069 getline(file, line); 00070 ss << line; 00071 ss >> max_num_n >> max_num_m; 00072 it_assert(!ss.fail(), 00073 "GF2mat_sparse_alist::read(): Wrong alist data (max_num_{n,m})"); 00074 it_assert((max_num_n >= 0) && (max_num_n <= N) && 00075 (max_num_m >= 0) && (max_num_m <= M), 00076 "GF2mat_sparse_alist::read(): Wrong alist data"); 00077 ss.seekg(0, std::ios::end); 00078 ss.clear(); 00079 00080 // parse weight of each column n 00081 num_nlist.set_size(N); 00082 num_nlist.clear(); 00083 getline(file, line); 00084 ss << line; 00085 for (int i = 0; i < N; i++) { 00086 ss >> num_nlist(i); 00087 it_assert(!ss.fail(), 00088 "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist(" 00089 << i << "))"); 00090 it_assert((num_nlist(i) >= 0) && (num_nlist(i) <= M), 00091 "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist(" 00092 << i << "))"); 00093 } 00094 ss.seekg(0, std::ios::end); 00095 ss.clear(); 00096 00097 // parse weight of each row m 00098 num_mlist.set_size(M); 00099 num_mlist.clear(); 00100 getline(file, line); 00101 ss << line; 00102 for (int i = 0; i < M; i++) { 00103 ss >> num_mlist(i); 00104 it_assert(!ss.fail(), 00105 "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist(" 00106 << i << "))"); 00107 it_assert((num_mlist(i) >= 0) && (num_mlist(i) <= N), 00108 "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist(" 00109 << i << "))"); 00110 } 00111 ss.seekg(0, std::ios::end); 00112 ss.clear(); 00113 00114 // parse coordinates in the n direction with non-zero entries 00115 nlist.set_size(N, max_num_n); 00116 nlist.clear(); 00117 for (int i = 0; i < N; i++) { 00118 getline(file, line); 00119 ss << line; 00120 for (int j = 0; j < num_nlist(i); j++) { 00121 ss >> nlist(i, j); 00122 it_assert(!ss.fail(), 00123 "GF2mat_sparse_alist::read(): Wrong alist data (nlist(" 00124 << i << "," << j << "))"); 00125 it_assert((nlist(i, j) >= 0) && (nlist(i, j) <= M), 00126 "GF2mat_sparse_alist::read(): Wrong alist data (nlist(" 00127 << i << "," << j << "))"); 00128 } 00129 ss.seekg(0, std::ios::end); 00130 ss.clear(); 00131 } 00132 00133 // parse coordinates in the m direction with non-zero entries 00134 mlist.set_size(M, max_num_m); 00135 mlist.clear(); 00136 for (int i = 0; i < M; i++) { 00137 getline(file, line); 00138 ss << line; 00139 for (int j = 0; j < num_mlist(i); j++) { 00140 ss >> mlist(i, j); 00141 it_assert(!ss.fail(), 00142 "GF2mat_sparse_alist::read(): Wrong alist data (mlist(" 00143 << i << "," << j << "))"); 00144 it_assert((mlist(i, j) >= 0) && (mlist(i, j) <= N), 00145 "GF2mat_sparse_alist::read(): Wrong alist data (mlist(" 00146 << i << "," << j << "))"); 00147 } 00148 ss.seekg(0, std::ios::end); 00149 ss.clear(); 00150 } 00151 00152 file.close(); 00153 data_ok = true; 00154 } 00155 00156 void GF2mat_sparse_alist::write(const std::string &fname) const 00157 { 00158 it_assert(data_ok, 00159 "GF2mat_sparse_alist::write(): alist data not ready for writing"); 00160 00161 std::ofstream file(fname.c_str(), std::ofstream::out); 00162 it_assert(file.is_open(), 00163 "GF2mat_sparse_alist::write(): Could not open file \"" 00164 << fname << "\" for writing"); 00165 00166 file << N << " " << M << std::endl; 00167 file << max_num_n << " " << max_num_m << std::endl; 00168 00169 for (int i = 0; i < num_nlist.length() - 1; i++) 00170 file << num_nlist(i) << " "; 00171 file << num_nlist(num_nlist.length() - 1) << std::endl; 00172 00173 for (int i = 0; i < num_mlist.length() - 1; i++) 00174 file << num_mlist(i) << " "; 00175 file << num_mlist(num_mlist.length() - 1) << std::endl; 00176 00177 for (int i = 0; i < N; i++) { 00178 for (int j = 0; j < num_nlist(i) - 1; j++) 00179 file << nlist(i, j) << " "; 00180 file << nlist(i, num_nlist(i) - 1) << std::endl; 00181 } 00182 00183 for (int i = 0; i < M; i++) { 00184 for (int j = 0; j < num_mlist(i) - 1; j++) 00185 file << mlist(i, j) << " "; 00186 file << mlist(i, num_mlist(i) - 1) << std::endl; 00187 } 00188 00189 file.close(); 00190 } 00191 00192 00193 GF2mat_sparse GF2mat_sparse_alist::to_sparse(bool transpose) const 00194 { 00195 GF2mat_sparse sbmat(M, N, max_num_m); 00196 00197 for (int i = 0; i < M; i++) { 00198 for (int j = 0; j < num_mlist(i); j++) { 00199 sbmat.set_new(i, mlist(i, j) - 1, bin(1)); 00200 } 00201 } 00202 sbmat.compact(); 00203 00204 if (transpose) { 00205 return sbmat.transpose(); 00206 } else { 00207 return sbmat; 00208 } 00209 } 00210 00211 00212 // ---------------------------------------------------------------------- 00213 // WARNING: This method is very slow. Sparse_Mat has to be extended with 00214 // some extra functions to improve the performance of this. 00215 // ---------------------------------------------------------------------- 00216 void GF2mat_sparse_alist::from_sparse(const GF2mat_sparse &sbmat, 00217 bool transpose) 00218 { 00219 if (transpose) { 00220 from_sparse(sbmat.transpose(), false); 00221 } 00222 else { 00223 // check matrix dimension 00224 M = sbmat.rows(); 00225 N = sbmat.cols(); 00226 00227 num_mlist.set_size(M); 00228 num_nlist.set_size(N); 00229 00230 // fill mlist matrix, num_mlist vector and max_num_m 00231 mlist.set_size(0, 0); 00232 for (int i = 0; i < M; i++) { 00233 ivec temp_row(0); 00234 for (int j = 0; j < N; j++) { 00235 if (sbmat(i, j) == bin(1)) { 00236 temp_row = concat(temp_row, j + 1); 00237 } 00238 } 00239 mlist.set_size(mlist.rows(), std::max(mlist.cols(), temp_row.size()), 00240 true); 00241 mlist.append_row(temp_row); 00242 num_mlist(i) = temp_row.length(); 00243 } 00244 max_num_m = max(num_mlist); 00245 00246 // fill nlist matrix, num_nlist vector and max_num_n 00247 nlist.set_size(0, 0); 00248 for (int j = 0; j < N; j++) { 00249 ivec temp_row(0); 00250 for (int i = 0; i < M; i++) { 00251 if (sbmat(i, j) == bin(1)) { 00252 temp_row = concat(temp_row, i + 1); 00253 } 00254 } 00255 nlist.set_size(nlist.rows(), std::max(nlist.cols(), temp_row.size()), 00256 true); 00257 nlist.append_row(temp_row); 00258 num_nlist(j) = temp_row.length(); 00259 } 00260 max_num_n = max(num_nlist); 00261 00262 data_ok = true; 00263 } 00264 } 00265 00266 00267 // ---------------------------------------------------------------------- 00268 // Implementation of a dense GF2 matrix class 00269 // ---------------------------------------------------------------------- 00270 00271 GF2mat::GF2mat(int i, int j): nrows(i), ncols(j), 00272 nwords((j >> shift_divisor) + 1) 00273 { 00274 data.set_size(nrows, nwords); 00275 data.clear(); 00276 } 00277 00278 GF2mat::GF2mat(): nrows(1), ncols(1), nwords(1) 00279 { 00280 data.set_size(nrows, nwords); 00281 data.clear(); 00282 } 00283 00284 GF2mat::GF2mat(const bvec &x, bool is_column) 00285 { 00286 if (is_column) { // create column vector 00287 nrows = length(x); 00288 ncols = 1; 00289 nwords = 1; 00290 data.set_size(nrows, nwords); 00291 data.clear(); 00292 for (int i = 0; i < nrows; i++) { 00293 set(i, 0, x(i)); 00294 } 00295 } else { // create row vector 00296 nrows = 1; 00297 ncols = length(x); 00298 nwords = (ncols >> shift_divisor) + 1; 00299 data.set_size(nrows, nwords); 00300 data.clear(); 00301 for (int i = 0; i < ncols; i++) { 00302 set(0, i, x(i)); 00303 } 00304 } 00305 } 00306 00307 00308 GF2mat::GF2mat(const bmat &X): nrows(X.rows()), ncols(X.cols()) 00309 { 00310 nwords = (ncols >> shift_divisor) + 1; 00311 data.set_size(nrows, nwords); 00312 data.clear(); 00313 for (int i = 0; i < nrows; i++) { 00314 for (int j = 0; j < ncols; j++) { 00315 set(i, j, X(i, j)); 00316 } 00317 } 00318 } 00319 00320 00321 GF2mat::GF2mat(const GF2mat_sparse &X) 00322 { 00323 nrows=X.rows(); 00324 ncols=X.cols(); 00325 nwords = (ncols >> shift_divisor) + 1; 00326 data.set_size(nrows,nwords); 00327 for (int i=0; i<nrows; i++) { 00328 for (int j=0; j<nwords; j++) { 00329 data(i,j) = 0; 00330 } 00331 } 00332 00333 for (int j=0; j<ncols; j++) { 00334 for (int i=0; i<X.get_col(j).nnz(); i++) { 00335 bin b = X.get_col(j).get_nz_data(i); 00336 set(X.get_col(j).get_nz_index(i),j,b); 00337 } 00338 } 00339 } 00340 00341 GF2mat::GF2mat(const GF2mat_sparse &X, int m1, int n1, int m2, int n2) 00342 { 00343 it_assert(X.rows()>m2,"GF2mat(): indexes out of range"); 00344 it_assert(X.cols()>n2,"GF2mat(): indexes out of range"); 00345 it_assert(m1>=0 && n1>=0 && m2>=m1 && n2>=n1, 00346 "GF2mat::GF2mat(): indexes out of range"); 00347 00348 nrows=m2-m1+1; 00349 ncols=n2-n1+1; 00350 nwords = (ncols >> shift_divisor) + 1; 00351 data.set_size(nrows,nwords); 00352 00353 for (int i=0; i<nrows; i++) { 00354 for (int j=0; j<nwords; j++) { 00355 data(i,j) = 0; 00356 } 00357 } 00358 00359 for (int i=0; i<nrows; i++) { 00360 for (int j=0; j<ncols; j++) { 00361 bin b = X(i+m1,j+n1); 00362 set(i,j,b); 00363 } 00364 } 00365 } 00366 00367 00368 GF2mat::GF2mat(const GF2mat_sparse &X, const ivec &columns) 00369 { 00370 it_assert(X.cols()>max(columns), 00371 "GF2mat::GF2mat(): index out of range"); 00372 it_assert(min(columns)>=0, 00373 "GF2mat::GF2mat(): column index must be positive"); 00374 00375 nrows=X.rows(); 00376 ncols=length(columns); 00377 nwords = (ncols >> shift_divisor) + 1; 00378 data.set_size(nrows,nwords); 00379 00380 for (int i=0; i<nrows; i++) { 00381 for (int j=0; j<nwords; j++) { 00382 data(i,j) = 0; 00383 } 00384 } 00385 00386 for (int j=0; j<ncols; j++) { 00387 for (int i=0; i<X.get_col(columns(j)).nnz(); i++) { 00388 bin b = X.get_col(columns(j)).get_nz_data(i); 00389 set(X.get_col(columns(j)).get_nz_index(i),j,b); 00390 } 00391 } 00392 } 00393 00394 00395 void GF2mat::set_size(int m, int n, bool copy) 00396 { 00397 nrows = m; 00398 ncols = n; 00399 nwords = (ncols >> shift_divisor) + 1; 00400 data.set_size(nrows, nwords, copy); 00401 if (!copy) 00402 data.clear(); 00403 } 00404 00405 00406 GF2mat_sparse GF2mat::sparsify() const 00407 { 00408 GF2mat_sparse Z(nrows,ncols); 00409 for (int i=0; i<nrows; i++) { 00410 for (int j=0; j<ncols; j++) { 00411 if (get(i,j)==1) { 00412 Z.set(i,j,1); 00413 } 00414 } 00415 } 00416 00417 return Z; 00418 } 00419 00420 bvec GF2mat::bvecify() const 00421 { 00422 it_assert(nrows==1 || ncols==1, 00423 "GF2mat::bvecify() matrix must be a vector"); 00424 int n = (nrows == 1 ? ncols : nrows); 00425 bvec result(n); 00426 if (nrows == 1) { 00427 for (int i = 0; i < n; i++) { 00428 result(i) = get(0, i); 00429 } 00430 } else { 00431 for (int i = 0; i < n; i++) { 00432 result(i) = get(i, 0); 00433 } 00434 } 00435 return result; 00436 } 00437 00438 00439 void GF2mat::set_row(int i, bvec x) 00440 { 00441 it_assert(length(x)==ncols, 00442 "GF2mat::set_row(): dimension mismatch"); 00443 for (int j=0; j<ncols; j++) { 00444 set(i,j,x(j)); 00445 } 00446 } 00447 00448 void GF2mat::set_col(int j, bvec x) 00449 { 00450 it_assert(length(x)==nrows, 00451 "GF2mat::set_col(): dimension mismatch"); 00452 for (int i=0; i<nrows; i++) { 00453 set(i,j,x(i)); 00454 } 00455 } 00456 00457 00458 GF2mat gf2dense_eye(int m) 00459 { 00460 GF2mat Z(m,m); 00461 for (int i=0; i<m; i++) { 00462 Z.set(i,i,1); 00463 } 00464 return Z; 00465 } 00466 00467 GF2mat GF2mat::get_submatrix(int m1, int n1, int m2, int n2) const 00468 { 00469 it_assert(m1>=0 && n1>=0 && m2>=m1 && n2>=n1 00470 && m2<nrows && n2<ncols, 00471 "GF2mat::get_submatrix() index out of range"); 00472 GF2mat result(m2-m1+1,n2-n1+1); 00473 00474 for (int i=m1; i<=m2; i++) { 00475 for (int j=n1; j<=n2; j++) { 00476 result.set(i-m1,j-n1,get(i,j)); 00477 } 00478 } 00479 00480 return result; 00481 } 00482 00483 00484 GF2mat GF2mat::concatenate_vertical(const GF2mat &X) const 00485 { 00486 it_assert(X.ncols==ncols, 00487 "GF2mat::concatenate_vertical(): dimension mismatch"); 00488 it_assert(X.nwords==nwords, 00489 "GF2mat::concatenate_vertical(): dimension mismatch"); 00490 00491 GF2mat result(nrows+X.nrows,ncols); 00492 for (int i=0; i<nrows; i++) { 00493 for (int j=0; j<nwords; j++) { 00494 result.data(i,j) = data(i,j); 00495 } 00496 } 00497 00498 for (int i=0; i<X.nrows; i++) { 00499 for (int j=0; j<nwords; j++) { 00500 result.data(i+nrows,j) = X.data(i,j); 00501 } 00502 } 00503 00504 return result; 00505 } 00506 00507 GF2mat GF2mat::concatenate_horizontal(const GF2mat &X) const 00508 { 00509 it_assert(X.nrows==nrows, 00510 "GF2mat::concatenate_horizontal(): dimension mismatch"); 00511 00512 GF2mat result(nrows,X.ncols+ncols); 00513 for (int i=0; i<nrows; i++) { 00514 for (int j=0; j<ncols; j++) { 00515 result.set(i,j,get(i,j)); 00516 } 00517 } 00518 00519 for (int i=0; i<nrows; i++) { 00520 for (int j=0; j<X.ncols; j++) { 00521 result.set(i,j+ncols,X.get(i,j)); 00522 } 00523 } 00524 00525 return result; 00526 } 00527 00528 bvec GF2mat::get_row(int i) const 00529 { 00530 bvec result(ncols); 00531 for (int j=0; j<ncols; j++) { 00532 result(j) = get(i,j); 00533 } 00534 00535 return result; 00536 } 00537 00538 bvec GF2mat::get_col(int j) const 00539 { 00540 bvec result(nrows); 00541 for (int i=0; i<nrows; i++) { 00542 result(i) = get(i,j); 00543 } 00544 00545 return result; 00546 } 00547 00548 00549 int GF2mat::T_fact(GF2mat &T, GF2mat &U, ivec &perm) const 00550 { 00551 T = gf2dense_eye(nrows); 00552 U = *this; 00553 00554 perm = zeros_i(ncols); 00555 for (int i=0; i<ncols; i++) { 00556 perm(i) = i; 00557 } 00558 00559 if (nrows>250) { // avoid cluttering output ... 00560 it_info_debug("Performing T-factorization of GF(2) matrix... rows: " 00561 << nrows << " cols: " << ncols << " .... " << std::endl); 00562 } 00563 int pdone=0; 00564 for (int j=0; j<nrows; j++) { 00565 // Now working on diagonal element j,j 00566 // First try find a row with a 1 in column i 00567 int i1,j1; 00568 for (i1=j; i1<nrows; i1++) { 00569 for (j1=j; j1<ncols; j1++) { 00570 if (U.get(i1,j1)==1) { goto found; } 00571 } 00572 } 00573 00574 return j; 00575 00576 found: 00577 U.swap_rows(i1,j); 00578 T.swap_rows(i1,j); 00579 U.swap_cols(j1,j); 00580 00581 int temp = perm(j); 00582 perm(j) = perm(j1); 00583 perm(j1) = temp; 00584 00585 // now subtract row i from remaining rows 00586 for (int i1=j+1; i1<nrows; i1++) { 00587 if (U.get(i1,j)==1) { 00588 int ptemp = floor_i(100.0*(i1+j*nrows)/(nrows*nrows)); 00589 if (nrows>250 && ptemp>pdone+10) { 00590 it_info_debug(ptemp << "% done."); 00591 pdone=ptemp; 00592 } 00593 U.add_rows(i1,j); 00594 T.add_rows(i1,j); 00595 } 00596 } 00597 } 00598 return nrows; 00599 } 00600 00601 00602 int GF2mat::T_fact_update_bitflip(GF2mat &T, GF2mat &U, 00603 ivec &perm, int rank, 00604 int r, int c) const 00605 { 00606 // First, update U (before re-triangulization) 00607 int c0=c; 00608 for (c=0; c<ncols; c++) { 00609 if (perm(c)==c0) { 00610 goto foundidx; 00611 } 00612 } 00613 it_error("GF2mat::T_fact_update_bitflip() - internal error"); 00614 00615 foundidx: 00616 for (int i=0; i<nrows; i++) { 00617 if (T.get(i,r)==1) { 00618 U.addto_element(i,c,1); 00619 } 00620 } 00621 00622 // first move column c to the end 00623 bvec lastcol = U.get_col(c); 00624 int temp_perm = perm(c); 00625 for (int j=c; j<ncols-1; j++) { 00626 U.set_col(j,U.get_col(j+1)); 00627 perm(j) = perm(j+1); 00628 } 00629 U.set_col(ncols-1,lastcol); 00630 perm(ncols-1) = temp_perm; 00631 00632 // then, if the matrix is tall, also move row c to the end 00633 if (nrows>=ncols) { 00634 bvec lastrow_U = U.get_row(c); 00635 bvec lastrow_T = T.get_row(c); 00636 for (int i=c; i<nrows-1; i++) { 00637 U.set_row(i,U.get_row(i+1)); 00638 T.set_row(i,T.get_row(i+1)); 00639 } 00640 U.set_row(nrows-1,lastrow_U); 00641 T.set_row(nrows-1,lastrow_T); 00642 00643 // Do Gaussian elimination on the last row 00644 for (int j=c; j<ncols; j++) { 00645 if (U.get(nrows-1,j)==1) { 00646 U.add_rows(nrows-1,j); 00647 T.add_rows(nrows-1,j); 00648 } 00649 } 00650 } 00651 00652 // Now, continue T-factorization from the point (rank-1,rank-1) 00653 for (int j=rank-1; j<nrows; j++) { 00654 int i1,j1; 00655 for (i1=j; i1<nrows; i1++) { 00656 for (j1=j; j1<ncols; j1++) { 00657 if (U.get(i1,j1)==1) { 00658 goto found; 00659 } 00660 } 00661 } 00662 00663 return j; 00664 00665 found: 00666 U.swap_rows(i1,j); 00667 T.swap_rows(i1,j); 00668 U.swap_cols(j1,j); 00669 00670 int temp = perm(j); 00671 perm(j) = perm(j1); 00672 perm(j1) = temp; 00673 00674 for (int i1=j+1; i1<nrows; i1++) { 00675 if (U.get(i1,j)==1) { 00676 U.add_rows(i1,j); 00677 T.add_rows(i1,j); 00678 } 00679 } 00680 } 00681 00682 return nrows; 00683 } 00684 00685 bool GF2mat::T_fact_update_addcol(GF2mat &T, GF2mat &U, 00686 ivec &perm, bvec newcol) const 00687 { 00688 int i0 = T.rows(); 00689 int j0 = U.cols(); 00690 it_assert(j0>0,"GF2mat::T_fact_update_addcol(): dimension mismatch"); 00691 it_assert(i0==T.cols(), 00692 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00693 it_assert(i0==U.rows(), 00694 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00695 it_assert(length(perm)==j0, 00696 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00697 it_assert(U.get(j0-1,j0-1)==1, 00698 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00699 // The following test is VERY TIME-CONSUMING 00700 it_assert_debug(U.row_rank()==j0, 00701 "GF2mat::T_fact_update_addcol(): factorization has incorrect rank"); 00702 00703 bvec z = T*newcol; 00704 GF2mat Utemp = U.concatenate_horizontal(GF2mat(z,1)); 00705 00706 // start working on position (j0,j0) 00707 int i; 00708 for (i=j0; i<i0; i++) { 00709 if (Utemp.get(i,j0)==1) { 00710 goto found; 00711 } 00712 } 00713 return (false); // adding the new column would not improve the rank 00714 00715 found: 00716 perm.set_length(j0+1,true); 00717 perm(j0) = j0; 00718 00719 Utemp.swap_rows(i,j0); 00720 T.swap_rows(i,j0); 00721 00722 for (int i1=j0+1; i1<i0; i1++) { 00723 if (Utemp.get(i1,j0)==1) { 00724 Utemp.add_rows(i1,j0); 00725 T.add_rows(i1,j0); 00726 } 00727 } 00728 00729 U = Utemp; 00730 return (true); // the new column was successfully added 00731 } 00732 00733 00734 00735 00736 GF2mat GF2mat::inverse() const 00737 { 00738 it_assert(nrows==ncols,"GF2mat::inverse(): Matrix must be square"); 00739 00740 // first compute the T-factorization 00741 GF2mat T,U; 00742 ivec perm; 00743 int rank = T_fact(T,U,perm); 00744 it_assert(rank==ncols,"GF2mat::inverse(): Matrix is not full rank"); 00745 00746 // backward substitution 00747 for (int i=ncols-2; i>=0; i--) { 00748 for (int j=ncols-1; j>i; j--) { 00749 if (U.get(i,j)==1) { 00750 U.add_rows(i,j); 00751 T.add_rows(i,j); 00752 } 00753 } 00754 } 00755 T.permute_rows(perm,1); 00756 return T; 00757 } 00758 00759 int GF2mat::row_rank() const 00760 { 00761 GF2mat T,U; 00762 ivec perm; 00763 int rank = T_fact(T,U,perm); 00764 return rank; 00765 } 00766 00767 bool GF2mat::is_zero() const 00768 { 00769 for (int i=0; i<nrows; i++) { 00770 for (int j=0; j<nwords; j++) { 00771 if (data(i,j)!=0) { 00772 return false; 00773 } 00774 } 00775 } 00776 return true; 00777 } 00778 00779 bool GF2mat::operator==(const GF2mat &X) const 00780 { 00781 if (X.nrows!=nrows) { return false; } 00782 if (X.ncols!=ncols) { return false; } 00783 it_assert(X.nwords==nwords,"GF2mat::operator==() dimension mismatch"); 00784 00785 for (int i=0; i<nrows; i++) { 00786 for (int j=0; j<nwords; j++) { 00787 // if (X.get(i,j)!=get(i,j)) { 00788 if (X.data(i,j)!=data(i,j)) { 00789 return false; 00790 } 00791 } 00792 } 00793 return true; 00794 } 00795 00796 00797 void GF2mat::add_rows(int i, int j) 00798 { 00799 it_assert(i>=0 && i<nrows,"GF2mat::add_rows(): out of range"); 00800 it_assert(j>=0 && j<nrows,"GF2mat::add_rows(): out of range"); 00801 for (int k=0; k<nwords; k++) { 00802 data(i,k) ^= data(j,k); 00803 } 00804 } 00805 00806 void GF2mat::swap_rows(int i, int j) 00807 { 00808 it_assert(i>=0 && i<nrows,"GF2mat::swap_rows(): index out of range"); 00809 it_assert(j>=0 && j<nrows,"GF2mat::swap_rows(): index out of range"); 00810 for (int k=0; k<nwords; k++) { 00811 int temp = data(i,k); 00812 data(i,k) = data(j,k); 00813 data(j,k) = temp; 00814 } 00815 } 00816 00817 void GF2mat::swap_cols(int i, int j) 00818 { 00819 it_assert(i>=0 && i<ncols,"GF2mat::swap_cols(): index out of range"); 00820 it_assert(j>=0 && j<ncols,"GF2mat::swap_cols(): index out of range"); 00821 bvec temp = get_col(i); 00822 set_col(i,get_col(j)); 00823 set_col(j,temp); 00824 } 00825 00826 00827 void GF2mat::operator=(const GF2mat &X) 00828 { 00829 nrows=X.nrows; 00830 ncols=X.ncols; 00831 nwords=X.nwords; 00832 data = X.data; 00833 } 00834 00835 GF2mat operator*(const GF2mat &X, const GF2mat &Y) 00836 { 00837 it_assert(X.ncols==Y.nrows,"GF2mat::operator*(): dimension mismatch"); 00838 it_assert(X.nwords>0,"Gfmat::operator*(): dimension mismatch"); 00839 it_assert(Y.nwords>0,"Gfmat::operator*(): dimension mismatch"); 00840 00841 /* 00842 // this can be done more efficiently? 00843 GF2mat result(X.nrows,Y.ncols); 00844 for (int i=0; i<X.nrows; i++) { 00845 for (int j=0; j<Y.ncols; j++) { 00846 bin b=0; 00847 for (int k=0; k<X.ncols; k++) { 00848 bin x = X.get(i,k); 00849 bin y = Y.get(k,j); 00850 b ^= (x&y); 00851 } 00852 result.set(i,j,b); 00853 } 00854 } 00855 return result; */ 00856 00857 // is this better? 00858 return mult_trans(X,Y.transpose()); 00859 } 00860 00861 bvec operator*(const GF2mat &X, const bvec &y) 00862 { 00863 it_assert(length(y)==X.ncols,"GF2mat::operator*(): dimension mismatch"); 00864 it_assert(X.nwords>0,"Gfmat::operator*(): dimension mismatch"); 00865 00866 /* 00867 // this can be done more efficiently? 00868 bvec result = zeros_b(X.nrows); 00869 for (int i=0; i<X.nrows; i++) { 00870 // do the nwords-1 data columns first 00871 for (int j=0; j<X.nwords-1; j++) { 00872 int ind = j<<shift_divisor; 00873 unsigned char r=X.data(i,j); 00874 while (r) { 00875 result(i) ^= (r & y(ind)); 00876 r >>= 1; 00877 ind++; 00878 } 00879 } 00880 // do the last column separately 00881 for (int j=(X.nwords-1)<<shift_divisor; j<X.ncols; j++) { 00882 result(i) ^= (X.get(i,j) & y(j)); 00883 } 00884 } 00885 return result; */ 00886 00887 // is this better? 00888 return (mult_trans(X,GF2mat(y,0))).bvecify(); 00889 } 00890 00891 GF2mat mult_trans(const GF2mat &X, const GF2mat &Y) 00892 { 00893 it_assert(X.ncols==Y.ncols,"GF2mat::mult_trans(): dimension mismatch"); 00894 it_assert(X.nwords>0,"GF2mat::mult_trans(): dimension mismatch"); 00895 it_assert(Y.nwords>0,"GF2mat::mult_trans(): dimension mismatch"); 00896 it_assert(X.nwords==Y.nwords,"GF2mat::mult_trans(): dimension mismatch"); 00897 00898 GF2mat result(X.nrows,Y.nrows); 00899 00900 for (int i=0; i<X.nrows; i++) { 00901 for (int j=0; j<Y.nrows; j++) { 00902 bin b=0; 00903 int kindx =i; 00904 int kindy =j; 00905 for (int k=0; k<X.nwords; k++) { 00906 unsigned char r = X.data(kindx) & Y.data(kindy); 00907 /* The following can be speeded up by using a small lookup 00908 table for the number of ones and shift only a few times (or 00909 not at all if table is large) */ 00910 while (r) { 00911 b ^= r&1; 00912 r>>=1; 00913 }; 00914 kindx += X.nrows; 00915 kindy += Y.nrows; 00916 } 00917 result.set(i,j,b); 00918 } 00919 } 00920 return result; 00921 } 00922 00923 GF2mat GF2mat::transpose() const 00924 { 00925 // CAN BE SPEEDED UP 00926 GF2mat result(ncols,nrows); 00927 00928 for (int i=0; i<nrows; i++) { 00929 for (int j=0; j<ncols; j++) { 00930 result.set(j,i,get(i,j)); 00931 } 00932 } 00933 return result; 00934 } 00935 00936 GF2mat operator+(const GF2mat &X, const GF2mat &Y) 00937 { 00938 it_assert(X.nrows==Y.nrows,"GF2mat::operator+(): dimension mismatch"); 00939 it_assert(X.ncols==Y.ncols,"GF2mat::operator+(): dimension mismatch"); 00940 it_assert(X.nwords==Y.nwords,"GF2mat::operator+(): dimension mismatch"); 00941 GF2mat result(X.nrows,X.ncols); 00942 00943 for (int i=0; i<X.nrows; i++) { 00944 for (int j=0; j<X.nwords; j++) { 00945 result.data(i,j) = X.data(i,j) ^ Y.data(i,j); 00946 } 00947 } 00948 00949 return result; 00950 } 00951 00952 void GF2mat::permute_cols(ivec &perm, bool I) 00953 { 00954 it_assert(length(perm)==ncols, 00955 "GF2mat::permute_cols(): dimensions do not match"); 00956 00957 GF2mat temp = (*this); 00958 for (int j=0; j<ncols; j++) { 00959 if (I==0) { 00960 set_col(j,temp.get_col(perm(j))); 00961 } else { 00962 set_col(perm(j),temp.get_col(j)); 00963 } 00964 } 00965 } 00966 00967 void GF2mat::permute_rows(ivec &perm, bool I) 00968 { 00969 it_assert(length(perm)==nrows, 00970 "GF2mat::permute_rows(): dimensions do not match"); 00971 00972 GF2mat temp = (*this); 00973 for (int i=0; i<nrows; i++) { 00974 if (I==0) { 00975 for (int j=0; j<nwords; j++) { 00976 data(i,j) = temp.data(perm(i),j); 00977 } 00978 } else { 00979 for (int j=0; j<nwords; j++) { 00980 data(perm(i),j) = temp.data(i,j); 00981 } 00982 } 00983 } 00984 } 00985 00986 00987 std::ostream &operator<<(std::ostream &os, const GF2mat &X) 00988 { 00989 int i,j; 00990 os << "---- GF(2) matrix of dimension " << X.nrows << "*" << X.ncols 00991 << " -- Density: " << X.density() << " ----" << std::endl; 00992 00993 for (i=0; i<X.nrows; i++) { 00994 os << " "; 00995 for (j=0; j<X.ncols; j++) { 00996 os << X.get(i,j) << " "; 00997 } 00998 os << std::endl; 00999 } 01000 01001 return os; 01002 } 01003 01004 double GF2mat::density() const 01005 { 01006 int no_of_ones=0; 01007 01008 for (int i=0; i<nrows; i++) { 01009 for (int j=0; j<ncols; j++) { 01010 no_of_ones += (get(i,j)==1 ? 1 : 0); 01011 } 01012 } 01013 01014 return ((double) no_of_ones)/(nrows*ncols); 01015 } 01016 01017 01018 it_file &operator<<(it_file &f, const GF2mat &X) 01019 { 01020 // 3 64-bit unsigned words for: nrows, ncols and nwords + rest for char data 01021 uint64_t bytecount = 3 * sizeof(uint64_t) 01022 + X.nrows * X.nwords * sizeof(char); 01023 f.write_data_header("GF2mat", bytecount); 01024 01025 f.low_level_write(static_cast<uint64_t>(X.nrows)); 01026 f.low_level_write(static_cast<uint64_t>(X.ncols)); 01027 f.low_level_write(static_cast<uint64_t>(X.nwords)); 01028 for (int i=0; i<X.nrows; i++) { 01029 for (int j=0; j<X.nwords; j++) { 01030 f.low_level_write(static_cast<char>(X.data(i,j))); 01031 } 01032 } 01033 return f; 01034 } 01035 01036 it_ifile &operator>>(it_ifile &f, GF2mat &X) 01037 { 01038 it_file::data_header h; 01039 01040 f.read_data_header(h); 01041 if (h.type == "GF2mat") { 01042 uint64_t tmp; 01043 f.low_level_read(tmp); X.nrows = static_cast<int>(tmp); 01044 f.low_level_read(tmp); X.ncols = static_cast<int>(tmp); 01045 f.low_level_read(tmp); X.nwords = static_cast<int>(tmp); 01046 X.data.set_size(X.nrows,X.nwords); 01047 for (int i=0; i<X.nrows; i++) { 01048 for (int j=0; j<X.nwords; j++) { 01049 char r; 01050 f.low_level_read(r); 01051 X.data(i,j) = static_cast<unsigned char>(r); 01052 } 01053 } 01054 } 01055 else { 01056 it_error("it_ifile &operator>>() - internal error"); 01057 } 01058 01059 return f; 01060 } 01061 01062 } // namespace itpp 01063
Generated on Sun Dec 9 17:30:59 2007 for IT++ by Doxygen 1.5.4